Compute the eigenvalues and eigenvectors of a real upper Hessenberg matrix using QR method.
This subroutine is a translation of the ALGOL procedure HQR2, NUM. MATH. 16, 181-204(1970) by Peters and Wilkinson. HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
This subroutine finds the eigenvalues and eigenvectors of a REAL UPPER Hessenberg matrix by the QR method. The eigenvectors of a REAL GENERAL matrix can also be found if ELMHES and ELTRAN or ORTHES and ORTRAN have been used to reduce this general matrix to Hessenberg form and to accumulate the similarity transformations.
NM must be set to the row dimension of the two-dimensional
array parameters, H and Z, 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. H is a two-dimensional
REAL array, dimensioned H(NM,N).
Z contains the transformation matrix produced by ELTRAN
after the reduction by ELMHES, or by ORTRAN after the
reduction by ORTHES, if performed. If the eigenvectors
of the Hessenberg matrix are desired, Z must contain the
identity matrix. Z is a two-dimensional REAL array,
dimensioned Z(NM,M).
H has been destroyed.
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).
Z contains the real and imaginary parts of the eigenvectors.
If the J-th eigenvalue is real, the J-th column of Z
contains its eigenvector. If the J-th eigenvalue is complex
with positive imaginary part, the J-th and (J+1)-th
columns of Z contain the real and imaginary parts of its
eigenvector. The eigenvectors are unnormalized. If an
error exit is made, none of the eigenvectors has been found.
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, but no eigenvectors are
computed.
Calls CDIV for complex division.
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(*) | ||||
| real(kind=wp) | :: | z(Nm,*) | ||||
| integer | :: | Ierr |
subroutine hqr2(Nm, n, Low, Igh, h, Wr, Wi, z, Ierr) implicit none integer :: i, j, k, l, m, n, en, ii, jj, ll, mm, na, Nm, & nn, gt integer :: Igh, itn, its, Low, mp2, enm2, Ierr real(wp) :: h(Nm, *), Wr(*), Wi(*), z(Nm, *) real(wp) :: p, q, r, s, t, w, x, y, ra, sa, vi, vr, 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 h(en, en) = x + t Wr(en) = h(en, en) 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, n 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 = 1, 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 ! ACCUMULATE TRANSFORMATIONS do i = Low, Igh ! p = x*z(i, k) + y*z(i, k + 1) if (notlas) then p = p + zz*z(i, k + 2) z(i, k + 2) = z(i, k + 2) - p*r endif z(i, k + 1) = z(i, k + 1) - p*q z(i, k) = z(i, k) - p ! enddo enddo endif endif enddo if (gt == 1) cycle ! TWO ROOTS FOUND p = (y - x)/2.0_wp q = p*p + w zz = sqrt(abs(q)) h(en, en) = x + t x = h(en, en) h(na, na) = y + 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 x = h(en, na) s = abs(x) + abs(zz) p = x/s q = zz/s r = sqrt(p*p + q*q) p = p/r q = q/r ! ROW MODIFICATION do j = na, n zz = h(na, j) h(na, j) = q*zz + p*h(en, j) h(en, j) = q*h(en, j) - p*zz enddo ! COLUMN MODIFICATION do i = 1, en zz = h(i, na) h(i, na) = q*zz + p*h(i, en) h(i, en) = q*h(i, en) - p*zz enddo ! ACCUMULATE TRANSFORMATIONS do i = Low, Igh zz = z(i, na) z(i, na) = q*zz + p*z(i, en) z(i, en) = q*z(i, en) - p*zz enddo endif en = enm2 enddo ! ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND ! VECTORS OF UPPER TRIANGULAR FORM if (norm /= 0.0_wp) then ! FOR EN=N STEP -1 UNTIL 1 DO -- do nn = 1, n en = n + 1 - nn p = Wr(en) q = Wi(en) na = en - 1 if (q < 0) then ! END COMPLEX VECTOR ! COMPLEX VECTOR m = na ! LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT ! EIGENVECTOR MATRIX IS TRIANGULAR if (abs(h(en, na)) <= abs(h(na, en))) then call cdiv(0.0_wp, -h(na, en), h(na, na) - p, q, h(na, na), h(na, en)) else h(na, na) = q/h(en, na) h(na, en) = -(h(en, en) - p)/h(en, na) endif h(en, na) = 0.0_wp h(en, en) = 1.0_wp enm2 = na - 1 if (enm2 /= 0) then ! FOR I=EN-2 STEP -1 UNTIL 1 DO -- do ii = 1, enm2 i = na - ii w = h(i, i) - p ra = 0.0_wp sa = h(i, en) ! do j = m, na ra = ra + h(i, j)*h(j, na) sa = sa + h(i, j)*h(j, en) enddo ! if (Wi(i) >= 0.0_wp) then m = i if (Wi(i) == 0.0_wp) then call cdiv(-ra, -sa, w, q, h(i, na), h(i, en)) else ! SOLVE COMPLEX EQUATIONS x = h(i, i + 1) y = h(i + 1, i) vr = (Wr(i) - p)*(Wr(i) - p) + Wi(i)*Wi(i) - q*q vi = (Wr(i) - p)*2.0_wp*q if (vr == 0.0_wp .and. vi == 0.0_wp) then s1 = norm*(abs(w) + abs(q) + abs(x) + abs(y) + abs(zz)) vr = s1 do while (.true.) vr = 0.5_wp*vr if (s1 + vr <= s1) exit enddo vr = 2.0_wp*vr endif call cdiv(x*r - zz*ra + q*sa, x*s - zz*sa - q*ra, vr, vi, h(i, na), & h(i, en)) if (abs(x) <= abs(zz) + abs(q)) then call cdiv(-r - y*h(i, na), -s - y*h(i, en), zz, q, h(i + 1, na), & h(i + 1, en)) else h(i + 1, na) = (-ra - w*h(i, na) + q*h(i, en))/x h(i + 1, en) = (-sa - w*h(i, en) - q*h(i, na))/x endif endif else zz = w r = ra s = sa endif enddo endif elseif (q == 0) then ! REAL VECTOR m = en h(en, en) = 1.0_wp if (na /= 0) then ! FOR I=EN-1 STEP -1 UNTIL 1 DO -- do ii = 1, na i = en - ii w = h(i, i) - p r = h(i, en) if (m <= na) then do j = m, na r = r + h(i, j)*h(j, en) enddo endif if (Wi(i) >= 0.0_wp) then ! END REAL VECTOR m = i if (Wi(i) == 0.0_wp) then t = w if (t == 0.0_wp) then t = norm do while (.true.) t = 0.5_wp*t if (norm + t <= norm) exit enddo t = 2.0_wp*t endif h(i, en) = -r/t else ! SOLVE REAL EQUATIONS x = h(i, i + 1) y = h(i + 1, i) q = (Wr(i) - p)*(Wr(i) - p) + Wi(i)*Wi(i) t = (x*s - zz*r)/q h(i, en) = t if (abs(x) <= abs(zz)) then h(i + 1, en) = (-s - y*t)/zz else h(i + 1, en) = (-r - w*t)/x endif endif else zz = w s = r endif enddo endif endif enddo ! END BACK SUBSTITUTION. ! VECTORS OF ISOLATED ROOTS do i = 1, n if (i < Low .or. i > Igh) then do j = i, n z(i, j) = h(i, j) enddo endif enddo ! MULTIPLY BY TRANSFORMATION MATRIX TO GIVE ! VECTORS OF ORIGINAL FULL MATRIX. ! FOR J=N STEP -1 UNTIL LOW DO -- do jj = Low, n j = n + Low - jj m = min(j, Igh) do i = Low, Igh zz = 0.0_wp do k = Low, m zz = zz + z(i, k)*h(k, j) enddo z(i, j) = zz enddo enddo endif end subroutine hqr2