C/ FIGURE 3.6.1 SUBROUTINE DEGENP(A,B,N,LAM0,POLY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DECLARATIONS FOR ARGUMENTS DOUBLE PRECISION A(N,N),B(N,N),LAM0,POLY(N+1) INTEGER N C DECLARATIONS FOR LOCAL VARIABLES DOUBLE PRECISION DET(N+1,N+1) C C SUBROUTINE DEGENP CALLS HESSQZ TO REDUCE THE GENERALIZED EIGENVALUE C PROBLEM C C A*X = LAMBDA*B*X C C TO A SIMILAR PROBLEM WITH THE SAME EIGENVALUES, WHERE A IS UPPER C HESSENBERG AND B IS UPPER TRIANGULAR, THEN CALCULATES THE COEFFICIENTS C OF THE POLYNOMIAL DET(A-LAMBDA*B). THE ROOTS OF THIS POLYNOMIAL WILL C BE THE EIGENVALUES OF THE GENERALIZED PROBLEM. C C ARGUMENTS C C ON INPUT ON OUTPUT C -------- --------- C C A - THE N BY N A MATRIX. A IS UPPER HESSENBERG. C C B - THE N BY N B MATRIX. B IS UPPER TRIANGULAR. C C N - THE SIZE OF MATRICES A AND B. C C LAM0 - A SCALAR, SEE POLY. (SET C LAM0=0 TO GET THE USUAL C POLYNOMIAL COEFFICIENTS) C C POLY - VECTOR OF LENGTH N+1, C CONTAINING POLYNOMIAL C COEFFICIENTS. C DET(A-LAMBDA*B)= C SUM FROM I=0 TO N OF C POLY(I+1)* C (LAMBDA-LAM0)**I C C----------------------------------------------------------------------- C C CALL HESSQZ TO REDUCE A TO UPPER HESSENBERG AND C B TO UPPER TRIANGULAR FORM CALL HESSQZ(A,B,N) C DEGENP USES A RECURRENCE RELATION TO CALCULATE C DET(A-LAMBDA*B) AND ALL N DERIVATIVES AT C LAMBDA = LAM0, FROM WHICH THE POLYNOMIAL C COEFFICIENTS CAN BE FOUND. DET(1,N+1) = 1.0 DO 5 I=2,N+1 DET(I,N+1) = 0.0 5 CONTINUE DET(1,N) = A(N,N)-LAM0*B(N,N) DET(2,N) = -B(N,N) DO 10 I=3,N+1 DET(I,N) = 0.0 10 CONTINUE DO 30 K=N-1,1,-1 DET(1,K) = (A(K,K)-LAM0*B(K,K))*DET(1,K+1) DO 15 I=1,N DET(I+1,K) = (A(K,K)-LAM0*B(K,K))*DET(I+1,K+1) & - B(K,K)*DET(I,K+1) 15 CONTINUE FACT = 1.0 DO 25 J=K+1,N FACT = -FACT*A(J,J-1) IF (A(K,J).EQ.0.0 .AND. B(K,J).EQ.0.0) GO TO 25 DET(1,K) = DET(1,K) + FACT*(A(K,J)-LAM0*B(K,J))*DET(1,J+1) DO 20 I=1,N DET(I+1,K) = DET(I+1,K) + FACT* & ((A(K,J)-LAM0*B(K,J))*DET(I+1,J+1) - B(K,J)*DET(I,J+1)) 20 CONTINUE 25 CONTINUE 30 CONTINUE DO 35 I=1,N+1 POLY(I) = DET(I,1) 35 CONTINUE RETURN END SUBROUTINE HESSQZ(A,B,N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DECLARATIONS FOR ARGUMENTS DOUBLE PRECISION A(N,N),B(N,N) INTEGER N C C SUBROUTINE HESSQZ REDUCES THE GENERALIZED EIGENVALUE PROBLEM C C A*X = LAMBDA*B*X C C TO A SIMILAR PROBLEM WITH THE SAME EIGENVALUES, WHERE A IS C UPPER HESSENBERG AND B IS UPPER TRIANGULAR. C C ARGUMENTS C C ON INPUT ON OUTPUT C -------- --------- C C A - THE N BY N A MATRIX. A IS UPPER HESSENBERG. C C B - THE N BY N B MATRIX. B IS UPPER TRIANGULAR. C C N - THE SIZE OF MATRICES A AND B. C C----------------------------------------------------------------------- C PREMULTIPLY A AND B BY ORTHOGONAL MATRIX C (PRODUCT OF GIVENS MATRICES) Q, SUCH C THAT QB IS UPPER TRIANGULAR. DO 20 I=1,N-1 DO 15 J=I+1,N IF (B(J,I).EQ.0.0) GO TO 15 DEN = SQRT(B(I,I)**2+B(J,I)**2) S = -B(J,I)/DEN C = B(I,I)/DEN DO 5 K=I,N BIK = C*B(I,K)-S*B(J,K) BJK = S*B(I,K)+C*B(J,K) B(I,K) = BIK B(J,K) = BJK 5 CONTINUE DO 10 K=1,N AIK = C*A(I,K)-S*A(J,K) AJK = S*A(I,K)+C*A(J,K) A(I,K) = AIK A(J,K) = AJK 10 CONTINUE 15 CONTINUE 20 CONTINUE C PREMULTIPLY A AND B BY ORTHOGONAL MATRIX C Q, AND POSTMULTIPLY BY ORTHOGONAL MATRIX C Z, SUCH THAT QAZ IS UPPER HESSENBERG AND C QBZ IS STILL UPPER TRIANGULAR DO 50 I=1,N-2 DO 45 J=N,I+2,-1 IF (A(J,I).EQ.0.0) GO TO 45 C PREMULTIPLY A TO ZERO A(J,I) DEN = SQRT(A(J-1,I)**2+A(J,I)**2) S = -A(J,I)/DEN C = A(J-1,I)/DEN DO 25 K=I,N A1K = C*A(J-1,K) - S*A(J,K) A2K = S*A(J-1,K) + C*A(J,K) A(J-1,K) = A1K A(J,K) = A2K 25 CONTINUE C PREMULTIPLY B BY SAME MATRIX, CREATING C NEW NONZERO B(J,J-1) DO 30 K=J-1,N B1K = C*B(J-1,K) - S*B(J,K) B2K = S*B(J-1,K) + C*B(J,K) B(J-1,K) = B1K B(J,K) = B2K 30 CONTINUE IF (B(J,J-1).EQ.0.0) GO TO 45 C POSTMULTIPLY B TO ZERO B(J,J-1) DEN = SQRT(B(J,J-1)**2+B(J,J)**2) S = -B(J,J-1)/DEN C = B(J,J)/DEN DO 35 K=1,J BK1 = C*B(K,J-1) + S*B(K,J) BK2 = -S*B(K,J-1) + C*B(K,J) B(K,J-1) = BK1 B(K,J) = BK2 35 CONTINUE C POSTMULTIPLY A BY SAME MATRIX DO 40 K=1,N AK1 = C*A(K,J-1) + S*A(K,J) AK2 = -S*A(K,J-1) + C*A(K,J) A(K,J-1) = AK1 A(K,J) = AK2 40 CONTINUE 45 CONTINUE 50 CONTINUE RETURN END SUBROUTINE DEGEN(A,B,N,EIG,NEIG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DECLARATIONS FOR ARGUMENTS COMPLEX*16 EIG(N) DOUBLE PRECISION A(N,N),B(N,N) INTEGER N,NEIG C DECLARATIONS FOR LOCAL VARIABLES DOUBLE PRECISION ALFR(N),ALFI(N),BETA(N) C C SUBROUTINE DEGEN FINDS ALL EIGENVALUES OF THE GENERALIZED C EIGENVALUE PROBLEM C C A*X = LAMBDA*B*X C C ARGUMENTS C C ON INPUT ON OUTPUT C -------- --------- C C A - THE N BY N A MATRIX. DESTROYED. C C B - THE N BY N B MATRIX. DESTROYED. C C N - THE SIZE OF MATRICES A AND B. C C EIG - A COMPLEX N-VECTOR C CONTAINING THE C EIGENVALUES IN ITS C FIRST NEIG ELEMENTS C C NEIG - NUMBER OF EIGENVALUES C (NEIG .LE. N) C C----------------------------------------------------------------------- C EPS = MACHINE FLOATING POINT RELATIVE C PRECISION C ***************************** DATA EPS/2.D-16/ C ***************************** EPS1 = 100*EPS C REDUCE A TO UPPER HESSENBERG, B TO C UPPER TRIANGULAR FORM CALL HESSQZ(A,B,N) C REDUCE A TO QUASITRIANGULAR FORM CALL QZIT(N,A,B,EPS1,IERR) IF (IERR.NE.0) GO TO 20 C EXTRACT EIGENVALUES CALL QZVAL(N,A,B,ALFR,ALFI,BETA) NEIG = 0 EPSB = B(N,1) DO 10 I=1,N IF (BETA(I).GT.EPSB) THEN NEIG = NEIG+1 EIG(NEIG) = (ALFR(I) + ALFI(I)*CMPLX(0.0,1.0))/BETA(I) ENDIF 10 CONTINUE RETURN C QZ ITERATION DOES NOT CONVERGE 20 PRINT 30 30 FORMAT (' ***** QZ ITERATION DOES NOT CONVERGE *****') RETURN END subroutine qzit(n,a,b,eps1,ierr) c integer i,j,k,l,n,en,k1,k2,ld,ll,l1,na,ish,itn,its,km1,lm1, x enm2,ierr,lor1,enorn double precision a(n,n),b(n,n) double precision r,s,t,a1,a2,a3,ep,sh,u1,u2,u3,v1,v2,v3,ani,a11, x a12,a21,a22,a33,a34,a43,a44,bni,b11,b12,b22,b33,b34, x b44,epsa,epsb,eps1,anorm,bnorm logical notlas c c this subroutine is the second step of the qz algorithm c for solving generalized matrix eigenvalue problems, c siam j. numer. anal. 10, 241-256(1973) by moler and stewart, c as modified in technical note nasa tn d-7305(1973) by ward. c c this subroutine accepts a pair of real matrices, one of them c in upper hessenberg form and the other in upper triangular form. c it reduces the hessenberg matrix to quasi-triangular form using c orthogonal transformations while maintaining the triangular form c of the other matrix. c c on input c c n is the order of the matrices. c c a contains a real upper hessenberg matrix. c c b contains a real upper triangular matrix. c c eps1 is a tolerance used to determine negligible elements. c c on output c c a has been reduced to quasi-triangular form. the elements c below the first subdiagonal are still zero and no two c consecutive subdiagonal elements are nonzero. c c b is still in upper triangular form, although its elements c have been altered. the location b(n,1) is used to store c eps1 times the norm of b for later use by qzval and qzvec. 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 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 ------------------------------------------------------------------ ierr = 0 c .......... compute epsa,epsb .......... anorm = 0.0d0 bnorm = 0.0d0 c do 30 i = 1, n ani = 0.0d0 if (i .ne. 1) ani = abs(a(i,i-1)) bni = 0.0d0 c do 20 j = i, n ani = ani + abs(a(i,j)) bni = bni + abs(b(i,j)) 20 continue c if (ani .gt. anorm) anorm = ani if (bni .gt. bnorm) bnorm = bni 30 continue c if (anorm .eq. 0.0d0) anorm = 1.0d0 if (bnorm .eq. 0.0d0) bnorm = 1.0d0 ep = eps1 50 epsa = ep * anorm epsb = ep * bnorm c .......... reduce a to quasi-triangular form, while c keeping b triangular .......... lor1 = 1 enorn = n en = n itn = 30*n c .......... begin qz step .......... 60 if (en .le. 2) go to 1001 enorn = en its = 0 na = en - 1 enm2 = na - 1 70 ish = 2 c .......... check for convergence or reducibility. c for l=en step -1 until 1 do -- .......... do 80 ll = 1, en lm1 = en - ll l = lm1 + 1 if (l .eq. 1) go to 95 if (abs(a(l,lm1)) .le. epsa) go to 90 80 continue c 90 a(l,lm1) = 0.0d0 if (l .lt. na) go to 95 c .......... 1-by-1 or 2-by-2 block isolated .......... en = lm1 go to 60 c .......... check for small top of b .......... 95 ld = l 100 l1 = l + 1 b11 = b(l,l) if (abs(b11) .gt. epsb) go to 120 b(l,l) = 0.0d0 s = abs(a(l,l)) + abs(a(l1,l)) u1 = a(l,l) / s u2 = a(l1,l) / s r = sign(sqrt(u1*u1+u2*u2),u1) v1 = -(u1 + r) / r v2 = -u2 / r u2 = v2 / v1 c do 110 j = l, enorn t = a(l,j) + u2 * a(l1,j) a(l,j) = a(l,j) + t * v1 a(l1,j) = a(l1,j) + t * v2 t = b(l,j) + u2 * b(l1,j) b(l,j) = b(l,j) + t * v1 b(l1,j) = b(l1,j) + t * v2 110 continue c if (l .ne. 1) a(l,lm1) = -a(l,lm1) lm1 = l l = l1 go to 90 120 a11 = a(l,l) / b11 a21 = a(l1,l) / b11 if (ish .eq. 1) go to 140 c .......... iteration strategy .......... if (itn .eq. 0) go to 1000 if (its .eq. 10) go to 155 c .......... determine type of shift .......... b22 = b(l1,l1) if (abs(b22) .lt. epsb) b22 = epsb b33 = b(na,na) if (abs(b33) .lt. epsb) b33 = epsb b44 = b(en,en) if (abs(b44) .lt. epsb) b44 = epsb a33 = a(na,na) / b33 a34 = a(na,en) / b44 a43 = a(en,na) / b33 a44 = a(en,en) / b44 b34 = b(na,en) / b44 t = 0.5d0 * (a43 * b34 - a33 - a44) r = t * t + a34 * a43 - a33 * a44 if (r .lt. 0.0d0) go to 150 c .......... determine single shift zeroth column of a .......... ish = 1 r = sqrt(r) sh = -t + r s = -t - r if (abs(s-a44) .lt. abs(sh-a44)) sh = s c .......... look for two consecutive small c sub-diagonal elements of a. c for l=en-2 step -1 until ld do -- .......... do 130 ll = ld, enm2 l = enm2 + ld - ll if (l .eq. ld) go to 140 lm1 = l - 1 l1 = l + 1 t = a(l,l) if (abs(b(l,l)) .gt. epsb) t = t - sh * b(l,l) if (abs(a(l,lm1)) .le. abs(t/a(l1,l)) * epsa) go to 100 130 continue c 140 a1 = a11 - sh a2 = a21 if (l .ne. ld) a(l,lm1) = -a(l,lm1) go to 160 c .......... determine double shift zeroth column of a .......... 150 a12 = a(l,l1) / b22 a22 = a(l1,l1) / b22 b12 = b(l,l1) / b22 a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11) x / a21 + a12 - a11 * b12 a2 = (a22 - a11) - a21 * b12 - (a33 - a11) - (a44 - a11) x + a43 * b34 a3 = a(l1+1,l1) / b22 go to 160 c .......... ad hoc shift .......... 155 a1 = 0.0d0 a2 = 1.0d0 a3 = 1.1605d0 160 its = its + 1 itn = itn - 1 lor1 = ld c .......... main loop .......... do 260 k = l, na notlas = k .ne. na .and. ish .eq. 2 k1 = k + 1 k2 = k + 2 km1 = max(k-1,l) ll = min(en,k1+ish) if (notlas) go to 190 c .......... zero a(k+1,k-1) .......... if (k .eq. l) go to 170 a1 = a(k,km1) a2 = a(k1,km1) 170 s = abs(a1) + abs(a2) if (s .eq. 0.0d0) go to 70 u1 = a1 / s u2 = a2 / s r = sign(sqrt(u1*u1+u2*u2),u1) v1 = -(u1 + r) / r v2 = -u2 / r u2 = v2 / v1 c do 180 j = km1, enorn t = a(k,j) + u2 * a(k1,j) a(k,j) = a(k,j) + t * v1 a(k1,j) = a(k1,j) + t * v2 t = b(k,j) + u2 * b(k1,j) b(k,j) = b(k,j) + t * v1 b(k1,j) = b(k1,j) + t * v2 180 continue c if (k .ne. l) a(k1,km1) = 0.0d0 go to 240 c .......... zero a(k+1,k-1) and a(k+2,k-1) .......... 190 if (k .eq. l) go to 200 a1 = a(k,km1) a2 = a(k1,km1) a3 = a(k2,km1) 200 s = abs(a1) + abs(a2) + abs(a3) if (s .eq. 0.0d0) go to 260 u1 = a1 / s u2 = a2 / s u3 = a3 / s r = sign(sqrt(u1*u1+u2*u2+u3*u3),u1) v1 = -(u1 + r) / r v2 = -u2 / r v3 = -u3 / r u2 = v2 / v1 u3 = v3 / v1 c do 210 j = km1, enorn t = a(k,j) + u2 * a(k1,j) + u3 * a(k2,j) a(k,j) = a(k,j) + t * v1 a(k1,j) = a(k1,j) + t * v2 a(k2,j) = a(k2,j) + t * v3 t = b(k,j) + u2 * b(k1,j) + u3 * b(k2,j) b(k,j) = b(k,j) + t * v1 b(k1,j) = b(k1,j) + t * v2 b(k2,j) = b(k2,j) + t * v3 210 continue c if (k .eq. l) go to 220 a(k1,km1) = 0.0d0 a(k2,km1) = 0.0d0 c .......... zero b(k+2,k+1) and b(k+2,k) .......... 220 s = abs(b(k2,k2)) + abs(b(k2,k1)) + abs(b(k2,k)) if (s .eq. 0.0d0) go to 240 u1 = b(k2,k2) / s u2 = b(k2,k1) / s u3 = b(k2,k) / s r = sign(sqrt(u1*u1+u2*u2+u3*u3),u1) v1 = -(u1 + r) / r v2 = -u2 / r v3 = -u3 / r u2 = v2 / v1 u3 = v3 / v1 c do 230 i = lor1, ll t = a(i,k2) + u2 * a(i,k1) + u3 * a(i,k) a(i,k2) = a(i,k2) + t * v1 a(i,k1) = a(i,k1) + t * v2 a(i,k) = a(i,k) + t * v3 t = b(i,k2) + u2 * b(i,k1) + u3 * b(i,k) b(i,k2) = b(i,k2) + t * v1 b(i,k1) = b(i,k1) + t * v2 b(i,k) = b(i,k) + t * v3 230 continue c b(k2,k) = 0.0d0 b(k2,k1) = 0.0d0 c .......... zero b(k+1,k) .......... 240 s = abs(b(k1,k1)) + abs(b(k1,k)) if (s .eq. 0.0d0) go to 260 u1 = b(k1,k1) / s u2 = b(k1,k) / s r = sign(sqrt(u1*u1+u2*u2),u1) v1 = -(u1 + r) / r v2 = -u2 / r u2 = v2 / v1 c do 250 i = lor1, ll t = a(i,k1) + u2 * a(i,k) a(i,k1) = a(i,k1) + t * v1 a(i,k) = a(i,k) + t * v2 t = b(i,k1) + u2 * b(i,k) b(i,k1) = b(i,k1) + t * v1 b(i,k) = b(i,k) + t * v2 250 continue c b(k1,k) = 0.0d0 260 continue c .......... end qz step .......... go to 70 c .......... set error -- all eigenvalues have not c converged after 30*n iterations .......... 1000 ierr = en c .......... save epsb for use by qzval and qzvec .......... 1001 if (n .gt. 1) b(n,1) = epsb return end subroutine qzval(n,a,b,alfr,alfi,beta) c integer i,j,n,en,na,nn,isw double precision a(n,n),b(n,n),alfr(n),alfi(n),beta(n) double precision c,d,e,r,s,t,an,a1,a2,bn,cq,cz,di,dr,ei,ti,tr,u1, x u2,v1,v2,a1i,a11,a12,a2i,a21,a22,b11,b12,b22,sqi,sqr, x ssi,ssr,szi,szr,a11i,a11r,a12i,a12r,a22i,a22r,epsb c c this subroutine is the third step of the qz algorithm c for solving generalized matrix eigenvalue problems, c siam j. numer. anal. 10, 241-256(1973) by moler and stewart. c c this subroutine accepts a pair of real matrices, one of them c in quasi-triangular form and the other in upper triangular form. c it reduces the quasi-triangular matrix further, so that any c remaining 2-by-2 blocks correspond to pairs of complex c eigenvalues, and returns quantities whose ratios give the c generalized eigenvalues. c c on input c c n is the order of the matrices. c c a contains a real upper quasi-triangular matrix. c c b contains a real upper triangular matrix. in addition, c location b(n,1) contains the tolerance quantity (epsb) c computed and saved in qzit. c c on output c c a has been reduced further to a quasi-triangular matrix c in which all nonzero subdiagonal elements correspond to c pairs of complex eigenvalues. c c b is still in upper triangular form, although its elements c have been altered. b(n,1) is unaltered. c c alfr and alfi contain the real and imaginary parts of the c diagonal elements of the triangular matrix that would be c obtained if a were reduced completely to triangular form c by unitary transformations. non-zero values of alfi occur c in pairs, the first member positive and the second negative. c c beta contains the diagonal elements of the corresponding b, c normalized to be real and non-negative. the generalized c eigenvalues are then the ratios ((alfr+i*alfi)/beta). 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 epsb = b(n,1) isw = 1 c .......... find eigenvalues of quasi-triangular matrices. c for en=n step -1 until 1 do -- .......... do 510 nn = 1, n en = n + 1 - nn na = en - 1 if (isw .eq. 2) go to 505 if (en .eq. 1) go to 410 if (a(en,na) .ne. 0.0d0) go to 420 c .......... 1-by-1 block, one real root .......... 410 alfr(en) = a(en,en) if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en) beta(en) = abs(b(en,en)) alfi(en) = 0.0d0 go to 510 c .......... 2-by-2 block .......... 420 if (abs(b(na,na)) .le. epsb) go to 455 if (abs(b(en,en)) .gt. epsb) go to 430 a1 = a(en,en) a2 = a(en,na) bn = 0.0d0 go to 435 430 an = abs(a(na,na)) + abs(a(na,en)) + abs(a(en,na)) x + abs(a(en,en)) bn = abs(b(na,na)) + abs(b(na,en)) + abs(b(en,en)) a11 = a(na,na) / an a12 = a(na,en) / an a21 = a(en,na) / an a22 = a(en,en) / an b11 = b(na,na) / bn b12 = b(na,en) / bn b22 = b(en,en) / bn e = a11 / b11 ei = a22 / b22 s = a21 / (b11 * b22) t = (a22 - e * b22) / b22 if (abs(e) .le. abs(ei)) go to 431 e = ei t = (a11 - e * b11) / b11 431 c = 0.5d0 * (t - s * b12) d = c * c + s * (a12 - e * b12) if (d .lt. 0.0d0) go to 480 c .......... two real roots. c zero both a(en,na) and b(en,na) .......... e = e + (c + sign(sqrt(d),c)) a11 = a11 - e * b11 a12 = a12 - e * b12 a22 = a22 - e * b22 if (abs(a11) + abs(a12) .lt. x abs(a21) + abs(a22)) go to 432 a1 = a12 a2 = a11 go to 435 432 a1 = a22 a2 = a21 c .......... choose and apply real z .......... 435 s = abs(a1) + abs(a2) u1 = a1 / s u2 = a2 / s r = sign(sqrt(u1*u1+u2*u2),u1) v1 = -(u1 + r) / r v2 = -u2 / r u2 = v2 / v1 c do 440 i = 1, en t = a(i,en) + u2 * a(i,na) a(i,en) = a(i,en) + t * v1 a(i,na) = a(i,na) + t * v2 t = b(i,en) + u2 * b(i,na) b(i,en) = b(i,en) + t * v1 b(i,na) = b(i,na) + t * v2 440 continue c if (bn .eq. 0.0d0) go to 475 if (an .lt. abs(e) * bn) go to 455 a1 = b(na,na) a2 = b(en,na) go to 460 455 a1 = a(na,na) a2 = a(en,na) c .......... choose and apply real q .......... 460 s = abs(a1) + abs(a2) if (s .eq. 0.0d0) go to 475 u1 = a1 / s u2 = a2 / s r = sign(sqrt(u1*u1+u2*u2),u1) v1 = -(u1 + r) / r v2 = -u2 / r u2 = v2 / v1 c do 470 j = na, n t = a(na,j) + u2 * a(en,j) a(na,j) = a(na,j) + t * v1 a(en,j) = a(en,j) + t * v2 t = b(na,j) + u2 * b(en,j) b(na,j) = b(na,j) + t * v1 b(en,j) = b(en,j) + t * v2 470 continue c 475 a(en,na) = 0.0d0 b(en,na) = 0.0d0 alfr(na) = a(na,na) alfr(en) = a(en,en) if (b(na,na) .lt. 0.0d0) alfr(na) = -alfr(na) if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en) beta(na) = abs(b(na,na)) beta(en) = abs(b(en,en)) alfi(en) = 0.0d0 alfi(na) = 0.0d0 go to 505 c .......... two complex roots .......... 480 e = e + c ei = sqrt(-d) a11r = a11 - e * b11 a11i = ei * b11 a12r = a12 - e * b12 a12i = ei * b12 a22r = a22 - e * b22 a22i = ei * b22 if (abs(a11r) + abs(a11i) + abs(a12r) + abs(a12i) .lt. x abs(a21) + abs(a22r) + abs(a22i)) go to 482 a1 = a12r a1i = a12i a2 = -a11r a2i = -a11i go to 485 482 a1 = a22r a1i = a22i a2 = -a21 a2i = 0.0d0 c .......... choose complex z .......... 485 cz = sqrt(a1*a1+a1i*a1i) if (cz .eq. 0.0d0) go to 487 szr = (a1 * a2 + a1i * a2i) / cz szi = (a1 * a2i - a1i * a2) / cz r = sqrt(cz*cz+szr*szr+szi*szi) cz = cz / r szr = szr / r szi = szi / r go to 490 487 szr = 1.0d0 szi = 0.0d0 490 if (an .lt. (abs(e) + ei) * bn) go to 492 a1 = cz * b11 + szr * b12 a1i = szi * b12 a2 = szr * b22 a2i = szi * b22 go to 495 492 a1 = cz * a11 + szr * a12 a1i = szi * a12 a2 = cz * a21 + szr * a22 a2i = szi * a22 c .......... choose complex q .......... 495 cq = sqrt(a1*a1+a1i*a1i) if (cq .eq. 0.0d0) go to 497 sqr = (a1 * a2 + a1i * a2i) / cq sqi = (a1 * a2i - a1i * a2) / cq r = sqrt(cq*cq+sqr*sqr+sqi*sqi) cq = cq / r sqr = sqr / r sqi = sqi / r go to 500 497 sqr = 1.0d0 sqi = 0.0d0 c .......... compute diagonal elements that would result c if transformations were applied .......... 500 ssr = sqr * szr + sqi * szi ssi = sqr * szi - sqi * szr i = 1 tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21 x + ssr * a22 ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22 dr = cq * cz * b11 + cq * szr * b12 + ssr * b22 di = cq * szi * b12 + ssi * b22 go to 503 502 i = 2 tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21 x + cq * cz * a22 ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21 dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22 di = -ssi * b11 - sqi * cz * b12 503 t = ti * dr - tr * di j = na if (t .lt. 0.0d0) j = en r = sqrt(dr*dr+di*di) beta(j) = bn * r alfr(j) = an * (tr * dr + ti * di) / r alfi(j) = an * t / r if (i .eq. 1) go to 502 505 isw = 3 - isw 510 continue b(n,1) = epsb c return end