subroutine hqr2(nm,n,low,igh,h,wr,wi,z,ierr) c integer i,j,k,l,m,n,en,ii,jj,ll,mm,na,nm,nn, x igh,itn,its,low,mp2,enm2,ierr double precision h(nm,n),wr(n),wi(n),z(nm,n) double precision p,q,r,s,t,w,x,y,ra,sa,vi,vr,zz,norm,tst1,tst2 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 burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 norm = 0.0d0 k = 1 c .......... store roots isolated by balanc c and compute matrix norm .......... do 50 i = 1, n c do 40 j = k, n 40 norm = norm + dabs(h(i,j)) c k = i if (i .ge. low .and. i .le. igh) go to 50 wr(i) = h(i,i) wi(i) = 0.0d0 50 continue c en = igh t = 0.0d0 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 ll = low, en l = en + low - ll if (l .eq. low) go to 100 s = dabs(h(l-1,l-1)) + dabs(h(l,l)) if (s .eq. 0.0d0) s = norm tst1 = s tst2 = tst1 + dabs(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 do 120 i = low, en 120 h(i,i) = h(i,i) - x c s = dabs(h(en,na)) + dabs(h(na,enm2)) x = 0.75d0 * s y = x w = -0.4375d0 * 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 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 = dabs(p) + dabs(q) + dabs(r) p = p / s q = q / s r = r / s if (m .eq. l) go to 150 tst1 = dabs(p)*(dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1))) tst2 = tst1 + dabs(h(m,m-1))*(dabs(q) + dabs(r)) if (tst2 .eq. tst1) go to 150 140 continue c 150 mp2 = m + 2 c do 160 i = mp2, en h(i,i-2) = 0.0d0 if (i .eq. mp2) go to 160 h(i,i-3) = 0.0d0 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.0d0 if (notlas) r = h(k+2,k-1) x = dabs(p) + dabs(q) + dabs(r) if (x .eq. 0.0d0) go to 260 p = p / x q = q / x r = r / x 170 s = dsign(dsqrt(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 .......... do 200 j = k, n p = h(k,j) + q * h(k+1,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y 200 continue c j = min0(en,k+3) c .......... column modification .......... do 210 i = 1, j p = x * h(i,k) + y * h(i,k+1) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q 210 continue c .......... accumulate transformations .......... do 220 i = low, igh p = x * z(i,k) + y * z(i,k+1) z(i,k) = z(i,k) - p z(i,k+1) = z(i,k+1) - p * q 220 continue go to 255 225 continue c .......... row modification .......... do 230 j = k, n p = h(k,j) + q * h(k+1,j) + r * h(k+2,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y h(k+2,j) = h(k+2,j) - p * zz 230 continue c j = min0(en,k+3) c .......... column modification .......... do 240 i = 1, j p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q h(i,k+2) = h(i,k+2) - p * r 240 continue c .......... accumulate transformations .......... do 250 i = low, igh p = x * z(i,k) + y * z(i,k+1) + zz * z(i,k+2) z(i,k) = z(i,k) - p z(i,k+1) = z(i,k+1) - p * q z(i,k+2) = z(i,k+2) - p * 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.0d0 en = na go to 60 c .......... two roots found .......... 280 p = (y - x) / 2.0d0 q = p * p + w zz = dsqrt(dabs(q)) h(en,en) = x + t x = h(en,en) h(na,na) = y + t if (q .lt. 0.0d0) go to 320 c .......... real pair .......... zz = p + dsign(zz,p) wr(na) = x + zz wr(en) = wr(na) if (zz .ne. 0.0d0) wr(en) = x - w / zz wi(na) = 0.0d0 wi(en) = 0.0d0 x = h(en,na) s = dabs(x) + dabs(zz) p = x / s q = zz / s r = dsqrt(p*p+q*q) p = p / r q = q / r c .......... row modification .......... do 290 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 290 continue c .......... column modification .......... do 300 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 300 continue c .......... accumulate transformations .......... do 310 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 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.0d0) go to 1001 c .......... for en=n step -1 until 1 do -- .......... do 800 nn = 1, n en = n + 1 - nn p = wr(en) q = wi(en) na = en - 1 if (q) 710, 600, 800 c .......... real vector .......... 600 m = en h(en,en) = 1.0d0 if (na .eq. 0) go to 800 c .......... for i=en-1 step -1 until 1 do -- .......... do 700 ii = 1, na i = en - ii w = h(i,i) - p r = 0.0d0 c do 610 j = m, en 610 r = r + h(i,j) * h(j,en) c if (wi(i) .ge. 0.0d0) go to 630 zz = w s = r go to 700 630 m = i if (wi(i) .ne. 0.0d0) go to 640 t = w if (t .ne. 0.0d0) go to 635 tst1 = norm t = tst1 632 t = 0.01d0 * 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 (dabs(x) .le. dabs(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 = dabs(h(i,en)) if (t .eq. 0.0d0) go to 700 tst1 = t tst2 = tst1 + 1.0d0/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 (dabs(h(en,na)) .le. dabs(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.0d0,-h(na,en),h(na,na)-p,q,h(na,na),h(na,en)) 730 h(en,na) = 0.0d0 h(en,en) = 1.0d0 enm2 = na - 1 if (enm2 .eq. 0) go to 800 c .......... for i=en-2 step -1 until 1 do -- .......... do 795 ii = 1, enm2 i = na - ii w = h(i,i) - p ra = 0.0d0 sa = 0.0d0 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.0d0) go to 770 zz = w r = ra s = sa go to 795 770 m = i if (wi(i) .ne. 0.0d0) 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.0d0 * q if (vr .ne. 0.0d0 .or. vi .ne. 0.0d0) go to 784 tst1 = norm * (dabs(w) + dabs(q) + dabs(x) x + dabs(y) + dabs(zz)) vr = tst1 783 vr = 0.01d0 * 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 (dabs(x) .le. dabs(zz) + dabs(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 = dmax1(dabs(h(i,na)), dabs(h(i,en))) if (t .eq. 0.0d0) go to 795 tst1 = t tst2 = tst1 + 1.0d0/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 840 i = 1, n if (i .ge. low .and. i .le. igh) go to 840 c do 820 j = i, n 820 z(i,j) = h(i,j) c 840 continue c .......... multiply by transformation matrix to give c vectors of original full matrix. c for j=n step -1 until low do -- .......... do 880 jj = low, n j = n + low - jj m = min0(j,igh) c do 880 i = low, igh zz = 0.0d0 c do 860 k = low, m 860 zz = zz + z(i,k) * h(k,j) c 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