Reduce a real general matrix to upper Hessenberg form using stabilized elementary similarity transformations.
This subroutine is a translation of the ALGOL procedure ELMHES, NUM. MATH. 12, 349-368(1968) by Martin and Wilkinson. HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
Given a REAL GENERAL matrix, this subroutine reduces a submatrix situated in rows and columns LOW through IGH to upper Hessenberg form by stabilized elementary similarity transformations.
NM must be set to the row dimension of the two-dimensional
array parameter, A, as declared in the calling program
dimension statement. NM is an INTEGER variable.
N is the order of the matrix, A. N is an INTEGER variable.
N must be less than or equal to NM.
LOW and IGH are two INTEGER variables determined by the
balancing subroutine BALANC. If BALANC has not been
used, set LOW=1 and IGH equal to the order of the matrix, N.
A contains the input matrix. A is a two-dimensional REAL
array, dimensioned A(NM,N).
A contains the upper Hessenberg matrix. The multipliers which
were used in the reduction are stored in the remaining
triangle under the Hessenberg matrix.
INTV contains information on the rows and columns interchanged
in the reduction. Only elements LOW through IGH are used.
INTV is a one-dimensional INTEGER array, dimensioned INTV(IGH).
Questions and comments should be directed to B. S. Garbow,
APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer | :: | Nm | ||||
integer | :: | n | ||||
integer | :: | Low | ||||
integer | :: | Igh | ||||
real(kind=wp) | :: | a(Nm,*) | ||||
integer | :: | Intv(*) |
subroutine elmhes(Nm, n, Low, Igh, a, Intv) implicit none integer :: i, j, m, n, la, Nm, Igh, kp1, Low, mm1, mp1 real(wp) :: a(Nm, *) real(wp) :: x, y integer :: Intv(*) la = Igh - 1 kp1 = Low + 1 if (la >= kp1) then do m = kp1, la mm1 = m - 1 x = 0.0_wp i = m do j = m, Igh if (abs(a(j, mm1)) > abs(x)) then x = a(j, mm1) i = j endif enddo Intv(m) = i if (i /= m) then ! INTERCHANGE ROWS AND COLUMNS OF A do j = mm1, n y = a(i, j) a(i, j) = a(m, j) a(m, j) = y enddo do j = 1, Igh y = a(j, i) a(j, i) = a(j, m) a(j, m) = y enddo endif ! END INTERCHANGE if (x /= 0.0_wp) then mp1 = m + 1 do i = mp1, Igh y = a(i, mm1) if (y /= 0.0_wp) then y = y/x a(i, mm1) = y do j = m, n a(i, j) = a(i, j) - y*a(m, j) enddo do j = 1, Igh a(j, m) = a(j, m) + y*a(j, i) enddo endif enddo endif enddo endif end subroutine elmhes