C #QREGB V1A 17-APR-72 C LAST UPDATE C SUBROUTINE QREGB(A,Z,D,E,N,M,EPM) DIMENSION A(M,M),Z(M,M),D(M),E(M) C C *PURPOSE. C TO DIAGONALIZE A REAL SYMMETRIC MATRIX. C C *PARAMETERS. C C A - REAL ARRAY OF DIMENSION "M" X "M" WHICH CONTAINS THE C SYMMETRIC REAL "N" X "N" MATRIX TO BE DIAGONALIZED IN C ITS TOP LEFT CORNER. IF "A" AND "U" ARE DISTINCT "A" C IS LEFT UNCHANGED BY THE CALL. C U - REAL ARRAY OF DIMENSION "M" X "M" WHICH IS RETURNED WITH C THE "N" EIGENVECTORS IN THE FIRST "N" COLUMNS. THAT C IS THE K'TH EIGENVECTOR IS U(*,K). IF "A" IS NOT REQUIRED C AFTER THE CALL, "A" AND "U" MAY BE THE SAME ARRAY IN C THE CALL THUS SAVING SPACE. C D - REAL ARRAY OF DIMENSION "N" OR MORE WHICH IS RETURNED WITH C THE "N" ORDERED EIGENVALUES. THE EIGENVALUE D(K) C CORRESPONDS TO THE EIGENVECTOR U(*,K). C WRK - REAL ARRAY OF AT LEAST "N" REAL ELEMENTS USED AS WORKSPACE. C N - THE LOGICAL SIZE OF THE MATRIX TO BE DIAGONALIZED. C M - THE PHYSICAL SIZE OF THE ARRAYS CONTAINING THE MATRIX. C EPM - A TOLERANCE PARAMETER WHICH SHOULD BE SET BY THE USER TO C 100 TO 1000 TIMES THE MACHINE SIGNIFICANCE C C INPUT PARAMETERS: C A, N, M, EPM C OUTPUT PARAMETERS: C U, D, WRK C C *METHOD. C A HOUSEHOLDER TRANSFORMATION TO HESSENBERG FORM IS CARRIED OUT C AND FOLLOWED BY Q-R ITERATIONS. A MAXIMUM NUMBER OF ITERATIONS OF C 10*N IS SET AND IF THIS EXCEEDED EPM IS SET TO -1. THIS CASE HAS C NEVER BEEN FOUND TO OCCUR IN PRACTICE. C C *ACCURACY. C THIS DEPENDS OF EPM BUT 5 SIGNIFICANT DIGIT ACCURACY CAN BE C OBTAINED FOR N AROUND 20 AND EPM AROUND 1.E-4. C C *RESTRICTIONS. C C *ERROR CONDITIONS. C EPM IS RETURNED AS -1 IF THE Q-R PROCESS DOES NOT CONVERGE. THIS C HAS BEVER BEEN FOUND TO HAPPEN. NORMALLY EPM IS LEFT UNCHANGED. C C *NON STANDARD ROUTINES CALLED. C C *TYPICAL TIMES. C THE TIME IS ROUGHLY PROPORTIONAL TO N**3. C C *ORIGIN. C THIS ROUTINE HAS BEEN TAKEN WITH ONLY SLIGHT CHANGES C FORM "EIGQLT" BY J.P.M.BAILEY OF CHEM DIV DSIR. NEW ZEALAND. C THE CHANGES BY M.R.MANNING. C C *COMMENTS. C THIS ROUTINE IS GENERALLY FASTER AND MORE ACCURATE THAN THE JACOBI C METHODS AND QREGB SHOULD BE USED IN PREFERENCE TO JACOB EXCEPT C FOR SMALL MATRICES AND TIGHT CORE REQUIREMENTS. C C C #END. C C #START C C COPY MATRIX A INTO Z. C DO 20 I=1,N DO 20 J=1,N 20 Z(I,J)=A(I,J) I=N DO 103 IP=2,N L=I-2 F=Z(I,I-1) G=0.0 IF(L) 104,104,105 105 DO 106 K=1,L X=Z(I,K) 106 G=G+X*X 104 H=G+F*F IF(G) 207,207,107 207 E(I)=F D(I)=0.0 GO TO 108 107 L=L+1 G=SQRT(H) IF(F) 111,110,110 110 G=-G 111 E(I)=G H=1.0/(H-F*G) Z(I,I-1)=F-G F=0.0 C REMEMBER TRANSFORMATION PARAMETERS DO 115 J=1,L Z(J,I)=H*Z(I,J) G=0.0 DO 116 K=1,J 116 G=G+Z(J,K)*Z(I,K) IF(J.EQ.L) GO TO 117 J1=J+1 DO 118 K=J1,L 118 G=G+Z(K,J)*Z(I,K) 117 E(J)=G*H F=F+G*Z(J,I) 115 CONTINUE HH=0.5*F*H DO 119 J=1,L F=Z(I,J) G=E(J)-HH*F E(J)=G DO 119 K=1,J 119 Z(J,K)=Z(J,K)-F*E(K)-G*Z(I,K) D(I)=H 108 I=I-1 103 CONTINUE C C RECOVER ORTHOGONAL TRANSFORMATION MATRIX FOR C COMPLETE HOUSEHOLDER TRANSFORMATION. C D(1)=0.0 E(1)=0.0 DO 125 I=1,N L=I-1 IF(D(I)) 124,126,124 124 DO 127 J=1,L G=0.0 DO 128 K=1,L 128 G=G+Z(I,K)*Z(K,J) DO 129 K=1,L 129 Z(K,J)=Z(K,J)-G*Z(K,I) 127 CONTINUE 126 D(I)=Z(I,I) Z(I,I)=1.0 IF(L) 125,125,130 130 DO 131 J=1,L Z(J,I)=0.0 131 Z(I,J)=0.0 125 CONTINUE DO 132 I=2,N 132 E(I-1)=E(I) C C D NOW CONTAINS DIAGONAL MATRIX ELEMENTS, AND E SUBDIAGONAL C Z CONTAINS ORTHOGONAL TRANSFORMATION TO BECOME EIGENVECTORS. C C C START Q-R ITERATIONS. C MAXST=10*N F=0.0 B=0.0 E(N)=0.0 DO 140 L=1,N EL=E(L) DL=D(L) J=0 H=EPM*(ABS(DL)+ABS(EL)) IF(B.LT.H) B=H DO 141 M=L,N IF(ABS(E(M)).LE.B) GO TO 142 141 CONTINUE 142 IF(M.EQ.L) GO TO 143 C JUST TO STOP INFINITE LOOPS...SHOULD NEVER BE NEEDED. 144 IF(J.GT.MAXST) GO TO 145 J=J+1 P=0.5*(D(L+1)-DL)/EL R=SQRT(P*P+1.0) IF(P) 146,147,147 146 H=P-R GO TO 148 147 H=P+R 148 H=DL-EL/H DO 149 I=L,N 149 D(I)=D(I)-H DL=DL-H F=F+H P=D(M) C=1.0 S=0.0 MM1=M-1 I=MM1 DO 150 IP=L,MM1 I1=I+1 EI=E(I) G=C*EI H=C*P IF(ABS(P).LT.ABS(EI)) GO TO 151 C=EI/P R=SQRT(C*C+1.0) E(I1)=S*P*R S=C/R C=1.0/R GO TO 152 151 C=P/EI R=SQRT(C*C+1.0) E(I1)=S*R*EI S=1.0/R C=C*S 152 P=C*D(I)-S*G D(I1)=H+S*(C*G+S*D(I)) DO 155 K=1,N H=Z(K,I1) Z(K,I1)=S*Z(K,I)+C*H 155 Z(K,I)=C*Z(K,I)-S*H 150 I=I-1 EL=S*P E(L)=EL DL=C*P D(L)=DL IF(ABS(EL).GT.B) GO TO 144 143 D(L)=D(L)+F 140 CONTINUE C C ORDER EIGENVALUES ALGEBRAICALLY. C MM1=N-1 DO 160 I=1,MM1 K=I P=D(I) I1=I+1 DO 165 J=I1,N IF(D(J).GE.P) GO TO 165 K=J P=D(J) 165 CONTINUE IF(K.EQ.I) GO TO 160 D(K)=D(I) D(I)=P DO 168 J=1,N P=Z(I,J) Z(I,J)=Z(K,J) 168 Z(K,J)=P 160 CONTINUE C C FINISHED...NORMAL RETURN C RETURN C C ALARM RETURN. C 145 EPM= -1. RETURN END