C*********************************************************************** C * C QR METHOD * C * C*********************************************************************** C C C C TO OBTAIN THE EIGENVALUES OF A SYMMETRIC, TRIDIAGONAL N BY N C MATRIX: C C A(1) B(2) C B(2) A(2) B(3) C . . . C . . . C . . . C B(N-1) A(N) B(N) C B(N) A(N) C C INPUT: N; A(1),...,A(N) (DIAGONAL OF A); B(2),...,B(N) C (OFF-DIAGONAL OF A); MAXIMUM NUMBER OF ITERATIONS M; C TOLERANCE TOL C C OUTPUT: EIGENVALUES OF A OR RECOMMENDED SPLITTING OF A C OR A MESSAGE THAT THE MAXIMUM NUMBER OF C ITERATIONS WERE EXCEEDED. C DIMENSION A(9),B(9),C(9),D(9),Q(9),R(9),S(9),X(9),Y(9),Z(9) C INPUT N, A(1),...,A(N), B(2),...,B(N), MM, TOL 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 the QR Method.' WRITE(6,*) 'The tridiagonal symmetric array A will be input' WRITE(6,*) 'from a text file in the order:' WRITE(6,*) '(diagonal): A(1), A(2), ..., A(n)' WRITE(6,*) '(subdiagonal): B(2), B(3), ..., B(n)' 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 READ(INP,*) (A(I),I=1,N), (B(I),I=2,N) 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 13 WRITE(6,*) 'Input the tolerance.' WRITE(6,*) ' ' READ(5,*) TOL IF (TOL .GT. 0.0) THEN OK = .TRUE. ELSE WRITE(6,*) 'Tolerance must be positive.' ENDIF GOTO 112 13 OK = .FALSE. 14 IF (OK) GOTO 15 WRITE(6,*) 'Input maximum number of iterations.' WRITE(6,*) READ(5,*) MM IF (MM .GT. 0) THEN OK = .TRUE. ELSE WRITE(6,*) 'Number must be a positive integer.' ENDIF GOTO 14 ELSE WRITE(6,*) 'The program will end so the input file can ' WRITE(6,*) 'be created. ' OK = .FALSE. ENDIF 15 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,*) 'QR METHOD' ZERO = 1.0E-30 C STEP 1 C SET THE ACCUMULATED SHIFT TO ZERO SHIFT = 0 K=1 C STEP 2 100 IF (K.GT.MM) GOTO 200 WRITE(OUP,3) K WRITE(OUP,4) (A(I),I=1,N), (B(I),I=2,N) C TEST FOR CONVERGENCE / DEFLATION C TEST FOR POSSIBLE SPLIT C STEP 3 IF (ABS(B(N)).LE.TOL) THEN A(N)=A(N)+SHIFT WRITE(OUP,6) A(N) N=N-1 END IF IF(ABS(B(2)).LE.TOL) THEN A(1)=A(1)+SHIFT WRITE(OUP,6) A(1) N=N-1 A(1)=A(2) DO 10 I=2,N A(I)=A(I+1) 10 B(I)=B(I+1) END IF M=N-1 IF(M.GE.2) THEN DO 20 I=3,M IF(ABS(B(I)).LE.TOL) THEN WRITE(OUP,5) GOTO 400 END IF 20 CONTINUE END IF C STEP 4 C COMPUTE SHIFT B1=-(A(N-1)+A(N)) C1=A(N)*A(N-1)-B(N)*B(N) IF (ABS(C1).LE.ZERO) THEN WRITE(OUP,7) C GOTO 400 END IF D1=B1*B1-4*C1 IF(D1.LT.0.0) THEN WRITE(OUP,8) GOTO 400 END IF D1 = SQRT(D1) IF(B1.GT.0.0) THEN X1=-2*C1/(B1+D1) X2=-(B1+D1)/2 ELSE X1=(D1-B1)/2 X2=2*C1/(D1-B1) END IF C IF N IS 2 THEN WE HAVE COMPUTED THE 2 EIGENVALUES IF(N.EQ.2) THEN X1=X1+SHIFT X2=X2+SHIFT WRITE(OUP,6) X1 WRITE(OUP,6) X2 GOTO 400 END IF C STEP 5 C ACCUMULATE SHIFT IF(ABS(A(N)-X1).GE.ABS(A(N)-X2)) X1=X2 SHIFT=SHIFT+X1 WRITE(OUP,9) X1 C STEP 6 C PERFORM SHIFT DO 30 I=1,N 30 D(I)=A(I)-X1 C STEP 7 C COMPUTE R(K) X(1)=D(1) Y(1)=B(2) DO 40 J=2,N Z(J-1)=SQRT(X(J-1)*X(J-1)+B(J)*B(J)) C(J)=X(J-1)/Z(J-1) S(J)=B(J)/Z(J-1) Q(J-1)=C(J)*Y(J-1)+S(J)*D(J) X(J)=C(J)*D(J)-S(J)*Y(J-1) IF(J.NE.N) THEN R(J-1)=S(J)*B(J+1) Y(J)=C(J)*B(J+1) END IF 40 CONTINUE M1=N-2 C STEP 8 C COMPUTE NEW A Z(N)=X(N) A(1)=S(2)*Q(1)+C(2)*Z(1) B(2)=S(2)*Z(2) IF(N.GT.2) THEN DO 50 J=2,M A(J)=S(J+1)*Q(J)+C(J+1)*C(J)*Z(J) 50 B(J+1)=S(J+1)*Z(J+1) END IF A(N)=C(N)*X(N) K=K+1 GOTO 100 200 WRITE(OUP,12) C STEP 9 C THE PROCESS IF COMPLETE 400 CLOSE(UNIT=5) CLOSE(UNIT=OUP) IF(OUP.NE.6) CLOSE(UNIT=6) STOP 3 FORMAT(1X,I2) 4 FORMAT(1X,'Vectors A-B',7(2X,E15.8)) 5 FORMAT(1X,'TIME TO SPLIT') 6 FORMAT(1X,'EIGENVALUE = ',E15.8,' DEFLATE AND MOVE ON') 7 FORMAT(1X,'NEARLY SINGULAR') 8 FORMAT(1X,'COMPLEX ROOTS') 9 FORMAT(1X,'SHIFT = ',E15.8) 11 FORMAT(7(2X,E15.8)) 12 FORMAT(1X,'FAILURE') END