function [ a, ind ] = elmhes ( nm, n, low, igh, a )

%*****************************************************************************80
%
%% ELMHES transforms a real general matrix to upper Hessenberg form.
%
%  Discussion:
%
%    Given a real general matrix, this routine reduces a submatrix
%    situated in rows and columns LOW through IGH to upper Hessenberg
%    form by stabilized elementary similarity transformations.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 October 2007
%
%  Author:
%
%    This version by John Burkardt.
%
%  Reference:
%
%    RS Martin, James Wilkinson,
%    ELMHES,
%    Numerische Mathematik,
%    Volume 12, pages 349-368, 1968.
%
%    James Wilkinson, Christian Reinsch,
%    Handbook for Automatic Computation,
%    Volume II, Linear Algebra, Part 2,
%    Springer Verlag, 1971.
%
%    Brian Smith, James Boyle, Jack Dongarra, Burton Garbow, Y Ikebe, 
%    Virginia Klema, Cleve Moler,
%    Matrix Eigensystem Routines, EISPACK Guide,
%    Lecture Notes in Computer Science, Volume 6,
%    Springer Verlag, 1976.
%
%  Input:
%
%    integer NM, the leading dimension of the array A.
%    NM must be at least N.
%
%    integer N, the order of the matrix.
%
%    integer LOW, IGH, are determined by the balancing 
%    routine BALANC.  If BALANC has not been used, set LOW = 1, IGH = N.
%
%    real A(NM,N), the matrix to be reduced.
%
%  Output:
%
%    real A(NM,N), the Hessenberg matrix.  The multipliers
%    which were used in the reduction are stored in the
%    remaining triangle under the Hessenberg matrix.
%
%    integer IND(N), contains information on the rows 
%    and columns interchanged in the reduction.  Only elements LOW through 
%    IGH are used.
%
  la = igh - 1;

  ind(1:n) = 0;

  for m = low + 1 : la

    x = 0.0;
    i = m;

    for j = m : igh
      if ( abs ( x ) < abs ( a(j,m-1) ) )
        x = a(j,m-1);
        i = j;
      end
    end

    ind(m) = i;
%
%  Interchange rows and columns of the matrix.
%
    if ( i ~= m )

      b(1,m-1:n) = a(i,m-1:n);
      a(i,m-1:n) = a(m,m-1:n);
      a(m,m-1:n) = b(1,m-1:n);

      c(1:igh,1) = a(1:igh,i);
      a(1:igh,i) = a(1:igh,m);
      a(1:igh,m) = c(1:igh,1);

    end

    if ( x ~= 0.0 )

      ffor i = m+1 : igh

        y = a(i,m-1)

        if ( y ~= 0.0 )

          y = y / x;

          a(i,m-1) = y;
          a(i,m:n) = a(i,m:n) - y * a(m,m:n)

          a(1:igh,m) = a(1:igh,m) + y * a(1:igh,i);

        end

      end

    end

  end

  return
end

