Compute the eigenvalues of a real upper Hessenberg matrix using the QR method.
This subroutine is a translation of the ALGOL procedure HQR, NUM. MATH. 14, 219-231(1970) by Martin, Peters, and Wilkinson. HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971).
This subroutine finds the eigenvalues of a REAL UPPER Hessenberg matrix by the QR method.
NM must be set to the row dimension of the two-dimensional
array parameter, H, as declared in the calling program
dimension statement. NM is an INTEGER variable.
N is the order of the matrix H. 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.
H contains the upper Hessenberg matrix. Information about
the transformations used in the reduction to Hessenberg
form by ELMHES or ORTHES, if performed, is stored
in the remaining triangle under the Hessenberg matrix.
H is a two-dimensional REAL array, dimensioned H(NM,N).
H has been destroyed. Therefore, it must be saved before
calling HQR if subsequent calculation and back
transformation of eigenvectors is to be performed.
WR and WI contain the real and imaginary parts, respectively,
of the eigenvalues. The eigenvalues are unordered except
that complex conjugate pairs of values appear consecutively
with the eigenvalue having the positive imaginary part first.
If an error exit is made, the eigenvalues should be correct
for indices IERR+1, IERR+2, ..., N. WR and WI are one-
dimensional REAL arrays, dimensioned WR(N) and WI(N).
IERR is an INTEGER flag set to
Zero for normal return,
J if the J-th eigenvalue has not been
determined after a total of 30*N iterations.
The eigenvalues should be correct for indices
IERR+1, IERR+2, ..., N.
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) | :: | h(Nm,*) | ||||
| real(kind=wp) | :: | Wr(*) | ||||
| real(kind=wp) | :: | Wi(*) | ||||
| integer | :: | Ierr |
subroutine hqr(Nm, n, Low, Igh, h, Wr, Wi, Ierr) implicit none integer :: i, j, k, l, m, n, en, ll, mm, na, Nm, Igh, & itn, its, Low, mp2, enm2, Ierr, gt real(wp) :: h(Nm, *), Wr(*), Wi(*) real(wp) :: p, q, r, s, t, w, x, y, zz, norm, s1, s2 logical :: notlas gt = 0 Ierr = 0 norm = 0.0_wp k = 1 ! STORE ROOTS ISOLATED BY BALANC ! AND COMPUTE MATRIX NORM do i = 1, n do j = k, n norm = norm + abs(h(i, j)) enddo k = i if (i < Low .or. i > Igh) then Wr(i) = h(i, i) Wi(i) = 0.0_wp endif enddo en = Igh t = 0.0_wp itn = 30*n ! SEARCH FOR NEXT EIGENVALUES do while (en >= Low) gt = 0 its = 0 na = en - 1 enm2 = na - 1 ! LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT ! FOR L=EN STEP -1 UNTIL LOW DO -- do while (.true.) do ll = Low, en l = en + Low - ll if (l == Low) exit s = abs(h(l - 1, l - 1)) + abs(h(l, l)) if (s == 0.0_wp) s = norm s2 = s + abs(h(l, l - 1)) if (s2 == s) exit enddo ! FORM SHIFT x = h(en, en) if (l == en) then ! ONE ROOT FOUND Wr(en) = x + t Wi(en) = 0.0_wp en = na gt = 1 exit else y = h(na, na) w = h(en, na)*h(na, en) if (l == na) exit if (itn == 0) then ! SET ERROR -- NO CONVERGENCE TO AN ! EIGENVALUE AFTER 30*N ITERATIONS Ierr = en return else if (its == 10 .or. its == 20) then ! FORM EXCEPTIONAL SHIFT t = t + x do i = Low, en h(i, i) = h(i, i) - x enddo s = abs(h(en, na)) + abs(h(na, enm2)) x = 0.75_wp*s y = x w = -0.4375_wp*s*s endif its = its + 1 itn = itn - 1 ! LOOK FOR TWO CONSECUTIVE SMALL ! SUB-DIAGONAL ELEMENTS. ! FOR M=EN-2 STEP -1 UNTIL L DO -- do mm = l, enm2 m = enm2 + l - mm zz = h(m, m) r = x - zz s = y - zz p = (r*s - w)/h(m + 1, m) + h(m, m + 1) q = h(m + 1, m + 1) - zz - r - s r = h(m + 2, m + 1) s = abs(p) + abs(q) + abs(r) p = p/s q = q/s r = r/s if (m == l) exit s1 = abs(p)*(abs(h(m - 1, m - 1)) + abs(zz) + abs(h(m + 1, m + 1))) s2 = s1 + abs(h(m, m - 1))*(abs(q) + abs(r)) if (s2 == s1) exit enddo mp2 = m + 2 do i = mp2, en h(i, i - 2) = 0.0_wp if (i /= mp2) h(i, i - 3) = 0.0_wp enddo ! DOUBLE QR STEP INVOLVING ROWS L TO EN AND ! COLUMNS M TO EN do k = m, na notlas = k /= na if (k /= m) then p = h(k, k - 1) q = h(k + 1, k - 1) r = 0.0_wp if (notlas) r = h(k + 2, k - 1) x = abs(p) + abs(q) + abs(r) if (x == 0.0_wp) cycle p = p/x q = q/x r = r/x endif s = sign(sqrt(p*p + q*q + r*r), p) if (k == m) then if (l /= m) h(k, k - 1) = -h(k, k - 1) else h(k, k - 1) = -s*x endif p = p + s x = p/s y = q/s zz = r/s q = q/p r = r/p ! ROW MODIFICATION do j = k, en p = h(k, j) + q*h(k + 1, j) if (notlas) then p = p + r*h(k + 2, j) h(k + 2, j) = h(k + 2, j) - p*zz endif h(k + 1, j) = h(k + 1, j) - p*y h(k, j) = h(k, j) - p*x enddo j = min(en, k + 3) ! COLUMN MODIFICATION do i = l, j p = x*h(i, k) + y*h(i, k + 1) if (notlas) then p = p + zz*h(i, k + 2) h(i, k + 2) = h(i, k + 2) - p*r endif h(i, k + 1) = h(i, k + 1) - p*q h(i, k) = h(i, k) - p enddo enddo endif endif enddo ! TWO ROOTS FOUND if (gt == 0) then p = (y - x)/2.0_wp q = p*p + w zz = sqrt(abs(q)) x = x + t if (q < 0.0_wp) then ! COMPLEX PAIR Wr(na) = x + p Wr(en) = x + p Wi(na) = zz Wi(en) = -zz else ! REAL PAIR zz = p + sign(zz, p) Wr(na) = x + zz Wr(en) = Wr(na) if (zz /= 0.0_wp) Wr(en) = x - w/zz Wi(na) = 0.0_wp Wi(en) = 0.0_wp endif en = enm2 end if enddo end subroutine hqr