SUBROUTINE CMXDHE(AR,AI,D,ZR,ZI,N,NM,IERR) C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/cmxdhe.for C C=====Solve Eigenvalue problem for Complex Hermetian matrix(AR,AI) C (ZR,ZI) are returned as eigenvectors. The real eigenvalues are C returned as double precision vector D. N is the order of the C eigensystem; NM is the actual dimension of the matrices (50 C here). C DOUBLE PRECISION AR(50,N),AI(50,N),ZR(50,N),ZI(50,N) DOUBLE PRECISION D(50),E(50),E2(50),TAU(2,50) CALL HTRIDI(NM,N,AR,AI,D,E,E2,TAU) DO 1 I=1,N DO 1 J=1,N 1 ZR(I,J)=0.D0 DO 2 I=1,N 2 ZR(I,I)=1.D0 CALL IMTQL2(NM,N,D,E,ZR,IERR) CALL HTRIBK(NM,N,AR,AI,TAU,N,ZR,ZI) DO 3 I=1,N AR(I,I)=D(I) 3 AI(I,I)=0.D0 RETURN END C C ------------------------------------------------------------------ C SUBROUTINE HTRIDI(NM,N,AR,AI,D,E,E2,TAU) C INTEGER I,J,K,L,N,II,NM,JP1 DOUBLE PRECISION AR(NM,N),AI(NM,N),D(N),E(N),E2(N),TAU(2,N) DOUBLE PRECISION F,FI,G,GI,H,HH,SI,SCALE DOUBLE PRECISION DSQRT,CDABS,DABS COMPLEX*16 DCMPLX C C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF C THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968) C BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX C TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING C UNITARY SIMILARITY TRANSFORMATIONS. C C ON INPUT: C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT; C C N IS THE ORDER OF THE MATRIX; C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX. C ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. C C ON OUTPUT: C C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- C FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER C TRIANGLES. THEIR STRICT UPPER TRIANGLES AND THE C DIAGONAL OF AR ARE UNALTERED; C C D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX; C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO; C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED; C C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. C C ARITHMETIC IS REAL EXCEPT FOR THE USE OF THE SUBROUTINES C CDABS AND DCMPLX IN COMPUTING COMPLEX ABSOLUTE VALUES. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C TAU(1,N) = 1.0D0 TAU(2,N) = 0.0D0 C DO 100 I = 1, N 100 D(I) = AR(I,I) C C :::::::::: FOR I=N STEP -1 UNTIL 1 DO -- :::::::::: DO 300 II = 1, N I = N + 1 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 1) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(AR(I,K)) + DABS(AI(I,K)) C IF (SCALE .NE. 0.0D0) GO TO 140 TAU(1,L) = 1.0D0 TAU(2,L) = 0.0D0 130 E(I) = 0.0D0 E2(I) = 0.0D0 GO TO 290 C 140 DO 150 K = 1, L AR(I,K) = AR(I,K) / SCALE AI(I,K) = AI(I,K) / SCALE H = H + AR(I,K) * AR(I,K) + AI(I,K) * AI(I,K) 150 CONTINUE C E2(I) = SCALE * SCALE * H G = DSQRT(H) E(I) = SCALE * G F = CDABS(DCMPLX(AR(I,L),AI(I,L))) C :::::::::: FORM NEXT DIAGONAL ELEMENT OF MATRIX T :::::::::: IF (F .EQ. 0.0D0) GO TO 160 TAU(1,L) = (AI(I,L) * TAU(2,I) - AR(I,L) * TAU(1,I)) / F SI = (AR(I,L) * TAU(2,I) + AI(I,L) * TAU(1,I)) / F H = H + F * G G = 1.0D0 + G / F AR(I,L) = G * AR(I,L) AI(I,L) = G * AI(I,L) IF (L .EQ. 1) GO TO 270 GO TO 170 160 TAU(1,L) = -TAU(1,I) SI = TAU(2,I) AR(I,L) = G 170 F = 0.0D0 C DO 240 J = 1, L G = 0.0D0 GI = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J G = G + AR(J,K) * AR(I,K) + AI(J,K) * AI(I,K) GI = GI - AR(J,K) * AI(I,K) + AI(J,K) * AR(I,K) 180 CONTINUE C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L G = G + AR(K,J) * AR(I,K) - AI(K,J) * AI(I,K) GI = GI - AR(K,J) * AI(I,K) - AI(K,J) * AR(I,K) 200 CONTINUE C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H TAU(2,J) = GI / H F = F + E(J) * AR(I,J) - TAU(2,J) * AI(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = AR(I,J) G = E(J) - HH * F E(J) = G FI = -AI(I,J) GI = TAU(2,J) - HH * FI TAU(2,J) = -GI C DO 260 K = 1, J AR(J,K) = AR(J,K) - F * E(K) - G * AR(I,K) X + FI * TAU(2,K) + GI * AI(I,K) AI(J,K) = AI(J,K) - F * TAU(2,K) - G * AI(I,K) X - FI * E(K) - GI * AR(I,K) 260 CONTINUE C 270 DO 280 K = 1, L AR(I,K) = SCALE * AR(I,K) AI(I,K) = SCALE * AI(I,K) 280 CONTINUE C TAU(2,L) = -SI 290 HH = D(I) D(I) = AR(I,I) AR(I,I) = HH AI(I,I) = SCALE * SCALE * H 300 CONTINUE C RETURN C :::::::::: LAST CARD OF HTRIDI :::::::::: END C C ------------------------------------------------------------------ C SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,NM,MML,IERR DOUBLE PRECISION D(N),E(N),Z(NM,N) DOUBLE PRECISION B,C,F,G,P,R,S,MACHEP DOUBLE PRECISION DSQRT,DABS,DSIGN C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2, C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS C FULL MATRIX TO TRIDIAGONAL FORM. C C ON INPUT: C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT; C C N IS THE ORDER OF THE MATRIX; C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX; C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY; C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN C THE IDENTITY MATRIX. C C ON OUTPUT: C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT C UNORDERED FOR INDICES 1,2,...,IERR-1; C C E HAS BEEN DESTROYED; C C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED C EIGENVALUES; C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C :::::::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C MACHEP = 16.0D0**(-13) FOR LONG FORM ARITHMETIC C ON S360 :::::::::: MACHEP=D1MACH(IDUM) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C E(N) = 0.0D0 C DO 240 L = 1, N J = 0 C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: 105 DO 110 M = L, N IF (M .EQ. N) GO TO 120 IF (DABS(E(M)) .LE. MACHEP * (DABS(D(M)) + DABS(D(M+1)))) X GO TO 120 110 CONTINUE C 120 P = D(L) IF (M .EQ. L) GO TO 240 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: G = (D(L+1) - P) / (2.0D0 * E(L)) R = DSQRT(G*G+1.0D0) G = D(M) - P + E(L) / (G + DSIGN(R,G)) S = 1.0D0 C = 1.0D0 P = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II F = S * E(I) B = C * E(I) IF (DABS(F) .LT. DABS(G)) GO TO 150 C = G / F R = DSQRT(C*C+1.0D0) E(I+1) = F * R S = 1.0D0 / R C = C * S GO TO 160 150 S = F / G R = DSQRT(S*S+1.0D0) E(I+1) = G * R C = 1.0D0 / R S = S * C 160 G = D(I+1) - P R = (D(I) - G) * S + 2.0D0 * C * B P = S * R D(I+1) = G + P G = C * R - B C :::::::::: FORM VECTOR :::::::::: DO 180 K = 1, N F = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * F Z(K,I) = C * Z(K,I) - S * F 180 CONTINUE C 200 CONTINUE C D(L) = D(L) - P E(L) = G E(M) = 0.0D0 GO TO 105 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .GE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN C :::::::::: LAST CARD OF IMTQL2 :::::::::: END C C ------------------------------------------------------------------ C SUBROUTINE HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI) C INTEGER I,J,K,L,M,N,NM DOUBLE PRECISION AR(NM,N),AI(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M) DOUBLE PRECISION H,S,SI C C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF C THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968) C BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY HTRIDI. C C ON INPUT: C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT; C C N IS THE ORDER OF THE MATRIX; C C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- C FORMATIONS USED IN THE REDUCTION BY HTRIDI IN THEIR C FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR. C C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS; C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED; C C ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT: C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS C IN THEIR FIRST M COLUMNS. C C NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR C IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C :::::::::: TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC C TRIDIAGONAL MATRIX TO THOSE OF THE HERMITIAN C TRIDIAGONAL MATRIX. :::::::::: DO 50 K = 1, N C DO 50 J = 1, M ZI(K,J) = - ZR(K,J) * TAU(2,K) ZR(K,J) = ZR(K,J) * TAU(1,K) 50 CONTINUE C IF (N .EQ. 1) GO TO 200 C :::::::::: RECOVER AND APPLY THE HOUSEHOLDER MATRICES :::::::::: DO 140 I = 2, N L = I - 1 H = AI(I,I) IF (H .EQ. 0.0D0) GO TO 140 C DO 130 J = 1, M S = 0.0D0 SI = 0.0D0 C DO 110 K = 1, L S = S + AR(I,K) * ZR(K,J) - AI(I,K) * ZI(K,J) SI = SI + AR(I,K) * ZI(K,J) + AI(I,K) * ZR(K,J) 110 CONTINUE C S = S / H SI = SI / H C DO 120 K = 1, L ZR(K,J) = ZR(K,J) - S * AR(I,K) - SI * AI(I,K) ZI(K,J) = ZI(K,J) - SI * AR(I,K) + S * AI(I,K) 120 CONTINUE C 130 CONTINUE C 140 CONTINUE C 200 RETURN C :::::::::: LAST CARD OF HTRIBK :::::::::: END DOUBLE PRECISION FUNCTION D1MACH (IDUM) INTEGER IDUM C===This routine computes the unit round-off of the machine in Double C Precision. This is defined as the smallest positive machine real C real number, EPS, such that (1+EPS .GT. 1) & (1-EPS .LT. 1). C This particular computation of EPS is the work of Alan C. Hindmarsh. C DOUBLE PRECISION EPS, COMP EPS = 1.0D0 10 EPS = EPS*0.5D0 C---Note from Richard Hanson: C Make sure the following statement isn't optimized out of the loop. C could do this by making it a (separate) function and call here. COMP = 1.0D0 + EPS IF (COMP .NE. 1.0D0) GO TO 10 D1MACH = EPS*2.0D0 RETURN END