@process directive('" (') subroutine hqr2(nm,n,low,igh,h,wr,wi,z,ierr) c integer i,j,k,l,m,n,en,na,nm, x igh,itn,its,low,enm2,ierr real h(nm,n),wr(n),wi(n),z(nm,n) real p,q,r,s,t,w,x,y,ra,sa,vi,vr,zz,norm,tst1,tst2,foo logical notlas c c this subroutine is a translation of the algol procedure hqr2, c num. math. 16, 181-204(1970) by peters and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). c c this subroutine finds the eigenvalues and eigenvectors c of a real upper hessenberg matrix by the qr method. the c eigenvectors of a real general matrix can also be found c if elmhes and eltran or orthes and ortran have c been used to reduce this general matrix to hessenberg form c and to accumulate the similarity transformations. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c low and igh are integers determined by the balancing c subroutine balanc. if balanc has not been used, c set low=1, igh=n. c c h contains the upper hessenberg matrix. c c z contains the transformation matrix produced by eltran c after the reduction by elmhes, or by ortran after the c reduction by orthes, if performed. if the eigenvectors c of the hessenberg matrix are desired, z must contain the c identity matrix. c c on output c c h has been destroyed. c c wr and wi contain the real and imaginary parts, c respectively, of the eigenvalues. the eigenvalues c are unordered except that complex conjugate pairs c of values appear consecutively with the eigenvalue c having the positive imaginary part first. if an c error exit is made, the eigenvalues should be correct c for indices ierr+1,...,n. c c z contains the real and imaginary parts of the eigenvectors. c if the i-th eigenvalue is real, the i-th column of z c contains its eigenvector. if the i-th eigenvalue is complex c with positive imaginary part, the i-th and (i+1)-th c columns of z contain the real and imaginary parts of its c eigenvector. the eigenvectors are unnormalized. if an c error exit is made, none of the eigenvectors has been found. c c ierr is set to c zero for normal return, c j if the limit of 30*n iterations is exhausted c while the j-th eigenvalue is being sought. c c calls cdiv for complex division. c c Questions and comments should be directed to Alan K. Cline, c Pleasant Valley Software, 8603 Altus Cove, Austin, TX 78759. c Electronic mail to cline@cs.utexas.edu. c c this version dated january 1989. (for the IBM 3090vf) c c ------------------------------------------------------------------ c call xuflow(0) ierr = 0 norm = 0.0e0 c .......... store roots isolated by balanc c and compute matrix norm .......... do 50 j = 1, n do 50 i = 1,min0(n,j+1) 50 norm = norm + abs(h(i,j)) c c" ( prefer vector do 52 i = 1, low-1 wr(i) = h(i,i) wi(i) = 0.0e0 52 continue c c" ( prefer vector do 54 i = igh+1,n wr(i) = h(i,i) wi(i) = 0.0e0 54 continue c en = igh t = 0.0e0 itn = 30*n c .......... search for next eigenvalues .......... 60 if (en .lt. low) go to 340 its = 0 na = en - 1 enm2 = na - 1 c .......... look for single small sub-diagonal element c for l=en step -1 until low do -- .......... 70 do 80 l = en, low+1, -1 s = abs(h(l-1,l-1)) + abs(h(l,l)) if (s .eq. 0.0e0) s = norm tst1 = s tst2 = tst1 + abs(h(l,l-1)) if (tst2 .eq. tst1) go to 100 80 continue c .......... form shift .......... 100 x = h(en,en) if (l .eq. en) go to 270 y = h(na,na) w = h(en,na) * h(na,en) if (l .eq. na) go to 280 if (itn .eq. 0) go to 1000 if (its .ne. 10 .and. its .ne. 20) go to 130 c .......... form exceptional shift .......... t = t + x c c" ( prefer vector do 120 i = low, en 120 h(i,i) = h(i,i) - x c s = abs(h(en,na)) + abs(h(na,enm2)) x = 0.75e0 * s y = x w = -0.4375e0 * s * s 130 its = its + 1 itn = itn - 1 c .......... look for two consecutive small c sub-diagonal elements. c for m=en-2 step -1 until l do -- .......... do 140 m = en-2, l, -1 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 .eq. l) go to 150 tst1 = abs(p)*(abs(h(m-1,m-1)) + abs(zz) + abs(h(m+1,m+1))) tst2 = tst1 + abs(h(m,m-1))*(abs(q) + abs(r)) if (tst2 .eq. tst1) go to 150 140 continue c 150 h(m+2,m) = 0.0e0 c" ( prefer vector do 160 i = m+3, en h(i,i-2) = 0.0e0 h(i,i-3) = 0.0e0 160 continue c .......... double qr step involving rows l to en and c columns m to en .......... do 260 k = m, na notlas = k .ne. na if (k .eq. m) go to 170 p = h(k,k-1) q = h(k+1,k-1) r = 0.0e0 if (notlas) r = h(k+2,k-1) x = abs(p) + abs(q) + abs(r) if (x .eq. 0.0e0) go to 260 p = p / x q = q / x r = r / x 170 s = sign(sqrt(p*p+q*q+r*r),p) if (k .eq. m) go to 180 h(k,k-1) = -s * x go to 190 180 if (l .ne. m) h(k,k-1) = -h(k,k-1) 190 p = p + s x = p / s y = q / s zz = r / s q = q / p r = r / p if (notlas) go to 225 c .......... row modification .......... c" ( prefer vector do 200 j = k, n foo = h(k,j) + q * h(k+1,j) h(k,j) = h(k,j) - foo * x h(k+1,j) = h(k+1,j) - foo * y 200 continue c j = min0(en,k+3) c .......... column modification .......... do 210 i = 1, j foo = x * h(i,k) + y * h(i,k+1) h(i,k) = h(i,k) - foo h(i,k+1) = h(i,k+1) - foo * q 210 continue c .......... accumulate transformations .......... do 220 i = low, igh foo = x * z(i,k) + y * z(i,k+1) z(i,k) = z(i,k) - foo z(i,k+1) = z(i,k+1) - foo * q 220 continue go to 255 225 continue c .......... row modification .......... c" ( prefer vector do 230 j = k, n foo = h(k,j) + q * h(k+1,j) + r * h(k+2,j) h(k,j) = h(k,j) - foo * x h(k+1,j) = h(k+1,j) - foo * y h(k+2,j) = h(k+2,j) - foo * zz 230 continue c j = min0(en,k+3) c .......... column modification .......... do 240 i = 1, j foo = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2) h(i,k) = h(i,k) - foo h(i,k+1) = h(i,k+1) - foo * q h(i,k+2) = h(i,k+2) - foo * r 240 continue c .......... accumulate transformations .......... do 250 i = low, igh foo = x * z(i,k) + y * z(i,k+1) + zz * z(i,k+2) z(i,k) = z(i,k) - foo z(i,k+1) = z(i,k+1) - foo * q z(i,k+2) = z(i,k+2) - foo * r 250 continue 255 continue c 260 continue c go to 70 c .......... one root found .......... 270 h(en,en) = x + t wr(en) = h(en,en) wi(en) = 0.0e0 en = na go to 60 c .......... two roots found .......... 280 p = (y - x) / 2.0e0 q = p * p + w zz = sqrt(abs(q)) h(en,en) = x + t x = h(en,en) h(na,na) = y + t if (q .lt. 0.0e0) go to 320 c .......... real pair .......... zz = p + sign(zz,p) wr(na) = x + zz wr(en) = wr(na) if (zz .ne. 0.0e0) wr(en) = x - w / zz wi(na) = 0.0e0 wi(en) = 0.0e0 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 c .......... row modification .......... c" ( prefer vector do 290 j = na, n foo = h(na,j) h(na,j) = q * foo + p * h(en,j) h(en,j) = q * h(en,j) - p * foo 290 continue c .......... column modification .......... do 300 i = 1, en foo = h(i,na) h(i,na) = q * foo + p * h(i,en) h(i,en) = q * h(i,en) - p * foo 300 continue c .......... accumulate transformations .......... do 310 i = low, igh foo = z(i,na) z(i,na) = q * foo + p * z(i,en) z(i,en) = q * z(i,en) - p * foo 310 continue c go to 330 c .......... complex pair .......... 320 wr(na) = x + p wr(en) = x + p wi(na) = zz wi(en) = -zz 330 en = enm2 go to 60 c .......... all roots found. backsubstitute to find c vectors of upper triangular form .......... 340 if (norm .eq. 0.0e0) go to 1001 c .......... for en=n step -1 until 1 do -- .......... do 800 en = n, 1, -1 p = wr(en) q = wi(en) na = en - 1 if (q) 710, 600, 800 c .......... real vector .......... 600 m = en h(en,en) = 1.0e0 if (na .eq. 0) go to 800 c .......... for i=en-1 step -1 until 1 do -- .......... do 700 i = en-1, 1, -1 w = h(i,i) - p r = 0.0e0 c do 610 j = m, en 610 r = r + h(i,j) * h(j,en) c if (wi(i) .ge. 0.0e0) go to 630 zz = w s = r go to 700 630 m = i if (wi(i) .ne. 0.0e0) go to 640 t = w if (t .ne. 0.0e0) go to 635 tst1 = norm t = tst1 632 t = 0.01e0 * t tst2 = norm + t if (tst2 .gt. tst1) go to 632 635 h(i,en) = -r / t go to 680 c .......... solve real equations .......... 640 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) .le. abs(zz)) go to 650 h(i+1,en) = (-r - w * t) / x go to 680 650 h(i+1,en) = (-s - y * t) / zz c c .......... overflow control .......... 680 t = abs(h(i,en)) if (t .eq. 0.0e0) go to 700 tst1 = t tst2 = tst1 + 1.0e0/tst1 if (tst2 .gt. tst1) go to 700 do 690 j = i, en h(j,en) = h(j,en)/t 690 continue c 700 continue c .......... end real vector .......... go to 800 c .......... complex vector .......... 710 m = na c .......... last vector component chosen imaginary so that c eigenvector matrix is triangular .......... if (abs(h(en,na)) .le. abs(h(na,en))) go to 720 h(na,na) = q / h(en,na) h(na,en) = -(h(en,en) - p) / h(en,na) go to 730 720 call cdiv(0.0e0,-h(na,en),h(na,na)-p,q,h(na,na),h(na,en)) 730 h(en,na) = 0.0e0 h(en,en) = 1.0e0 enm2 = na - 1 if (enm2 .eq. 0) go to 800 c .......... for i=en-2 step -1 until 1 do -- .......... do 795 i = en-2, 1, -1 w = h(i,i) - p ra = 0.0e0 sa = 0.0e0 c do 760 j = m, en ra = ra + h(i,j) * h(j,na) sa = sa + h(i,j) * h(j,en) 760 continue c if (wi(i) .ge. 0.0e0) go to 770 zz = w r = ra s = sa go to 795 770 m = i if (wi(i) .ne. 0.0e0) go to 780 call cdiv(-ra,-sa,w,q,h(i,na),h(i,en)) go to 790 c .......... solve complex equations .......... 780 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.0e0 * q if (vr .ne. 0.0e0 .or. vi .ne. 0.0e0) go to 784 tst1 = norm * (abs(w) + abs(q) + abs(x) x + abs(y) + abs(zz)) vr = tst1 783 vr = 0.01e0 * vr tst2 = tst1 + vr if (tst2 .gt. tst1) go to 783 784 call cdiv(x*r-zz*ra+q*sa,x*s-zz*sa-q*ra,vr,vi, x h(i,na),h(i,en)) if (abs(x) .le. abs(zz) + abs(q)) go to 785 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 go to 790 785 call cdiv(-r-y*h(i,na),-s-y*h(i,en),zz,q, x h(i+1,na),h(i+1,en)) c c .......... overflow control .......... 790 t = amax1(abs(h(i,na)), abs(h(i,en))) if (t .eq. 0.0e0) go to 795 tst1 = t tst2 = tst1 + 1.0e0/tst1 if (tst2 .gt. tst1) go to 795 do 792 j = i, en h(j,na) = h(j,na)/t h(j,en) = h(j,en)/t 792 continue c 795 continue c .......... end complex vector .......... 800 continue c .......... end back substitution. c vectors of isolated roots .......... do 820 i = 1, low-1 c" ( prefer vector do 820 j = i, n 820 z(i,j) = h(i,j) c do 840 i = igh+1, n c" ( prefer vector do 840 j = i, n 840 z(i,j) = h(i,j) c c .......... multiply by transformation matrix to give c vectors of original full matrix. c for j=n step -1 until low do -- .......... do 880 j = n, low, -1 do 880 i = low, igh zz = 0.0e0 c" ( prefer vector do 860 k = low, min0(j,igh) 860 zz = zz + z(i,k) * h(k,j) z(i,j) = zz 880 continue c go to 1001 c .......... set error -- all eigenvalues have not c converged after 30*n iterations .......... 1000 ierr = en 1001 return end