C/ FIGURE 3.5.2 SUBROUTINE DPOWER(A,N,EIG,V,IUPDAT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DECLARATIONS FOR ARGUMENTS COMPLEX*16 EIG,V(N) DOUBLE PRECISION A(N,N) INTEGER N,IUPDAT C DECLARATIONS FOR LOCAL VARIABLES COMPLEX*16 VN(N),VNP1(N),B(N,N),PN,R,RNUM,RDEN INTEGER IPERM(N) C C SUBROUTINE DPOWER FINDS ONE EIGENVALUE OF A, AND A CORRESPONDING C EIGENVECTOR, USING THE SHIFTED INVERSE POWER METHOD. C C ARGUMENTS C C ON INPUT ON OUTPUT C -------- --------- C C A - THE N BY N MATRIX. C C N - THE SIZE OF MATRIX A. C C EIG - A (COMPLEX) INITIAL GUESS AT AN EIGENVALUE OF A, C AN EIGENVALUE. NORMALLY THE ONE CLOSEST C TO THE INITIAL GUESS. C C V - A (COMPLEX) STARTING VECTOR AN EIGENVECTOR OF A, C FOR THE SHIFTED INVERSE POWER CORRESPONDING TO THE C METHOD. IF ALL COMPONENTS OF COMPUTED EIGENVALUE. C V ARE ZERO ON INPUT, A RANDOM C STARTING VECTOR WILL BE USED. C C IUPDAT - THE NUMBER OF SHIFTED INVERSE C POWER ITERATIONS TO BE DONE C BETWEEN UPDATES OF P. IF C IUPDAT=1, P WILL BE UPDATED EVERY C ITERATION. IF IUPDAT > 1000, C P WILL NEVER BE UPDATED. C C----------------------------------------------------------------------- C EPS = MACHINE FLOATING POINT RELATIVE C PRECISION C ***************************** DATA EPS/2.D-16/ C ***************************** DO 5 I=1,N IF (V(I).NE.0.0) GO TO 15 5 CONTINUE C IF V = 0, GENERATE A RANDOM STARTING C VECTOR SEED = N+10000 DEN = 2.0**31-1.0 DO 10 I=1,N SEED = MOD(7**5*SEED,DEN) V(I) = SEED/(DEN+1.0) 10 CONTINUE 15 CONTINUE C NORMALIZE V, AND SET VN=V VNORM = 0.0 DO 20 I=1,N VNORM = VNORM + ABS(V(I))**2 20 CONTINUE VNORM = SQRT(VNORM) DO 25 I=1,N V(I) = V(I)/VNORM VN(I) = V(I) 25 CONTINUE C BEGIN SHIFTED INVERSE POWER ITERATION NITER = 1000 DO 60 ITER=0,NITER IF (MOD(ITER,IUPDAT).EQ.0) THEN C EVERY IUPDAT ITERATIONS, UPDATE PN C AND SOLVE (A-PN*I)*VNP1 = VN PN = EIG DO 30 I=1,N DO 30 J=1,N IF (I.EQ.J) THEN B(I,J) = A(I,J) - PN ELSE B(I,J) = A(I,J) ENDIF 30 CONTINUE CALL CLINEQ(B,N,VNP1,V,IPERM) ELSE C BETWEEN UPDATES, WE CAN USE THE LU C DECOMPOSITION OF B=A-PN*I CALCULATED C EARLIER, TO SOLVE B*VNP1=VN FASTER CALL CRESLV(B,N,VNP1,V,IPERM) ENDIF C CALCULATE NEW EIGENVALUE ESTIMATE, C PN + (VN*VN)/(VN*VNP1) RNUM = 0.0 RDEN = 0.0 DO 35 I=1,N RNUM = RNUM + VN(I)*VN(I) RDEN = RDEN + VN(I)*VNP1(I) 35 CONTINUE R = RNUM/RDEN EIG = PN + R C SET V = NORMALIZED VNP1 VNORM = 0.0 DO 40 I=1,N VNORM = VNORM + ABS(VNP1(I))**2 40 CONTINUE VNORM = SQRT(VNORM) DO 45 I=1,N V(I) = VNP1(I)/VNORM 45 CONTINUE C IF R*VNP1 = VN (R = (VN*VN)/(VN*VNP1) ), C ITERATION HAS CONVERGED. ERRMAX = 0.0 DO 50 I=1,N ERRMAX = MAX(ERRMAX,ABS(R*VNP1(I)-VN(I))) 50 CONTINUE IF (ERRMAX.LE.SQRT(EPS)) RETURN C SET VN = V = NORMALIZED VNP1 DO 55 I=1,N VN(I) = V(I) 55 CONTINUE 60 CONTINUE PRINT 65 65 FORMAT (' ***** INVERSE POWER METHOD DOES NOT CONVERGE *****') RETURN END SUBROUTINE CLINEQ(A,N,X,B,IPERM) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DECLARATIONS FOR ARGUMENTS COMPLEX*16 A(N,N),X(N),B(N) INTEGER N,IPERM(N) C DECLARATIONS FOR LOCAL VARIABLES COMPLEX*16 B_(N),LJI,TEMP,SUM C C SUBROUTINE CLINEQ SOLVES THE COMPLEX LINEAR SYSTEM A*X=B C C ARGUMENTS C C ON INPUT ON OUTPUT C -------- --------- C C A - THE N BY N COEFFICIENT MATRIX. THE DIAGONAL AND UPPER C TRIANGLE OF A CONTAINS U C AND THE LOWER TRIANGLE C OF A CONTAINS THE LOWER C TRIANGLE OF L, WHERE C PA = LU, P BEING THE C PERMUTATION MATRIX C DEFINED BY IPERM. C C N - THE SIZE OF MATRIX A. C C X - AN N-VECTOR CONTAINING C THE SOLUTION. C C B - THE RIGHT HAND SIDE N-VECTOR. C C IPERM - AN N-VECTOR CONTAINING C A RECORD OF THE ROW C INTERCHANGES MADE. IF C J = IPERM(K), THEN ROW C J ENDED UP AS THE K-TH C ROW. C C----------------------------------------------------------------------- DO 10 K=1,N C INITIALIZE IPERM = (1,2,3,...,N) IPERM(K) = K C COPY B TO B_, SO B WILL NOT BE ALTERED B_(K) = B(K) 10 CONTINUE C BEGIN FORWARD ELIMINATION DO 35 I=1,N-1 C SEARCH FROM A(I,I) ON DOWN FOR C LARGEST POTENTIAL PIVOT, A(L,I) BIG = ABS(A(I,I)) L = I DO 15 J=I+1,N IF (ABS(A(J,I)).GT.BIG) THEN BIG = ABS(A(J,I)) L = J ENDIF 15 CONTINUE C IF LARGEST POTENTIAL PIVOT IS ZERO, C MATRIX IS SINGULAR IF (BIG.EQ.0.0) GO TO 50 C SWITCH ROW I WITH ROW L, TO BRING C UP LARGEST PIVOT DO 20 K=1,N TEMP = A(L,K) A(L,K) = A(I,K) A(I,K) = TEMP 20 CONTINUE C SWITCH B_(I) AND B_(L) TEMP = B_(L) B_(L) = B_(I) B_(I) = TEMP C SWITCH IPERM(I) AND IPERM(L) ITEMP = IPERM(L) IPERM(L) = IPERM(I) IPERM(I) = ITEMP DO 30 J=I+1,N C CHOOSE MULTIPLIER TO ZERO A(J,I) LJI = A(J,I)/A(I,I) IF (LJI.NE.0.0) THEN C SUBTRACT LJI TIMES ROW I FROM ROW J DO 25 K=I+1,N A(J,K) = A(J,K) - LJI*A(I,K) 25 CONTINUE C SUBTRACT LJI TIMES B_(I) FROM B_(J) B_(J) = B_(J) - LJI*B_(I) ENDIF C SAVE LJI IN A(J,I). IT IS UNDERSTOOD, C HOWEVER, THAT A(J,I) IS REALLY ZERO. 30 A(J,I) = LJI 35 CONTINUE IF (A(N,N).EQ.0.0) GO TO 50 C SOLVE U*X = B_ USING BACK SUBSTITUTION. X(N) = B_(N)/A(N,N) DO 45 I=N-1,1,-1 SUM = 0.0 DO 40 J=I+1,N SUM = SUM + A(I,J)*X(J) 40 CONTINUE X(I) = (B_(I)-SUM)/A(I,I) 45 CONTINUE RETURN C MATRIX IS NUMERICALLY SINGULAR. 50 PRINT 55 55 FORMAT (' ***** THE MATRIX IS SINGULAR *****') RETURN END SUBROUTINE CRESLV(A,N,X,C,IPERM) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DECLARATIONS FOR ARGUMENTS COMPLEX*16 A(N,N),X(N),C(N) INTEGER N,IPERM(N) C DECLARATIONS FOR LOCAL VARIABLES COMPLEX*16 C_(N),LJI,SUM C C SUBROUTINE CRESLV SOLVES THE COMPLEX LINEAR SYSTEM A*X=C IN O(N**2) C TIME, AFTER CLINEQ HAS PRODUCED AN LU DECOMPOSITION OF PA. C C ARGUMENTS C C ON INPUT ON OUTPUT C -------- --------- C C A - THE N BY N COEFFICIENT MATRIX C AFTER PROCESSING BY CLINEQ. C AS OUTPUT BY CLINEQ, A CONTAINS C AN LU DECOMPOSITION OF PA. C C N - THE SIZE OF MATRIX A. C C X - AN N-VECTOR CONTAINING C THE SOLUTION. C C C - THE RIGHT HAND SIDE N-VECTOR. C C IPERM - THE PERMUTATION VECTOR OF C LENGTH N OUTPUT BY CLINEQ. C C----------------------------------------------------------------------- C CALCULATE C_=P*C, WHERE P IS PERMUTATION C MATRIX DEFINED BY IPERM. DO 10 K=1,N J = IPERM(K) C_(K) = C(J) 10 CONTINUE C BEGIN FORWARD ELIMINATION, TO CALCULATE C C_ = L**(-1)*C_ DO 20 I=1,N-1 DO 15 J=I+1,N C RETRIEVE MULTIPLIER SAVED IN A(J,I) LJI = A(J,I) C SUBTRACT LJI TIMES C_(I) FROM C_(J) C_(J) = C_(J) - LJI*C_(I) 15 CONTINUE 20 CONTINUE C SOLVE U*X = C_ USING BACK SUBSTITUTION. X(N) = C_(N)/A(N,N) DO 30 I=N-1,1,-1 SUM = 0.0 DO 25 J=I+1,N SUM = SUM + A(I,J)*X(J) 25 CONTINUE X(I) = (C_(I)-SUM)/A(I,I) 30 CONTINUE RETURN END