PROGRAM POWRMETHOD C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 11.1 (Power Method). C Section 11.2, The Power Method, Page 557 C C Algorithm 11.2 (Shifted Inverse Power Method). C Section 11.2, The Power Method, Page 558 PARAMETER(MaxR = 5) CHARACTER Ans*1,Ach*1,Mess*80 INTEGER Count,CountR,CountS INTEGER I,InRC,Inum,N,P,Q REAL Alpha,Apq,C,Det,Epsilon REAL Err,Lambda,MaxA,S,Rnum,T INTEGER Row,MaxIter DIMENSION Row(1:MaxR) REAL X,XP,XQ,YP,YQ DIMENSION X(1:MaxR),XP(1:MaxR),XQ(1:MaxR) DIMENSION YP(1:MaxR),YQ(1:MaxR) REAL A,A1,V DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) DIMENSION V(1:MaxR,1:MaxR) INTEGER Meth,Power,ShiftedInverse INTEGER Proc,Auto,Observe EXTERNAL XMaxEle,DIST Auto = 0 Observe = 1 Power = 1 ShiftedInverse = 2 CALL MESSAGE(Meth) CALL INPUTS(A,A1,N,InRC,MaxIter) CALL PROCESSES(Proc) IF (Meth.EQ.Power) THEN CALL POWMETH(A,N,Epsilon,MaxIter,X,Lambda,Err,Count,Proc) ENDIF IF (Meth.EQ.ShiftedInverse) THEN WRITE(9,*)' ' WRITE(9,*)'Enter the value Alpha = ' READ(9,*) Alpha WRITE(9,*)' ' CALL SHIFTEDINVPOWER(A,N,Epsilon,MaxIter, +Alpha,X,Lambda,Err,Count,Proc) ENDIF CALL PowrResults(A1,N,Epsilon,X,Lambda,Err,Count) END SUBROUTINE MATPRINT(A,N) PARAMETER(MaxR = 5) INTEGER Col,N,Row REAL A DIMENSION A(1:MaxR,1:MaxR) DO Row=1,N WRITE(9,999) ( A(Row,Col), Col=1,N ) ENDDO RETURN 999 FORMAT(1X,5F15.7,4X) END SUBROUTINE Voutput(X,N) PARAMETER(MaxR = 5) INTEGER N,Row REAL X DIMENSION X(1:MaxR) WRITE(9,*) DO Row=1,N WRITE(9,*) 'X(',Row,') = ',X(Row) ENDDO RETURN END REAL FUNCTION XMaxEle(X,N) PARAMETER(MaxR = 5) REAL X,MaxE DIMENSION X(1:MaxR) INTEGER J,N MaxE = 0 DO J=1,N IF ABS(X(J)) > Abs(MaxE) THEN MaxE = X(J) ENDIF ENDDO XMaxEle = MaxE RETURN END REAL FUNCTION DIST(X,Y,N) PARAMETER(MaxR = 5) REAL X,Y,Sum DIMENSION X(1:MaxR),Y(1:MaxR) INTEGER J,N Sum = 0 DO J=1,N Sum = Sum + (Y(J) - X(J))*(Y(J) - X(J)) ENDDO DIST = SQRT(Sum) RETURN END SUBROUTINE POWMETH(A,N,Epsilon,MaxIter,X,Lambda,Err,Count,Proc) PARAMETER(MaxR = 5) REAL A,X,Y,Epsilon,Lambda,Err,C1,DC,DV,Sum DIMENSION X(1:MaxR),Y(1:MaxR) DIMENSION A(1:MaxR,1:MaxR) INTEGER Count,I,J,N,MaxIter INTEGER State,Done,Iterating,Proc,Observe EXTERNAL XMaxEle,DIST DO J=1,N X(J) = 1 ENDDO Lambda = 0 Count = 0 Err = 1 Done = 0 Iterating = 1 State = Iterating Observe = 1 WHILE ((Count.LE.MaxIter).AND.(State.EQ.Iterating)) DO I=1,N !Matrix multiplication Y = AX} Sum = 0 DO J=1,N Sum = Sum + A(I,J)*X(J) ENDDO Y(I) = Sum ENDDO C1 = 0 C1 = XMaxEle(Y,N) !Find largest element of Y DC = ABS(Lambda - C1) DO J=1,N !Perform scalar multiplication Y(J) = (1/C1)*Y(J) ENDDO DV = DIST(X,Y,N) Err = MAX(DC,DV) DO J=1,N !Update vector X X(J) = Y(J) ENDDO Lambda = C1 !Update scalar Lambda State = Done IF (Err.GT.Epsilon) THEN !Check for convergence State = Iterating ENDIF Count = Count+1 IF (Proc.EQ.Observe) THEN WRITE(9,*)' ' WRITE(9,*)'Lambda = ', Lambda CALL Voutput(X,N) ENDIF REPEAT RETURN !OUTPUT vector X and scalar Lambda END SUBROUTINE FACTOR(A,Row,N,Det) PARAMETER(MaxR = 5) INTEGER Row,N DIMENSION Row(1:MaxR) REAL A,Det DIMENSION A(1:MaxR,1:MaxR) INTEGER C,J,K,P,RowK,RowP,T Det = 1 DO J=1,N !Initialize Pointer Vector Row(J) = J ENDDO DO P=1,(N-1) !Upper Triangularization Loop DO K=(P+1),N IF ABS(A(Row(K),P)) > ABS(A(Row(P),P)) THEN T = Row(P) Row(P) = Row(K) Row(K) = T Det = -Det ENDIF ENDDO Det = Det*A(Row(P),P) IF (Det.EQ.0) RETURN !Check Singular Matrix DO K=(P+1),N !Gaussian Elimination RowK = Row(K) RowP = Row(P) A(RowK,P) = A(RowK,P)/A(RowP,P) DO C=(P+1),N A(RowK,C) = A(RowK,C) - A(RowK,P)*A(RowP,C) ENDDO ENDDO ENDDO Det = Det*A(Row(N),N) RETURN END SUBROUTINE SOLVE(A,Row,B,X,N) PARAMETER(MaxR = 5) INTEGER Row,N DIMENSION Row(1:MaxR) REAL A,B,X DIMENSION B(1:MaxR),X(1:MaxR) DIMENSION A(1:MaxR,1:MaxR) INTEGER C,K,RowK REAL Sum DO K=1,N IF (A(Row(K),K).EQ.0) RETURN !Check Singular Matrix ENDDO X(1) = B(Row(1)) !Forward Substitution DO K=2,N Sum = 0; RowK = Row(K) DO C=1,(K-1) Sum = Sum + A(RowK,C)*X(C) ENDDO X(K) = B(RowK) - Sum !End Forward Substitution ENDDO X(N) = X(N)/A(Row(N),N) !Back Substitution DO K=(N-1),1,-1 Sum = 0 RowK = Row(K) DO C=(K+1),N Sum = Sum + A(RowK,C)*X(C) ENDDO X(K) = (X(K) - Sum)/A(RowK,K) ENDDO RETURN END SUBROUTINE SHIFTEDINVPOWER(A,N,Epsilon,MaxIter,Alpha,X, +Lambda,Err,Count,Proc) PARAMETER(MaxR = 5) REAL A,A1,X,Y DIMENSION X(1:MaxR),Y(1:MaxR) DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) INTEGER Row,Count,J,K,N,MaxIter DIMENSION Row(1:MaxR) INTEGER State,Done,Iterating,Proc,Observe REAL Epsilon,Alpha,Lambda,Err REAL C1,DC,DV,Det EXTERNAL XMaxEle,DIST Done = 0 Iterating = 1 Observe = 1 DO J=1,N !Initialize the matrix V X(J) = 1 ENDDO DO J=1,N !Form the matrix A - aIpha I DO K=1,N A1(J,K) = A(J,K) ENDDO ENDDO DO J=1,N A1(J,J) = A1(J,J) - Alpha ENDDO Lambda = 0 Count = 0 Err = 1 CALL FACTOR(A1,Row,N,Det) State = Iterating WHILE ((Count.LE.MaxIter).AND.(State.EQ.Iterating)) CALL SOLVE(A1,Row,X,Y,N) !Solve system A1 Y = X C1 = XMaxEle(Y,N) !Find largest element of Y DC = ABS(Lambda - C1) DO J=1,N !Perform scalar multiplication Y(J) = (1/C1)*Y(J) ENDDO DV = DIST(X,Y,N); Err = MAX(DC,DV) DO J=1,N !Update vector X X(J) = Y(J) ENDDO Lambda = C1 !Update scalar Lambda State = Done IF (Err.GT.Epsilon) THEN !Check for convergence State = Iterating ENDIF Count = Count+1 IF (Proc.EQ.Observe) THEN WRITE(9,*)' ' Lambda = Alpha+1/C1 WRITE(9,*)'Alpha = ',Alpha, ' ','C',Count,' = ',C1 WRITE(9,*)'Lambda = ', Lambda CALL Voutput(X,N); ENDIF REPEAT Lambda = Alpha+1/C1 !OUTPUT vector X and scalar Lambda RETURN END SUBROUTINE MATIN(Ach,A,N,InRC) PARAMETER(MaxR = 5) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaxR) DO 20 Row=1,N DO 10 Col=1,N A(Row,Col)=0 10 CONTINUE 20 CONTINUE WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach WRITE(9,*)' ' IF (InRC.EQ.0) THEN DO 40 Row=1,N WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row WRITE(9,*)' ' DO 30 Col=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 30 CONTINUE 40 CONTINUE ELSE DO 60 Col=1,N WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col WRITE(9,*)' ' DO 50 Row=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 50 CONTINUE 60 CONTINUE ENDIF RETURN END SUBROUTINE PowrResults(A1,N,Epsilon,X,Lambda,Err,Count) PARAMETER(MaxR = 5) REAL A1,X DIMENSION X(1:MaxR) DIMENSION A1(1:MaxR,1:MaxR) INTEGER C,R,N,Count REAL A0,Epsilon,Lambda,Err WRITE(9,*)' ' WRITE(9,*)'The matrix A is:' WRITE(9,*)' ' DO 10 R=1,N WRITE(9,999) ( A1(R,C), C=1,N ) ENDDO WRITE(9,*)' ' WRITE(9,*)'The dominant eigenvalue of A is:' WRITE(9,*)' ' WRITE(9,*)'Lambda = ', Lambda WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'The dominant eigenvector of A is:' WRITE(9,*)' ' CALL Voutput(X,N) WRITE(9,*)' ' WRITE(9,*)' ' IF (Count.EQ.Max) THEN WRITE(9,*)'The maximum number of iterations was used.' WRITE(9,*)' ' ENDIF PAUSE RETURN 999 FORMAT(1X,5F15.7,4X) END SUBROUTINE CHOOSEMETH(Meth) INTEGER I,Meth,Power,ShiftedInverse Power = 1 ShiftedInverse = 2 WRITE(9,*)' ' WRITE(9,*)'Choose a strategy for finding an eigen-pair.' WRITE(9,*)' ' WRITE(9,*)'<1> The power method.' WRITE(9,*)' ' WRITE(9,*)'<2> The shifted inverse power method.' WRITE(9,*)' ' WRITE(9,*)' SELECT the strategy <1 or 2> ? ' READ(9,*) I IF ((I.LT.1).OR.(2.LT.I)) I = 1 IF (I.EQ.1) Meth = Power IF (I.EQ.2) Meth = ShiftedInverse RETURN END SUBROUTINE MESSAGE(Meth) INTEGER Meth WRITE(9,*)' POWER METHOD FOR FINDING AN EIGEN-PAIR' WRITE(9,*)' ' WRITE(9,*)'Assume that A is an N by N real matrix and that it' WRITE(9,*)' ' WRITE(9,*)'has a full set of eigenvectors V , V ,..., V . ' WRITE(9,*)' 1 2 N ' WRITE(9,*)' ' WRITE(9,*)'The powermethod of iteration is used to find an' WRITE(9,*)' ' WRITE(9,*)'eigenvalue and its corresponding eigenvector.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Press the key. ' PAUSE CALL CHOOSEMETH(Meth) RETURN END SUBROUTINE INPUTS(A,A1,N,InRC,MaxIter) PARAMETER(MaxR = 5) REAL A,A1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) INTEGER N,InRC,MaxIter,Row,Col WRITE(9,*)' ' WRITE(9,*)'We will now proceed with the power method' WRITE(9,*)' ' WRITE(9,*)'A must be a matrix of dimension N by N.' WRITE(9,*)' ' WRITE(9,*)'{N must be an integer between 2 and 10}' WRITE(9,*)' ' WRITE(9,*)'ENTER N = ' READ(9,*) N IF (N.LT.2) N = 2 IF (N.GT.10) N = 10 WRITE(9,*)'Iteration will continue until each coordinate' WRITE(9,*)' ' WRITE(9,*)'of the eigenvector has converged with an ' WRITE(9,*)' ' WRITE(9,*)'error less than Epsilon or the maximum' WRITE(9,*)' ' WRITE(9,*)'number of iterations Max has been reached.' WRITE(9,*)' ' WRITE(9,*)'ENTER Epsilon = ' READ(9,*) Epsilon IF (Epsilon.LT. 0.000000001) Epsilon = 0.000000001 WRITE(9,*)' ' WRITE(9,*)'ENTER Max = ' READ(9,*) MaxIter IF (50.LT.MaxIter) MaxIter = 50 InRC = 0 WRITE(9,*)'WANT TO INPUT COLUMNS ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN InRC = 1 ENDIF WRITE(9,*)' ' WRITE(9,*)'ENTER THE COEFFICIENT MATRIX A:' CALL MATIN('A',A,N,InRC) DO Row=1,N DO Col=1,N A1(Row,Col)=A(Row,Col) ENDDO ENDDO RETURN END SUBROUTINE PROCESSES(Proc) INTEGER I,Proc,Auto,Observe Auto = 0 Observe = 1 WRITE(9,*)' ' WRITE(9,*)'To what extent do you want to control the program?' WRITE(9,*)' ' WRITE(9,*)'<1> The computer does it all automatically.' WRITE(9,*)' ' WRITE(9,*)'<2> Computer selects, but we observe each step.' WRITE(9,*)' ' WRITE(9,*)'SELECT your choice for input <1 - 2> ? ' READ(9,*) I IF ((I.LT.1).OR.(2.LT.I)) I = 2 IF (I.EQ.1) Proc = Auto IF (I.EQ.2) Proc = Observe RETURN END PROGRAM JACOBIE C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 11.3 (Jacobi Iteration for Eigenvalues and Eigenvectors). C Section 11.3, Jacobi's Method, Page 571 C PARAMETER(MaxR = 4, MaxS = 50) CHARACTER Ans*1,Ach*1,Mess*80 INTEGER CountR,CountS INTEGER I,InRC,Inum,N,P,Q,Row REAL Apq,C,Epsilon,MaxA,S,Rnum,T DIMENSION Row(1:MaxR) REAL XP,XQ,YP,YQ DIMENSION XP(1:MaxR),XQ(1:MaxR) DIMENSION YP(1:MaxR),YQ(1:MaxR) REAL A,A1,V DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) DIMENSION V(1:MaxR,1:MaxR) INTEGER Method,Cyclic,Original INTEGER Proc,Auto,Manual,Observe EXTERNAL RMS Original = 1 Cyclic = 2 Auto = 1 Manual = 2 Observe = 3 CALL MESSAGE(Meth) CALL INPUTS(A,A1,N,Epsilon,InRC) CALL PROCESSES(Proc) CALL INITIALV(V,N) CALL MAXPQ(A,N,P,Q,MaxA) IF (Meth.EQ.Original) THEN CALL OJacob(A,V,N,Epsilon,CountR,Proc) ENDIF IF (Meth.EQ.Cyclic) THEN CALL CJacob(A,V,N,Epsilon,CountS,Proc) ENDIF CALL PrintResults(A,A1,V,N,Epsilon,MaxA,CountS,Meth) END SUBROUTINE Routput(P,Q,C,S,Apq) INTEGER P,Q REAL Apq,C,S WRITE(9,*)' ' WRITE(9,*)'The rotation matrix for zeroing out' WRITE(9,*)'A(',P, ',',Q, ') =',Apq,' is' WRITE(9,*)' ' WRITE(9,*)' ( ) ( )' WRITE(9,*)' ( c s ) (', C, ' ', S, ' )' WRITE(9,*)'R = ( ) = ( )' WRITE(9,*)' ( -s c ) (', -S, ' ', C , ' )' WRITE(9,*)' ( ) ( )' RETURN END SUBROUTINE MATPRINT(A,N) PARAMETER(MaxR = 4) INTEGER Col,N,Row REAL A DIMENSION A(1:MaxR,1:MaxR) WRITE(9,*)' ' DO 10 Row=1,N WRITE(9,999) ( A(Row,Col), Col=1,N ) 10 CONTINUE 999 FORMAT(1X,5F15.7,4X) RETURN END SUBROUTINE ZEROPQ(A,V,C,S,N,P,Q) PARAMETER(MaxR = 4) REAL A,V DIMENSION A(1:MaxR,1:MaxR),V(1:MaxR,1:MaxR) INTEGER I,J,N,P,Q,Pm,Qm REAL C,S,T,Theta REAL XP,XQ,YP,YQ DIMENSION XP(1:MaxR),XQ(1:MaxR) DIMENSION YP(1:MaxR),YQ(1:MaxR) Theta = (A(Q,Q) - A(P,P))/(2*A(P,Q)) IF (Theta.LT.1E19) THEN T = 1/(ABS(Theta) + SQRT(Theta*Theta+1)) ELSE T = 1/ABS(2*Theta) ENDIF IF (Theta.LT.0) THEN T = -T ENDIF C = 1/SQRT(T*T+1) S = C*T DO I=1,N XP(I) = A(I,P) XQ(I) = A(I,Q) ENDDO A(P,Q) = 0 A(Q,P) = 0 A(P,P) = C*C*XP(P) + S*S*XQ(Q) - 2*C*S*XQ(P) A(Q,Q) = S*S*XP(P) + C*C*XQ(Q) + 2*C*S*XQ(P) DO I=1,N IF ((I.NE.P).AND.(I.NE.Q)) THEN A(I,P) = C*XP(I) - S*XQ(I) A(I,Q) = C*XQ(I) + S*XP(I) A(P,I) = A(I,P) A(Q,I) = A(I,Q) ENDIF ENDDO DO I=1,N YP(I) = V(I,P) YQ(I) = V(I,Q) ENDDO DO I=1,N V(I,P) = C*YP(I) - S*YQ(I) V(I,Q) = S*YP(I) + C*YQ(I) ENDDO RETURN END SUBROUTINE MAXPQ(A,N,P,Q,MaxA) PARAMETER(MaxR = 4) REAL A,MaxA DIMENSION A(1:MaxR,1:MaxR) INTEGER I,J,N,P,Q,Pm,Qm MaxA = ABS(A(1,N)) Pm = 1 Qm = N DO I=1,N DO J=(I+1),N IF (ABS(A(I,J)).GT.MaxA) THEN Pm = I Qm = J MaxA = ABS(A(I,J)) ENDIF ENDDO ENDDO P = Pm Q = Qm RETURN END SUBROUTINE INITIALV(V,N) PARAMETER(MaxR = 4) REAL V DIMENSION V(1:MaxR,1:MaxR) INTEGER I,J,N DO I=1,N DO J=1,N IF (I.NE.J) THEN V(I,J) = 0 ELSE V(I,J) = 1 ENDIF ENDDO ENDDO RETURN END REAL FUNCTION RMS(A,N) PARAMETER(MaxR = 4) REAL A DIMENSION A(1:MaxR,1:MaxR) INTEGER J,N REAL Sum Sum = 0 DO J=1,N Sum = Sum + A(J,J) ENDDO RMS = SQRT(Sum/N) RETURN END SUBROUTINE SELECTPQ(A,N,P,Q,MaxA,Proc) PARAMETER(MaxR = 4) REAL A,Apq,MaxA DIMENSION A(1:MaxR,1:MaxR) INTEGER C,I,J,N,P,Q,Ptemp,Qtemp,Proc CHARACTER Resp*1 Manual = 2 CALL MAXPQ(A,N,P,Q,MaxA) Apq = A(P,Q) IF (Proc.NE.Manual) RETURN CALL MATPRINT(A,N) WRITE(9,*)' ' WRITE(9,*)' Choose an off diagonal element to "zero-out."' WRITE(9,*)' ' IF (N.EQ.2) THEN WRITE(9,*)'Select the row P = 1,2' WRITE(9,*)' and column Q = 1,2' ELSEIF (N.EQ.3) THEN WRITE(9,*)'Select the row P = 1,2,3' WRITE(9,*)' and column Q = 1,2,3' ELSE WRITE(9,*)'Select the row P = 1,2,...,', N WRITE(9,*)' and column Q = 1,2,...,', N ENDIF WRITE(9,*)' ' WRITE(9,*)' ENTER the row P = ' READ(9,*) Ptemp IF ((1.LE.Ptemp).AND.(Ptemp.LE.N)) P = Ptemp WRITE(9,*)'N = ',N WRITE(9,*)'Ptemp = ',Ptemp WRITE(9,*)'P = ',P WRITE(9,*)' ENTER column Q = ' READ(9,*) Qtemp IF ((1.LE.Qtemp).AND.(Qtemp.LE.N)) Q = Qtemp WRITE(9,*)'N = ',N WRITE(9,*)'Qtemp = ',Qtemp WRITE(9,*)'Q = ',Q Apq = A(P,Q) WRITE(9,*)' ' WRITE(9,*)'The element selected was A(',P,',',Q,') =',A(P, Q) WRITE(9,*)' ' PAUSE RETURN END SUBROUTINE OJacob(A,V,N,Epsilon,CountR,Proc) PARAMETER(MaxR = 4) REAL A,V,C,S,MaxA,Apq,Epsilon DIMENSION A(1:MaxR,1:MaxR),V(1:MaxR,1:MaxR) INTEGER CountR,N,P,Q,Stat,Done,Working INTEGER Proc,Auto,Manual,Observe CHARACTER Ans*1 Done = 0 Working = 1 Stat = Working Auto = 1 Manual = 2 Observe = 3 CountR = 0 WHILE(Stat.EQ.Working) CALL SELECTPQ(A,N,P,Q,MaxA,Proc) Apq = MaxA CALL ZEROPQ(A,V,C,S,N,P,Q) CountR = CountR+1 IF (Proc.EQ.Auto) THEN WRITE(9,*)'I am computing rotation # ', CountR ENDIF IF ((Proc.EQ.Manual).OR.(Proc.EQ.Observe)) THEN CALL Routput(P,Q,C,S,Apq) CALL MATPRINT(A,N) ENDIF CALL MAXPQ(A,N,P,Q,MaxA) IF (MaxA.LT.Epsilon*RMS(A,N)) Stat = Done IF (Proc.EQ.Manual) THEN WRITE(9,*)' ' WRITE(9,*)'Want to zero out another element ? ' READ(9,*) Ans IF ((Ans.EQ.'N').OR.(Ans.EQ.'n')) Stat = Done ENDIF REPEAT RETURN END SUBROUTINE CJacob(A,V,N,Epsilon,CountS,Proc) PARAMETER(MaxR = 4, MaxS = 50) REAL A,V,C,S,Apq,Epsilon DIMENSION A(1:MaxR,1:MaxR),V(1:MaxR,1:MaxR) INTEGER CountR,CountS,N,P,Q,Stat,Done,Working INTEGER Proc,Auto,Manual,Observe EXTERNAL RMS CountR = 0 CountS = 0 Done = 0 Working = 1 Stat = Working Auto = 1 Manual = 2 Observe = 3 WHILE ((Stat.EQ.Working).AND.(CountS.LT.MaxS)) CountS = CountS+1 T = RMS(A,N) Stat = Done DO P=1,(N-1) DO Q=(P+1),N IF ((ABS(A(P,Q))/T).GT.Epsilon) THEN Apq = A(P,Q) CALL ZEROPQ(A,V,C,S,N,P,Q) CountR = CountR+1 IF (Proc.EQ.Auto) THEN WRITE(9,*)'I am computing sweep # ', + CountS,' , rotation # ',CountR ENDIF IF (Proc.EQ.Observe) THEN CALL Routput(P,Q,C,S,Apq) CALL MATPRINT(A,N) ENDIF Stat = Working ENDIF ENDDO ENDDO REPEAT RETURN END SUBROUTINE MATIN(Ach,A,N,InRC) PARAMETER(MaxR = 4) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaxR) DO 20 Row=1,N DO 10 Col=1,N A(Row,Col)=0 10 CONTINUE 20 CONTINUE WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach WRITE(9,*)' ' IF (InRC.EQ.0) THEN DO 40 Row=1,N WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row WRITE(9,*)' ' DO 30 Col=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 30 CONTINUE 40 CONTINUE ELSE DO 60 Col=1,N WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col WRITE(9,*)' ' DO 50 Row=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 50 CONTINUE 60 CONTINUE ENDIF RETURN END SUBROUTINE PrintResults(A,A1,V,N,Epsilon,MaxA,CountS,Meth) PARAMETER(MaxR = 4, MaxS = 50) REAL A0,A1,V,Epsilon,MaxA DIMENSION A(1:MaxR,1:MaxR), A1(1:MaxR,1:MaxR),V(1:MaxR,1:MaxR) INTEGER N,C,R,P,Q,CountS INTEGER I,Meth,Cyclic,Original EXTERNAL RMS Original = 1 Cyclic = 2 WRITE(9,*)' ' WRITE(9,*)'The matrix A is:' WRITE(9,*)' ' DO R=1,N WRITE(9,999) ( A1(R,C), C=1,N ) ENDDO WRITE(9,*)' ' WRITE(9,*)'The eigenvalues of A are:' DO R=1,N WRITE(9,*)'X(',R,') = ',A(R,R) ENDDO WRITE(9,*)' ' WRITE(9,*)'The matrix of eigenvectors V is:' DO R=1,N WRITE(9,999) ( V(R,C), C=1,N ) ENDDO WRITE(9,*)' ' CALL MAXPQ(A,N,P,Q,MaxA) IF ((MaxA.GT.RMS(A,N)).AND.(Meth.EQ.Original)) THEN WRITE(9,*)' ' WRITE(9,*)'The iteration did NOT converge!' ENDIF IF ((MaxS.LT.CountS).AND.(Meth.EQ.Cyclic)) THEN WRITE(9,*)' ' WRITE(9,*)'The iteration did NOT converge!' ENDIF PAUSE RETURN 999 FORMAT(1X,5F15.7,4X) END SUBROUTINE CHOOSEMETH(Meth) INTEGER I,Meth,Cyclic,Original Original = 1 Cyclic = 2 WRITE(9,*)'Choose a strategy for selecting the off diagonal' WRITE(9,*)' ' WRITE(9,*)'element to annihilate in the construction process.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'<1> Jacobi`s original strategy.' WRITE(9,*)' ' WRITE(9,*)'Select p,q so that |a | = max |a |.' WRITE(9,*)' pq i The cyclic Jacobi method.' WRITE(9,*)' ' WRITE(9,*)'Now sweep through the matrix and annihilate' WRITE(9,*)' ' WRITE(9,*)'elements in the strict order' WRITE(9,*)' ' WRITE(9,*)'a ,a ,...,a ;a ,a ,...,a ;...;a .' WRITE(9,*)' 12 13 1N 23 24 2N N-1,N ' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' SELECT the strategy < 1 or 2 > ? ' READ(9,*) I IF ((I.LT.1).OR.(2.LT.I)) I = 1 IF (I.EQ.1) Meth = Original IF (I.EQ.2) Meth = Cyclic RETURN END SUBROUTINE MESSAGE(Meth) INTEGER I,Meth CHARACTER Ans*1 WRITE(9,*)' JACOBI`S METHOD FOR EIGENVECTORS' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' Assume that A is an N by N real symmetric matrix.' WRITE(9,*)' ' WRITE(9,*)'Then A has a full set of eigenvectors' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'V , V ,..., V . Jacobi`s method of iteration is used' WRITE(9,*)' 1 2 N ' WRITE(9,*)' ' WRITE(9,*)'to find the eigenvalues and eigenvectors of A. Let' WRITE(9,*)' ' WRITE(9,*)'A = A , and construct a sequence of orthogonal' WRITE(9,*)' 1 ' WRITE(9,*)' T ' WRITE(9,*)'matrices { R } such that D = R A R .' WRITE(9,*)' j j j j j ' WRITE(9,*)' ' WRITE(9,*)'{ D } converges to the diagonal matrix D ,' WRITE(9,*)' j ' WRITE(9,*)' ' WRITE(9,*)'of eigenvalues, and the sequence {V = R R ...R }' WRITE(9,*)' j 1 2 j ' WRITE(9,*)' ' WRITE(9,*)'converges to the matrix of eigenvectors.' WRITE(9,*)' ' WRITE(9,*)' Press the key. ' READ(9,*) Ans CALL CHOOSEMETH(Meth) RETURN END SUBROUTINE INPUTS(A,A1,N,Epsilon,InRC) PARAMETER(MaxR = 4) REAL A,A1,Epsilon DIMENSION A(1:MaxR,1:MaxR), A1(1:MaxR,1:MaxR) INTEGER InRC,N,C,Row,Col CHARACTER Ans*1 WRITE(9,*)' We will now proceed with Jacobi`s method' WRITE(9,*)' ' WRITE(9,*)'of iteration to find the set of N eigenvectors' WRITE(9,*)' ' WRITE(9,*)'V , V ,..., V for the symmetric matrix A.' WRITE(9,*)' 1 2 N ' WRITE(9,*)' ' WRITE(9,*)'A must be a symmetric matrix of dimension N by N.' WRITE(9,*)' ' WRITE(9,*)'{N must be an integer between 2 and 10}' WRITE(9,*)' ' WRITE(9,*)' ENTER N = ' READ(9,*) N IF (N.LT.2) N = 2 IF (N.GT.10) N = 10 WRITE(9,*)' ' WRITE(9,*)' For each sweep throughout the matrix,' WRITE(9,*)' ' WRITE(9,*)'the computed value T is the [R.M.S.] average' WRITE(9,*)' ' WRITE(9,*)'of the diagonal elements of the matrix A.' WRITE(9,*)' ' WRITE(9,*)'Give the error criterion for annihilating' WRITE(9,*)' ' WRITE(9,*)'the off diagonal elements' WRITE(9,*)' ' WRITE(9,*)'|d | > T*Epsilon for all q > p.' WRITE(9,*)' p,q' WRITE(9,*)' ' WRITE(9,*)' ENTER Epsilon = ' READ(9,*) Epsilon IF (Epsilon.LT.0.000000001) Epsilon = 0.000000001 InRC = 0 WRITE(9,*)'WANT TO INPUT COLUMNS ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN InRC = 1 ENDIF WRITE(9,*)' ' WRITE(9,*)'ENTER THE COEFFICIENT MATRIX A:' CALL MATIN('A',A,N,InRC) DO 30 Row=1,N DO 20 Col=1,N A1(Row,Col)=A(Row,Col) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE PROCESSES(Proc) INTEGER I,Proc,Auto,Manual,Observe Auto = 1 Manual = 2 Observe = 3 WRITE(9,*)'To what extent do you want to control the program?' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'<1> The computer does it all automatically.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'<2> The user selects elements to zero-out manually.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'<3> Computer selects, but we observe each step.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'SELECT your choice for input < 1 - 3 > ? ' READ(9,*) I IF ((I.LT.1).OR.(3.LT.I)) I = 3 IF (I.EQ.1) Proc = Auto IF (I.EQ.2) Proc = Manual IF (I.EQ.3) Proc = Observe RETURN END PROGRAM HOUSE C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 11.4 (Reduction to Tridiagonal Form). C Section 11.4, Eigenvalues of Symmetric Matrices, Page 581 C PARAMETER(MaxR = 5, MaxS = 50) CHARACTER Ans*1,Ach*1,Mess*80 INTEGER CountR,CountS INTEGER I,InRC,Inum,K,N REAL C,MaxA,Rnum,R,S,T REAL W,V,Q DIMENSION W(1:MaxR),V(1:MaxR),Q(1:MaxR) REAL A,A1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) INTEGER MatType INTEGER Proc,Auto,Observe Auto = 1 Observe = 2 CALL MESSAGE(Meth) CALL INPUTS(A,A1,N,InRC) CALL PROCESSES(Proc) CALL HOUSEHOLDER(A,W,V,Q,S,R,C,N,Proc) CALL PrintResults(A,A1,N) END SUBROUTINE MATPRINT(A,K,N) PARAMETER(MaxR = 5) INTEGER Col,K,N,Row REAL A DIMENSION A(1:MaxR,1:MaxR) WRITE(9,*) 'The matrix A is:' WRITE(9,*) ' ',K DO 10 Row=1,N WRITE(9,999) ( A(Row,Col), Col=1,N ) 10 CONTINUE 999 FORMAT(1X,5F15.7,4X) RETURN END SUBROUTINE StepOutput(W,V,Q,S,R,C,N) PARAMETER(MaxR = 5, MaxS = 50) REAL W,V,Q,S,R,C DIMENSION W(1:MaxR),V(1:MaxR),Q(1:MaxR) INTEGER N,J WRITE(9,*)' ' WRITE(9,*)' Constant S Constant + R Constant C' WRITE(9,*)' ',S,' ', + R,' ', C WRITE(9,*)' ' WRITE(9,*)' Vector W Vector + V Vector Q' DO J=1,N WRITE(9,*)' ',W(J), ' ', + V(J), ' ', Q(J) ENDDO RETURN END SUBROUTINE Form_W(W,A,S,R,K,N) PARAMETER(MaxR = 5) REAL W,A,S,R,Sum DIMENSION W(1:MaxR) DIMENSION A(1:MaxR,1:MaxR) INTEGER J,K,N Sum = 0 DO J=(K+1),N Sum = Sum + A(J,K)**2 ENDDO S = SQRT(Sum) IF A(K+1,K) < 0 THEN S = -S ENDIF R = SQRT(2*(A(K+1,K) + S)*S) DO J=1,K !First k elements are 0 W(J) = 0; ENDDO W(K+1) = (A(K+1,K) + S)/R DO J=(K+2),N W(J) = A(J,K)/R ENDDO RETURN END SUBROUTINE Form_V(V,A,W,K,N) PARAMETER(MaxR = 5) REAL V,W DIMENSION V(1:MaxR),W(1:MaxR) REAL A DIMENSION A(1:MaxR,1:MaxR) INTEGER I,J,K,N REAL Sum DO I=1,K !First k elements are 0 V(I) = 0 ENDDO DO I=(K+1),N Sum = 0 DO J=(K+1),N Sum = Sum + A(I,J)*W(J) ENDDO V(I) = Sum ENDDO RETURN END SUBROUTINE Form_Q(Q,V,W,C,K,N) PARAMETER(MaxR = 5) REAL C,V,W,Q DIMENSION V(1:MaxR),W(1:MaxR),Q(1:MaxR) INTEGER J,K,N C = 0 DO J=(K+1),N C = C + W(J)*V(J) ENDDO DO J=1,K !First k elements are 0} Q(J) = 0 ENDDO DO J=(K+1),N Q(J) = V(J) - C*W(J) ENDDO RETURN END SUBROUTINE Form_A(A,Q,W,S,K,N) PARAMETER(MaxR = 5) REAL A,Q,W,S DIMENSION Q(1:MaxR),W(1:MaxR) DIMENSION A(1:MaxR,1:MaxR) INTEGER I,J,K,N DO J=(K+2),N A(J,K) = 0 A(K,J) = 0 ENDDO A(K+1,K) = -S A(K,K+1) = -S DO J=K,N A(J,J) = A(J,J) - 4*Q(J)*W(J) ENDDO DO I=(K+1),N DO J=(I+1),N A(I,J) = A(I,J) - 2*W(I)*Q(J) - 2*Q(I)*W(J) A(J,I) = A(I,J); ENDDO ENDDO RETURN END SUBROUTINE HOUSEHOLDER(A,W,V,Q,S,R,C,N,Proc) PARAMETER(MaxR = 5) REAL A,W,V,Q,S,R,C DIMENSION W(1:MaxR),V(1:MaxR),Q(1:MaxR) DIMENSION A(1:MaxR,1:MaxR) INTEGER K,N,Proc,Observe Observe = 2 DO K=1,(N-2) CALL Form_W(W,A,S,R,K,N) CALL Form_V(V,A,W,K,N) CALL Form_Q(Q,V,W,C,K,N) CALL Form_A(A,Q,W,S,K,N) IF (Proc.EQ.Observe) THEN CALL StepOutput(W,V,Q,S,R,C,N) CALL MATPRINT(A,K,N) WRITE(9,*)' ' WRITE(9,*)' Press the key. ' PAUSE ENDIF ENDDO RETURN END SUBROUTINE MATIN(Ach,A,N,InRC) PARAMETER(MaxR = 5) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaxR) DO 20 Row=1,N DO 10 Col=1,N A(Row,Col)=0 10 CONTINUE 20 CONTINUE WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach WRITE(9,*)' ' IF (InRC.EQ.0) THEN DO 40 Row=1,N WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row WRITE(9,*)' ' DO 30 Col=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 30 CONTINUE 40 CONTINUE ELSE DO 60 Col=1,N WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col WRITE(9,*)' ' DO 50 Row=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 50 CONTINUE 60 CONTINUE ENDIF RETURN END SUBROUTINE MESSAGE(Meth) INTEGER Meth WRITE(9,*)' HOUSEHOLDER`S METHOD' WRITE(9,*)' ' WRITE(9,*)' Assume that A is an N by N real symmetric matrix.' WRITE(9,*)' ' WRITE(9,*)'Then A is similar to a tridiagonal matrix.' WRITE(9,*)' ' WRITE(9,*)'Starting with A = A, Householder`s method will' WRITE(9,*)' 0 ' WRITE(9,*)' ' WRITE(9,*)'construct a sequence of orthogonal symmetric matrices' WRITE(9,*)' ' WRITE(9,*)'{P } such that A = P A P ,' WRITE(9,*)' k k k k-1 k ' WRITE(9,*)' ' WRITE(9,*)' for k = 1,2,...,N-2.' WRITE(9,*)' ' WRITE(9,*)'Then A is the desired tridiagonal matrix.' WRITE(9,*)' N-2 ' WRITE(9,*)' ' PAUSE RETURN END SUBROUTINE INPUTS(A,A1,N,InRC) PARAMETER(MaxR = 5) REAL A,A1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) INTEGER N,InRC,Row,Col WRITE(9,*)' We will now proceed with Householder`s' WRITE(9,*)' ' WRITE(9,*)'method of iteration to find the tridiagonal' WRITE(9,*)' ' WRITE(9,*)'matrix A that is similar to A.' WRITE(9,*)' N-2 ' WRITE(9,*)' ' WRITE(9,*)'A must be a symmetric N by N matrix.' WRITE(9,*)' ' WRITE(9,*)'{N must be an integer between 2 and 10}' WRITE(9,*)' ' WRITE(9,*)'ENTER N = ' READ(9,*) N IF (N.LT.2) N = 2 IF (N.GT.10) N = 10 InRC = 0 WRITE(9,*)'WANT TO INPUT COLUMNS ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN InRC = 1 ENDIF WRITE(9,*)' ' WRITE(9,*)'ENTER THE COEFFICIENT MATRIX A:' CALL MATIN('A',A,N,InRC) DO 30 Row=1,N DO 20 Col=1,N A1(Row,Col)=A(Row,Col) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE PROCESSES(Proc) INTEGER I,Proc,Auto,Observe Auto = 1 Observe = 2 WRITE(9,*)'To what extent do you want to control the program?' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'<1> The computer does it all automatically.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'<3> Computer selects, but we observe each step.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'SELECT your choice for input <1 or 2> ? ' READ(9,*) I IF ((I.LT.1).OR.(3.LT.I)) I = 3 IF (I.EQ.1) Proc = Auto IF (I.EQ.2) Proc = Observe RETURN END SUBROUTINE PrintResults(A,A1,N) PARAMETER(MaxR = 5) REAL A,A1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) INTEGER N,R,C WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'The matrix A is:' WRITE(9,*)' ' DO 10 R=1,N WRITE(9,999) ( A1(R,C), C=1,N ) ENDDO WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'The similar tridiagonal matrix is:' DO 10 R=1,N WRITE(9,999) ( A(R,C), C=1,N ) ENDDO PAUSE RETURN 999 FORMAT(1X,5F15.7,4X) END PROGRAM QLALGOR C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 11.5 (The QL Method with Shifts). C Section 11.4, Eigenvalues of Symmetric Matrices, Page 587 C PARAMETER(MaxR = 10,MaxS = 50) CHARACTER Ans*1,Ach*1,Mess*80 INTEGER CountR,CountS,I,InRC,Inum,K,N REAL MaxA,Rnum,T REAL C,D,E,P,R,S,Q DIMENSION C(1:MaxR),D(1:MaxR),E(1:MaxR) DIMENSION P(1:MaxR),R(1:MaxR),S(1:MaxR),Q(1:MaxR) DIMENSION YP(1:MaxR),YQ(1:MaxR) REAL A,A1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) CALL MESSAGE(InRC,Mtype) CALL INPUTS(A,A1,N,InRC) CALL PROCESSES CALL QL_ITERATION(A,D,N) CALL PrintResults(A,A1,D,N) END SUBROUTINE MATPRINT(A,K,N) PARAMETER(MaxR=10) INTEGER Col,K,N,Row REAL A DIMENSION A(1:MaxR,1:MaxR) WRITE(9,*) 'The matrix A is:' WRITE(9,*) ' ',K DO 10 Row=1,N WRITE(9,999) ( A(Row,Col), Col=1,N ) 10 CONTINUE 999 FORMAT(1X,5F15.7,4X) RETURN END SUBROUTINE A_TO_TRI(D,E,A,N) REAL D,E,A PARAMETER(MaxR = 10) DIMENSION D(1:MaxR),E(1:MaxR) DIMENSION A(1:MaxR,1:MaxR) INTEGER N,J DO J=1,N D(J) = A(J,J) ENDDO DO J=1,(N-1) E(J) = A(J,J+1) ENDDO RETURN END SUBROUTINE TRI_TO_A(D,E,A,N) REAL D,E,A PARAMETER(MaxR = 10) DIMENSION C(1:MaxR),D(1:MaxR),E(1:MaxR) DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) INTEGER N,J DO J=1,N A(J,J) = D(J) ENDDO DO J=1,(N-1) A(J,J+1) = E(J) A(J+1,J) = E(J) ENDDO RETURN END SUBROUTINE Find_Shift(D,E,SH,R1,R2,M) PARAMETER(MaxR = 10) REAL D,E,SH,B1,C1,D1,R1,R2 DIMENSION D(1:MaxR),E(1:MaxR) INTEGER M B1 = -(D(M+1) + D(M)) C1 = D(M+1)*D(M) - E(M)*E(M) D1 = SQRT(ABS(B1*B1 - 4*C1)) IF (B1.GT.0) THEN R1 = (-B1 - D1)/2; R2 = 2*C1/(-B1 - D1) ELSE R1 = (-B1 + D1)/2 R2 = 2*C1/(-B1 + D1) ENDIF SH = R1 IF (ABS(D(M) - R2).LT.ABS(D(M) - R1)) THEN SH = R2 ENDIF WRITE(9,*)' ' WRITE(9,*)'the shift is ',SH WRITE(9,*)' ' PAUSE RETURN END SUBROUTINE Form_L(C,D,E,P,S,Q,N,M) PARAMETER(MaxR = 10) REAL C,D,E,P,S,Q DIMENSION C(1:MaxR),D(1:MaxR),E(1:MaxR) DIMENSION P(1:MaxR),S(1:MaxR),Q(1:MaxR) INTEGER N,M,J REAL P_J1,Q_J DO J=1,N P(J) = D(J) ENDDO DO J=1,(N-1) Q(J) = E(J) ENDDO P_J1 = P(N) !P_J1 old value of P(J+1) Q_J = Q(N-1) !Q_J old value of Q(J) DO J=(N-1),M,-1 P(J+1) = SQRT(P_J1**2 + Q(J)**2) C(J) = P_J1/P(J+1) S(J) = Q(J)/P(J+1) Q(J) = C(J)*Q_J + S(J)*D(J) P_J1 = C(J)*D(J) - S(J)*Q_J IF (J.GT.M) THEN ! R(J-1) = S(J)*Q(J-1) dead code Q_J = C(J)*Q(J-1) ENDIF ENDDO P(M) = P_J1 RETURN END SUBROUTINE Form_A(C,D,E,P,S,Q,N,M) PARAMETER(MaxR = 10) REAL C,D,E,P,S,Q DIMENSION C(1:MaxR),D(1:MaxR),E(1:MaxR) DIMENSION P(1:MaxR),S(1:MaxR),Q(1:MaxR) INTEGER N,M,J D(N) = S(N-1)*Q(N-1) + C(N-1)*P(N) E(N-1) = S(N-1)*P(N-1) DO J=(N-2),M,-1 D(J+1) = S(J)*Q(J) + C(J)*C(J+1)*P(J+1) E(J) = S(J)*P(J) ENDDO D(M) = C(M)*P(M) RETURN END SUBROUTINE QL_ITERATION(A,D,N) PARAMETER(DELTA = 0.00000001) PARAMETER(MaxR = 10) REAL A,C,D,E,P,R,S,Q DIMENSION C(1:MaxR),D(1:MaxR),E(1:MaxR) DIMENSION P(1:MaxR),R(1:MaxR),S(1:MaxR),Q(1:MaxR) DIMENSION A(1:MaxR,1:MaxR) INTEGER Cond,J,M,N,Proc,Observe REAL SH,SHIFT,R1,R2 Observe = 1 IF (Proc.EQ.Observe) THEN CALL Aoutput(A,0,N) ENDIF CALL A_TO_TRI(D,E,A,N) SHIFT = 0 M = 1 WHILE (M.LE.(N-2)) K = 1 Cond = 0 WHILE ((K.LT.50).AND.(Cond.EQ.0)) CALL Find_Shift(D,E,SH,R1,R2,M) IF (ABS(E(M)).GT.DELTA) THEN SHIFT = SHIFT + SH DO J=M,N D(J) = D(J) - SH ENDDO ELSE WRITE(9,*)' ' WRITE(9,*)'the eigenvalue is ' D(M) = D(M) + SHIFT; WRITE(9,*) D(M) WRITE(9,*)' ' WRITE(9,*)' Press the key. ' PAUSE M = M+1 Cond = 1 ENDIF CALL Form_L(C,D,E,P,S,Q,N,M) CALL Form_A(C,D,E,P,S,Q,N,M) IF (Proc.EQ.Observe) THEN CALL TRI_TO_A(D,E,A,N) CALL Aoutput(A,K,N) WRITE(9,*)' ' WRITE(9,*)' Press the key. ' PAUSE ENDIF K = K+1 REPEAT REPEAT CALL Find_Shift(D,E,SH,R1,R2,M) D(N-1) = R1 + SHIFT D(N) = R2 + SHIFT DO J=1,(N-1) E(J) = 0 ENDDO CALL TRI_TO_A(D,E,A,N) RETURN END SUBROUTINE MATIN(Ach,A,N,InRC) PARAMETER(MaxR = 10) INTEGER Col,InRC,N,Row REAL A CHARACTER Ach*1 DIMENSION A(1:MaxR,1:MaxR) DO 20 Row=1,N DO 10 Col=1,N A(Row,Col)=0 10 CONTINUE 20 CONTINUE WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach WRITE(9,*)' ' IF (InRC.EQ.0) THEN DO 40 Row=1,N WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row WRITE(9,*)' ' DO 30 Col=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 30 CONTINUE 40 CONTINUE ELSE DO 60 Col=1,N WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col WRITE(9,*)' ' DO 50 Row=1,N WRITE(9,*) Ach,'(',Row,',',Col,') = ' READ(9,*) A(Row,Col) 50 CONTINUE 60 CONTINUE ENDIF RETURN END SUBROUTINE PrintResults(A,A1,D,N) PARAMETER(MaxR = 10) REAL A,A1,D DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) DIMENSION D(1:MaxR) INTEGER J,N,R,C WRITE(9,*)' ' WRITE(9,*)'The matrix A is:' WRITE(9,*)' ' DO 10 R=1,N WRITE(9,999) ( A1(R,C), C=1,N ) ENDDO WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'The similar diagonal matrix is:' DO 10 R=1,N WRITE(9,999) ( A(R,C), C=1,N ) ENDDO WRITE(9,*)' ' WRITE(9,*)' ' DO J=1,N WRITE(9,*)'the ',J,'th eigenvalue is ',D(J) ENDDO WRITE(9,*)' ' PAUSE RETURN 999 FORMAT(1X,5F15.7,4X) END SUBROUTINE MESSAGE(Meth) WRITE(9,*)' THE QL ALGORITHM' WRITE(9,*)' ' WRITE(9,*)' Assume that A is real, symmetric' WRITE(9,*)' ' WRITE(9,*)'and tridiagonal. Start with A = A.' WRITE(9,*)' 1 ' WRITE(9,*)' ' WRITE(9,*)'The QL method constructs a sequence ' WRITE(9,*)' ' WRITE(9,*)'of orthogonal matrices {Q } such that' WRITE(9,*)' k ' WRITE(9,*)' ' WRITE(9,*)'A = Q L where L is lower-triangular.' WRITE(9,*)' k k k k' WRITE(9,*)'Then A is generated by' WRITE(9,*)' k+1 ' WRITE(9,*)' ' WRITE(9,*)'A = Q A Q for k = 1,2, ... .' WRITE(9,*)' k+1 k k k' WRITE(9,*)' ' WRITE(9,*)'The success of the method depends on' WRITE(9,*)' ' WRITE(9,*)'the result that lim A = D, where' WRITE(9,*)' k->oo k ' WRITE(9,*)' ' WRITE(9,*)'D is a diagonal matrix with the same ' WRITE(9,*)'eigenvalues as the original matrix A.' WRITE(9,*)' ' PAUSE RETURN END SUBROUTINE INPUTS(A,A1,N,InRC) PARAMETER(MaxR = 10) REAL A,A1 DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR) INTEGER N,InRC,Row,Col WRITE(9,*)' We will now proceed with the QL' WRITE(9,*)' ' WRITE(9,*)'method of iteration to reduce the' WRITE(9,*)' ' WRITE(9,*)'tridiagonal matrix A to diagonal form.' WRITE(9,*)' ' WRITE(9,*)'A must be a tridiagonal matrix of dimension N by N.' WRITE(9,*)' ' WRITE(9,*)'{N must be an integer between 2 and 10}' WRITE(9,*)' ' WRITE(9,*)'ENTER N = ' READ(9,*) N IF (N.LT.2) N = 2; IF (N.GT.10) N = 10 InRC = 0 WRITE(9,*)'WANT TO INPUT COLUMNS ? ' READ(9,'(A)') Ans IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN InRC = 1 ENDIF WRITE(9,*)' ' WRITE(9,*)'ENTER THE COEFFICIENT MATRIX A:' CALL MATIN('A',A,N,InRC) DO 30 Row=1,N DO 20 Col=1,N A1(Row,Col)=A(Row,Col) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE PROCESSES(Proc) INTEGER I,Proc,Auto,Manual,Observe Auto = 0 Manual = 1 Observe =2 WRITE(9,*)'To what extent do you want to control the program?' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'<1> The computer does it all automatically.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'<2> Computer selects, but we observe each step.' WRITE(9,*)' ' WRITE(9,*)' ' WRITE(9,*)'SELECT your choice for input < 1 - 2 > ? ' READ(9,*) I IF ((I.LT.1).OR.(3.LT.I)) I = 3 IF (I.EQ.1) Proc = Auto IF (I.EQ.2) Proc = Observe IF (I.EQ.3) Proc = Manual RETURN END