C*********************************************************************** C * C WIELANDT'S DEFLATION * C * C*********************************************************************** C C C C TO APPROXIMATE THE SECOND MOST DOMINANT EIGENVALUE AND AN C ASSOCIATED EIGENVECTOR OF THE N BY N MATRIX A GIVEN AN C APPROXIMATION LAMBDA TO THE DOMINANT EIGENVALUE, AN C APPROXIMATION V TO A CORRESPONDING EIGENVECTOR AND A VECTOR X C BELONGING TO R**(N-1): C C INPUT: DIMENSION N; MATRIX A; APPROXIMATE EIGENVALUE LAMBDA; C APPROXIMATE EIGENVECTOR V BELONGING TO R**N; VECTOR X C BELONGING TO R**(N-1). C C OUTPUT: APPROXIMATE EIGENVALUE MU; APPROXIMATE EIGENVECTOR U OR C A MESSAGE THAT THE METHOD FAILS. C C INITIALIZATION DIMENSION A(3,3),B(2,2),X(2),V(3),W(3),Y(2),VV(3) CHARACTER NAME*14,NAME1*14,AA*1 INTEGER INP,OUP,FLAG LOGICAL OK OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL') OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL') WRITE(6,*) 'This is Wielandt Deflation.' WRITE(6,*) 'The array will be input from a text file in the' WRITE(6,*) ' order: A(1,1), A(1,2), ..., A(1,N), A(2,1),' WRITE(6,*) ' A(2,2), ..., A(2,N)..., A(N,1), A(N,2),' WRITE(6,*) ' ..., A(N,N) ' WRITE(6,*) 'Next place the approximate eigenvector' WRITE(6,*) 'V(1) ... V(n), follow it with an' WRITE(6,*) 'approximate eigenvalue and finally with an' WRITE(6,*) 'initial approximate eigenvector of dimension' WRITE(6,*) 'n-1: X(1), ..., X(n-1)' WRITE(6,*) 'Place as many entries as desired on each line,' WRITE(6,*) ' but separate entries with at least one blank.' OK = .FALSE. WRITE(6,*) 'Has the input file been created?' WRITE(6,*) 'Enter Y or N - letter within quotes ' WRITE(6,*) ' ' READ(5,*) AA IF (( AA .EQ. 'Y' ) .OR.( AA .EQ. 'y' )) THEN WRITE(6,*) 'Input the file name in the form - ' WRITE(6,*) 'drive:name.ext contained in quotes' WRITE(6,*) 'as example: ''A:DATA.DTA'' ' WRITE(6,*) ' ' READ(5,*) NAME INP = 4 OPEN(UNIT=INP,FILE=NAME,ACCESS='SEQUENTIAL') OK = .FALSE. 19 IF (OK) GOTO 111 WRITE(6,*) 'Input the dimension n.' WRITE(6,*) READ(5,*) N IF (N .GT. 0) THEN M = N-1 READ(INP,*) ((A(I,J),J=1,N),I=1,N) READ(INP,*) (V(I),I=1,N) READ(INP,*) XMU READ(INP,*) (X(I),I=1,M) OK = .TRUE. CLOSE(UNIT=INP) ELSE WRITE(6,*) 'The number must be a positive integer' ENDIF GOTO 19 111 OK = .FALSE. 112 IF (OK) GOTO 113 WRITE(6,*) 'Input the tolerance for the power method.' WRITE(6,*) ' ' READ(5,*) TOL IF (TOL .GT. 0.0) THEN OK = .TRUE. ELSE WRITE(6,*) 'Tolerance must be positive.' ENDIF GOTO 112 113 OK = .FALSE. 114 IF (OK) GOTO 115 WRITE(6,*) 'Input maximum number of iterations.' WRITE(6,*) READ(5,*) NN IF (NN .GT. 0) THEN OK = .TRUE. ELSE WRITE(6,*) 'Number must be a positive integer.' ENDIF GOTO 114 ELSE WRITE(6,*) 'The program will end so the input file can ' WRITE(6,*) 'be created. ' OK = .FALSE. ENDIF 115 IF(.NOT. OK) GOTO 400 WRITE(6,*) 'Select output destinations: ' WRITE(6,*) '1. Screen ' WRITE(6,*) '2. Text file ' WRITE(6,*) 'Enter 1 or 2 ' WRITE(6,*) ' ' READ(5,*) FLAG IF ( FLAG .EQ. 2 ) THEN WRITE(6,*) 'Input the file name in the form - ' WRITE(6,*) 'drive:name.ext' WRITE(6,*) 'with the name contained within quotes' WRITE(6,*) 'as example: ''A:OUTPUT.DTA'' ' WRITE(6,*) ' ' READ(5,*) NAME1 OUP = 3 OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW') ELSE OUP = 6 ENDIF WRITE(OUP,*) 'WIELANDT DEFLATION' WRITE(OUP,5) WRITE(OUP,6) ((A(I,J),J=1,N),I=1,N) WRITE(OUP,7) WRITE(OUP,6) (V(I),I=1,N) WRITE(OUP,8) XMU WRITE(OUP,9) (X(I),I=1,M) C STEP 1 I=1 AMAX=ABS(V(1)) DO 100 J=2,N IF(ABS(V(J)).GT.AMAX) THEN I=J AMAX=ABS(V(J)) END IF 100 CONTINUE IF (AMAX.LT.1.0E-20) THEN WRITE(OUP,10) GOTO 400 END IF I1=I-1 N1=N-1 C STEP 2 IF (I .NE. 1) THEN DO 20 K=1,I1 DO 20 J=1,I1 20 B(K,J) = A(K,J)-V(K)*A(I,J)/V(I) END IF C STEP 3 IF(I.NE.1 .AND. I.NE.N) THEN DO 30 K=I,N1 DO 30 J=1,I1 B(K,J) = A(K+1,J)-V(K+1)*A(I,J)/V(I) 30 B(J,K) = A(J,K+1)-V(J)*A(I,K+1)/V(I) END IF C STEP 4 IF (I.NE.N) THEN DO 40 K=I,N1 DO 40 J=I,N1 40 B(K,J) = A(K+1,J+1)-V(K+1)*A(I,J+1)/V(I) END IF WRITE(OUP,11) WRITE(OUP,12) ((B(L1,L2),L2=1,M),L1=1,M) C STEP 5 CALL POWER(M,B,X,Y,YMU,NN,TOL) C Y IS USED IN PLACE OF W' C STEP 6 IF (I.NE.1) THEN DO 50 K=1,I1 50 W(K) = Y(K) END IF C STEP 7 W(I) = 0 C STEP 8 IF (I.NE.N) THEN I2 = I+1 DO 60 K=I2,N 60 W(K) = Y(K-1) END IF C STEP 9 DO 70 K=1,N S = 0 DO 80 J=1,N 80 S = S+A(I,J)*W(J) S = S/V(I) C COMPUTE EIGENVECTOR C VV IS USED IN PLACE OF U HERE 70 VV(K) = (YMU-XMU)*W(K)+S*V(K) WRITE(OUP,13) YMU WRITE(OUP,14) (Y(I),I=1,M) WRITE(OUP,15) (VV(I),I=1,N) 400 CLOSE(UNIT=5) CLOSE(UNIT=OUP) IF(OUP.NE.6) CLOSE(UNIT=6) STOP 5 FORMAT(1X,'ORIGINAL MATRIX A') 6 FORMAT(1X,3(3X,F10.5)) 7 FORMAT(1X,'APPROXIMATE EIGENVECTOR') 8 FORMAT(1X,'APPROX. DOMINANT EIGENVALUE',3X,E15.8) 9 FORMAT(1X,'INITIAL APPROX.',2(3X,F10.5)) 10 FORMAT(1X,'INPUT ERROR') 11 FORMAT(1X,'REDUCED MATRIX B') 12 FORMAT(1X,2(3X,E15.8)) 13 FORMAT(1X,'APPROX. TO SECOND EIGENVALUE',3X,E15.8) 14 FORMAT(1X,'EIGENVECTOR FOR B',2(3X,E15.8)) 15 FORMAT(1X,'EIGENVECTOR FOR A',3(3X,E15.8)) END C SUBROUTINE POWER(N,A,X,Y,YMU,NN,TOL) C SEE ALGORITHM 9.1 C NOTATION: N X N MATRIX A, ITERATION VECTOR Y, NORMALIZED C APPROX. EIGENVECTOR X, APPROX. EIGENVALUE XMU DIMENSION A(N,N),X(N),Y(N) K = 1 LP = 1 AMAX = ABS(X(1)) DO 10 I=2,N IF(ABS(X(I)).GT.AMAX) THEN AMAX = ABS(X(I)) LP = I END IF 10 CONTINUE DO 20 I=1,N 20 X(I)=X(I)/AMAX 110 IF (K.GT.NN) GOTO 200 DO 30 I=1,N Y(I) = 0.0 DO 30 J=1,N 30 Y(I) = Y(I)+A(I,J)*X(J) YMU=Y(LP) LP = 1 AMAX = ABS(Y(1)) DO 40 I=2,N IF(ABS(Y(I)).GT.AMAX) THEN AMAX = ABS(Y(I)) LP = I END IF 40 CONTINUE IF (AMAX.LE.TOL) THEN WRITE(6,100) STOP ENDIF ERR=0.0 DO 50 I=1,N T=Y(I)/Y(LP) IF(ABS(X(I)-T).GT.ERR) ERR=ABS(X(I)-T) 50 X(I) = T IF(ERR.LE.TOL) THEN DO 60 I=1,N 60 Y(I)=X(I) RETURN END IF K=K+1 GOTO 110 200 WRITE(6,101) STOP 100 FORMAT(1X,'ZERO VECTOR-ERROR') 101 FORMAT(1X,'DIVERGENCE-ERROR') END