SUBROUTINE CDIAG(A,VALUE,VEC,N, NEED) C C TO FIND THE EIGENVALUES AND EIGENVECTORS OF A HERMITIAN MATRIX. REAL VALUE(*),H COMPLEX W(6000) COMPLEX A(N,*),VEC(N,*) COMPLEX FM06AS IA=N IV=N C C REDUCE MATRIX TO A TRI-DIAGONAL HERMITIAN MATRIX. CALL ME08A(A,W,W(N+1),N,IA,W(2*N+1)) C C FIND THE EIGENVALUES AND EIGENVECTORS OF THE TRI-DIAGONAL MATRIX CALL EC08C(W,W(N+1),VALUE,VEC,N,IV,W(2*N+1)) IF(NEED.EQ.0) GOTO 50 IF(N.LT.2)RETURN C C THE EIGENVECTORS OF THE ORIGINAL MATRIX ARE NOW FOUND BY C BACK TRANSFORMATION USING INFORMATION STORE IN THE UPPER C TRIANGLE OF MATRIX A (BY ME08) DO 40 II=3,N I=N-II+1 H=W(N+I+1)*CONJG(A(I,I+1)) IF(H)10,40,10 10 DO30L=1,N I1=I+1 S=FM06AS(N-I,A(I,I+1),IA,VEC(I+1,L),1) S=S/H DO20K=I1,N 20 VEC(K,L)=VEC(K,L)+CONJG(A(I,K))*S 30 CONTINUE 40 CONTINUE 50 CALL SORT(VALUE,VEC,N) RETURN END SUBROUTINE EA08C(A,B,VALUE,VEC,M,IV,W) C STANDARD FORTRAN 66 (A VERIFIED PFORT SUBROUTINE) DIMENSION A(*),B(*),VALUE(*),VEC(*),W(*) DATA EPS/1.E-6/,A34/0.0/ C THIS USES QR ITERATION TO FIND THE EIGENVALUES AND EIGENVECTORS C OF THE SYMMETRIC TRIDIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS ARE C A(I),I=1,M AND OFF-DIAGONAL ELEMENTS ARE B(I),I=2,M. THE ARRAY C W IS USED FOR WORKSPACE AND MUST HAVE DIMENSION AT LEAST 2*M. C WE TREAT VEC AS IF IT HAD DIMENSIONS (IV,M). SML=EPS*FLOAT(M) CALL EA09C(A,B,W(M+1),M,W) C SET VEC TO THE IDENTITY MATRIX. DO 20 I=1,M VALUE(I)=A(I) W(I)=B(I) K=(I-1)*IV+1 L=K+M-1 DO 10 J=K,L 10 VEC(J)=0. KI=K+I-1 20 VEC(KI)=1. K=0 IF(M.EQ.1)RETURN DO 100 N3=2,M N2=M+2-N3 C EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2 MN2=M+N2 ROOT=W(MN2) DO 80 ITER=1,20 BB=(VALUE(N2)-VALUE(N2-1))*0.5 CC=W(N2)*W(N2) A22=VALUE(N2) IF(CC.NE.0.0)A22=A22+CC/(BB+SIGN(1.0,BB)*SQRT(BB*BB+CC)) DO 30 I=1,N2 MI=M+I IF(ABS(ROOT-A22).LE.ABS(W(MI)-A22))GO TO 30 ROOT=W(MI) MN=M+N2 W(MI)=W(MN) W(MN)=ROOT 30 CONTINUE DO 40 II=2,N2 N1=2+N2-II IF(ABS(W(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*SML)GO 1 TO 50 40 CONTINUE N1=1 50 IF(N2.EQ.N1) GO TO 100 N2M1=N2-1 IF(ITER.GE.3)ROOT=A22 K=K+1 A22=VALUE(N1) A12=A22-ROOT A23=W(N1+1) A13=A23 DO70 I=N1,N2M1 A33=VALUE(I+1) IF(I.NE.N2M1)A34=W (I+2) S=SIGN(SQRT(A12*A12+A13*A13),A12) SI=A13/S CO=A12/S JK=I*IV+1 J1=JK-IV J2=J1+MIN0(M,I+K)-1 DO 60 JI=J1,J2 V1=VEC(JI) V2=VEC(JK) VEC(JI)=V1*CO+V2*SI VEC(JK)=V2*CO-V1*SI 60 JK=JK+1 IF(I.NE.N1) W(I)=S A11=CO*A22+SI*A23 A12=CO*A23+SI*A33 A13=SI*A34 A21=CO*A23-SI*A22 A22=CO*A33-SI*A23 A23=CO*A34 VALUE(I)=A11*CO+A12*SI A12=-A11*SI+A12*CO W(I+1)=A12 70 A22=A22*CO-A21*SI 80 VALUE(N2)=A22 WRITE(6,90) 90 FORMAT(48H1CYCLE DETECTED IN SUBROUTINE EA08 -STOPPING NOW) STOP 100 CONTINUE C RAYLEIGH QUOTIENT DO 120 J=1,M K=(J-1)*IV XX=VEC(K+1)**2 XAX=XX*A(1) DO 110 I=2,M KI=K+I XX=XX+VEC(KI)**2 110 XAX=XAX+VEC(KI)*(2.*B(I)*VEC(KI-1)+A(I)*VEC(KI)) 120 VALUE(J)=XAX/XX RETURN END SUBROUTINE EA09C(A,B,VALUE,M,OFF) C STANDARD FORTRAN 66 (A VERFIED PFORT SUBROUTINE) DIMENSION A(*),B(*),VALUE(*),OFF(*) DATA A34/0.0/,EPS/1.0E-6/ SML=EPS*FLOAT(M) VALUE(1)=A(1) IF(M.EQ.1)RETURN DO 10 I=2,M VALUE(I)=A(I) 10 OFF(I)=B(I) C EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2 MAXIT=10*M DO 80 ITER=1,MAXIT DO 40 N3=2,M N2=M+2-N3 DO 20 II=2,N2 N1=2+N2-II IF(ABS(OFF(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*SML) 1GO TO 30 20 CONTINUE N1=1 30 IF(N2.NE.N1) GO TO 50 40 CONTINUE RETURN C ROOT IS THE EIGENVALUE OF THE BOTTOM 2*2 MATRIX THAT IS NEAREST C TO THE LAST MATRIX ELEMENT AND IS USED TO ACCELERATE THE C CONVERGENCE 50 BB=(VALUE(N2)-VALUE(N2-1))*0.5 CC=OFF(N2)*OFF(N2) SBB=1. IF(BB.LT.0.)SBB=-1. ROOT=VALUE(N2)+CC/(BB+SBB*SQRT(BB*BB+CC)) N2M1=N2-1 60 A22=VALUE(N1) A12=A22-ROOT A23=OFF(N1+1) A13=A23 DO 70 I=N1,N2M1 A33=VALUE(I+1) IF(I.NE.N2M1)A34=OFF(I+2) S=SQRT(A12*A12+A13*A13) SI=A13/S CO=A12/S IF(I.NE.N1)OFF(I)=S A11=CO*A22+SI*A23 A12=CO*A23+SI*A33 A13=SI*A34 A21=CO*A23-SI*A22 A22=CO*A33-SI*A23 A23=CO*A34 VALUE(I)=A11*CO+A12*SI A12=-A11*SI+A12*CO OFF(I+1)=A12 70 A22=A22*CO-A21*SI 80 VALUE(N2)=A22 WRITE(6,90) 90 FORMAT(39H1LOOPING DETECTED IN EA09-STOPPING NOW ) STOP END SUBROUTINE EC08C(A,B,VALUE,VEC,N,IV,W) C C TO FIND THE EIGENVALUES AND VECTORS OF A TRI-DIAGONAL C HERMITIAN MATRIX. REAL VALUE(*),W(*),PCK(2),ONE,ZERO,VEC(*) COMPLEX A(*),B(*),DN,UPCK EQUIVALENCE (PCK(1),UPCK) C WE TREAT VEC AS IF IT IS DEFINED AS COMPLEX VEC(IV,N) C IN THE CALLING PROGRAM. DATA ONE, ZERO/1.0,0.0/ IV2=IV+IV N2=N+N IL=IV2*(N-1)+1 W(1)=A(1) C C THE HERMITIAN FORM IS TRANSFORMED INTO A REAL FORM IF(N-1)80,80,10 10 DO 20 I=2,N W(I)=A(I) 20 W(I+N)=CABS(B(I)) C C FIND THE EIGENVALUES AND VECTORS OF THE REAL FORM 30 CALL EA08C(W,W(N+1),VALUE,VEC,N,IV2,W(N2+1)) C C THE VECTORS IN VEC AT THIS POINT ARE REAL,WE NOW EXPAND THEM C INTO VEC AS THOUGH THEY WERE COMPLEX. DO 50 I=1,IL,IV2 DO 40 J=1,N K=N-J L=K+K 40 VEC(I+L)=VEC(I+K) 50 VEC(I+1)=ZERO IF(N.LE.1)GO TO 80 DN=ONE K=1 C C TRANSFORM VECTORS OF REAL FORM TO THOSE OF COMPLEX FORM. DO 70 I=4,N2,2 K=K+1 UPCK=ONE IF(W(K+N).NE.ZERO)UPCK=DN*CONJG(B(K))/W(K+N) I1=IL-1+I DO 60 J=I,I1,IV2 VEC(J)=VEC(J-1)*PCK(2) 60 VEC(J-1)=VEC(J-1)*PCK(1) 70 DN=UPCK 80 RETURN END COMPLEX FUNCTION FM06AS(N,A,IA,B,IB) IMPLICIT COMPLEX (A-H,O-Z) COMPLEX A(*), B(*) ******************************************************* * * FM06AS - A FUNCTION ROUTINE TO COMPUTE THE VALUE OF THE * INNER PRODUCT, OR DOT PRODUCT, OF TWO SINGLE PRECISION * COMPLEX VECTORS, ACCUMULATING THE INTERMEDIATE PRODUCTS * DOUBLE PRECISION. THE ELEMENTS OF EACH VECTOR CAN BE * STORED IN ANY FIXED DISPLACEMENT FROM NEIGHBOURING * ELEMENTS. * * COMPUTES: SUM(J=1,N) A((J-1)*IA+1)*B((J-1)*IB+1) * * W = FM06AS(N,A,IA,B,IB) * * N INTEGER SCALAR; (USER:*); LENGTH OF THE VECTORS A AND B. * IF N <= 0 THE INNER PRODUCT VALUE IS DEFINED TO BE ZERO. * A COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY * CONTAINING THE 1ST VECTOR. THE FORTRAN CONVENTION OF STORING * REAL AND IMAGINARY PARTS IN ADJACENT WORDS IS ASSUMED. * IA INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF * AN ELEMENT IN THE ARRAY A TO ITS NEIGHBOUR, I.E. THE VECTOR * ELEMENTS ARE IN A(1), A(IA+1), A(2*IA+1),... * IF IA < 0 THE ELEMENTS ARE ASSUMED TO BE STORED IN * A(1-(N-1)*IA), A(1-(N-2)*IA),..., A(1-IA), A(1). * B COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY * CONTAINING THE SECOND VECTOR. TREAT LIKE A. * IB INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF * AN ELEMENT IN B TO ITS NEIGHBOUR. TREAT LIKE IA. * FM06AS COMPLEX FUNCTION; (*:FM06AS); THE INNER PRODUCT VALUE. * IT IS RETURNED DOUBLE PRECISION, THE REAL PART IN FLT PNT * REG 0 AND THE IMAGINARY PART IN FLT PNT REG 2. * * THIS ROUTINE IS WRITTEN IN FORTRAN. * *--------------------------------------------------------* SUM=(0.0,0.0) DO 10 I=1,N 10 SUM=SUM+A((I-1)*IA+1)*B((I-1)*IB+1) FM06AS=SUM RETURN END COMPLEX FUNCTION FM06BS(N,A,IA,B,IB) IMPLICIT COMPLEX (A-H,O-Z) COMPLEX A(*), B(*) ******************************************************* * * FM06BS - A FUNCTION ROUTINE TO COMPUTE THE VALUE OF THE * INNER PRODUCT, OR DOT PRODUCT, OF A SIGLE PRECISION * COMPLEX VECTORS, ACCUMULATING THE INTERMEDIATE PRODUCTS * DOUBLE PRECISION. THE ELEMENTS OF EACH VECTOR CAN BE * STORED IN ANY FIXED DISPLACEMENT FROM NEIGHBOURING * ELEMENTS. * * COMPUTES: SUM(J=1,N) A((J-1)*IA+1)*B((J-1)*IB+1) * * W = FM06BS(N,A,IA,B,IB) * * N INTEGER SCALAR; (USER:*); LENGTH OF THE VECTORS A AND B. * IF N <= 0 THE INNER PRODUCT VALUE IS DEFINED TO BE ZERO. * A COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY * CONTAINING THE 1ST VECTOR. THE FORTRAN CONVENTION OF STORING * REAL AND IMAGINARY PARTS IN ADJACENT WORDS IS ASSUMED. * IA INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF * AN ELEMENT IN THE ARRAY A TO ITS NEIGHBOUR, I.E. THE VECTOR * ELEMENTS ARE IN A(1), A(IA+1), A(2*IA+1),... * IF IA < 0 THE ELEMENTS ARE ASSUMED TO BE STORED IN * A(1-(N-1)*IA), A(1-(N-2)*IA),..., A(1-IA), A(1). * B COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY * CONTAINING THE SECOND VECTOR. TREAT LIKE A. * IB INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF * AN ELEMENT IN B TO ITS NEIGHBOUR. TREAT LIKE IA. * FM06BS COMPLEX FUNCTION; (*:FM06BS); THE INNER PRODUCT VALUE. * IT IS RETURNED DOUBLE PRECISION, THE REAL PART IN FLT PNT * REG 0 AND THE IMAGINARY PART IN FLT PNT REG 2. * * THIS ROUTINE IS WRITTEN IN FORTRAN. * *--------------------------------------------------------* SUM=(0.0,0.0) DO 10 I=1,N 10 SUM=SUM+A((I-1)*IA+1)*CONJG(B((I-1)*IB+1)) FM06BS=SUM RETURN END SUBROUTINE ME08A(A,ALPHA,BETA,N,IA,Q) COMPLEX A(IA,*),ALPHA(*),BETA(*),Q(*),CW,QJ,CZERO COMPLEX FM06AS,FM06BS REAL PP,ZERO,P5,H REAL S1,PP1 DATA ZERO/0.0/, P5 /0.50/ CZERO=ZERO C************************************************************** C THE REDUCTION OF FULL HERMITIAN MATRIX INTO TRI-DIAGONAL HERMITIAN C FORM IS DONE IN N-2 STEPS.AT THE I TH STEP ZEROS ARE INTRODUCED IN C THE I TH ROW AND COLUMNS WITHOUT DESTROYING PREVIOUSLY PRODUCED ZEROS C C H C AT THE I TH STEP WE HAVE A =P A P WITH P =I-U U /K C I I I-1 I I I I I C C WHERE U =(0,0,...,A (1+S /T ),A ,...,A ) C I I,I+1 I I I,I+2 I,N C 2 N 2 2 2 C S = SUM ] A ] K =S +T AND T = SQRT(]A ] S ) C I J=I+1 I,J I I I I I,I+1 I C C COMPUTATIONAL DETAILS AT THE I TH STAGE ARE (1) FORM S ,K THEN C I I C C H H H C (2) Q =A U /K (3) Q =Q -.5U (U Q /K ) (4) A =A -U Q -Q U C I I-1 I I I I I I I I I I-1 I I I I C C THE VECTORS U BEING APPROPIATELY IN A. IF(N.LE.0)GO TO 90 DO 10 J=1,N ALPHA(J)=A(J,J) DO 10 I=1,J 10 A(I,J)=CONJG(A(J,I)) IF(N-2)90,80,20 20 N2=N-2 DO 60 I=1,N2 I1=I+1 C (1) CW=FM06BS(N-I,A(I,I+1),IA,A(I,I+1),IA) PP=CW PP1=SQRT(PP) BETA(I+1)=CMPLX(-PP1,ZERO) S1=CABS(A(I,I+1)) IF(S1.GT.ZERO)BETA(I+1)=BETA(I+1)*A(I,I+1)/S1 IF(PP.LE.1.D-15)GO TO 60 H=PP+PP1*S1 A(I,I+1)=A(I,I+1)-BETA(I+1) C (2) DO 30 K=I1,N QJ=FM06AS(-(I-K),A(I+1,K),1,A(I,I+1),IA) QJ=FM06BS(N-K,A(K,K+1),IA,A(I,K+1),IA)+CONJG(QJ) 30 Q(K)=QJ/H C (3) CW=FM06AS(N-I,A(I,I+1),IA,Q(I+1),1) PP=CW*P5/H DO 40 K=I1,N 40 Q(K)=Q(K)-PP*CONJG(A(I,K)) C (4) DO 50 K=I1,N 50 CALL ME08B (A(K,K),Q(K),A(I,K),N-K+1,IA*2) 60 CONTINUE DO 70 I=2,N QJ=ALPHA(I) ALPHA(I)=A(I,I) 70 A(I,I)=QJ 80 BETA(N)=A(N-1,N) 90 RETURN END SUBROUTINE ME08B (A,Q,B,N,IA) REAL A(IA,*),Q(2,*),B(IA,*) DO 10 I=1,N A(1,I)=A(1,I) -Q(1,1)*B(1,I)+Q(2,1)*B(2,I) 1 -Q(1,I)*B(1,1)+Q(2,I)*B(2,1) 10 A(2,I)=A(2,I)-Q(2,1)*B(1,I)-Q(1,1)*B(2,I) 1 +Q(2,I)*B(1,1)+Q(1,I)*B(2,1) RETURN END SUBROUTINE SORT(VAL,VEC,N) COMPLEX VEC(N,*), SUM REAL VAL(*) DO 30 I=1,N X=1.E9 DO 10 J=I,N IF(VAL(J).LT.X) THEN K=J X=VAL(J) ENDIF 10 CONTINUE DO 20 J=1,N SUM=VEC(J,K) VEC(J,K)=VEC(J,I) 20 VEC(J,I)=SUM VAL(K)=VAL(I) VAL(I)=X 30 CONTINUE RETURN END