SUBROUTINE EIGCAB(A,B,N,N1,NE,NV,EPS,W,E,V,IER, lw,wd) parameter(MORB=200) integer n, n1, ne, nv, ier - ,lw(n1) - ,i,j,k,l double precision eps,e(ne) - ,wd(n1,7) complex*16 a(n1,n),b(n1,n),v(n1,nv),w(n1,7) - ,d,t,t1(MORB,MORB),t2(MORB,MORB) - ,sum,sum01,sum02,sum03,sum1,sum2 20 DO 24 K=N,1,-1 SUM = 0 DO 21 I=K+1,N SUM = SUM + W(I,1) * B(I,K)*DCONJG(B(I,K)) 21 CONTINUE D = B(K,K) - SUM if(real(d).lt.0) then write(*,*) '"D" became imaginary number !!' ier=5 end if W(K,1) = D D = 1 / D B(K,K) = cdSQRT(D) DO 23 J=1,K-1 SUM = 0 DO 22 I=K+1,N SUM = SUM + W(I,1) * DCONJG(B(I,K)) * B(I,J) 22 CONTINUE B(K,J) = ( B(K,J) - SUM ) * D 23 CONTINUE 24 CONTINUE C *************************************************************** C TRANSFORM THE MATRIX C (A NEW) = (L TRANSPOSED)**(-1) (A) (L)**(-1) C *************************************************************** do 40 k=n-1,1,-1 do 31 i=1,n do 30 j=1,i a(j,i)=dconjg(a(i,j)) 30 continue 31 continue do 33 j=1,k-1 sum1=0 do 32 i=k+1,n sum1=sum1+dconjg(b(i,k))*a(i,j) 32 continue t1(k,j)=a(k,j)-sum1 33 continue do 35 j=k+1,n sum2=0 do 34 i=k+1,n sum2=sum2+b(i,k)*a(j,i) 34 continue t2(j,k)=a(j,k)-sum2 35 continue do 81 i=1,n do 80 j=1,i a(j,i)=dconjg(a(i,j)) 80 continue 81 continue sum01=0 sum02=0 sum03=0 do 37 j=k+1,n sum01=sum01+dconjg(b(j,k))*a(j,k) sum02=sum02+ b(j,k) *a(k,j) do 36 l=k+1,n sum03=sum03+dconjg(b(l,k))*b(j,k)*a(l,j) 36 continue 37 continue a(k,k)=a(k,k)-sum01-sum02+sum03 do 38 j=1,k-1 a(k,j)=t1(k,j) 38 continue do 39 j=k+1,n a(j,k)=t2(j,k) 39 continue 40 continue C *************************************************************** C PRE-MULTIPLY AND POST-MULTIPLY (A) WITH (D)**(-1/2) C *************************************************************** DO 42 I=1,N A(I,I) = A(I,I) / W(I,1) T=B(I,I) DO 41 J=1,I-1 A(I,J) = A(I,J) * B(J,J) * T 41 CONTINUE 42 CONTINUE C *************************************************************** C FIND EIGENVALUES AND EIGENVECTORS OF THE C TRANSFORMED MATRIX C *************************************************************** C CALL EIGCH( A, N, N1, NE, NV, EPS, WD, LW, E, V, IER ) C C *************************************************************** C BACK-TRANSFORM EIGENVECTORS. C *************************************************************** IF ( NV .EQ. 0 ) RETURN DO 55 J=1,NV DO 51 I=1,N V(I,J) = V(I,J) * B(I,I) 51 CONTINUE DO 53 K=1,N-1 T = V(K,J) DO 52 I=K+1,N V(I,J) = V(I,J) - T * B(I,K) 52 CONTINUE 53 CONTINUE 55 CONTINUE RETURN END SUBROUTINE EIGCH(A, N, N1, NE, NV, EPS, W, LW, E, V, IER) ************************************************************************ * EIGENVALUES AND EIGENVECTORS OF A HERMITE MATRIX * * BY HOUSEHOLDER-BISECTION-INVERSE ITERATION METHOD. * * PARAMETERS * * (1) A: 2-DIM. ARRAY CONTAINING A HERMITE MATRIX * * (2) N: ORDER OF THE MATRIX * * (3) N1: ROW SIZE OF THE 2-DIM. ARRAY (A) * * (4) NE: NUMBER OF NEEDED EIGENVALUES * * (5) NV: NUMBER OF NEEDED EIGENVECTORS * * (6) EPS: TOLERANCE FOR CONVERGENCE * * (7) W: 2-DIM. WORKING ARRAY * * (8) LW: 1-DIM. WORKING ARRAY * * (9) E: 1-DIM. ARRAY CONTAINING COMPUTED EIGENVALUES * * (10) V: 2-DIM. ARRAY CONTAINING COMPUTED EIGENVALUES * * (11) IER: ERROR CODE * * COPYRIGHT T. OGUNI JUNE 30 1989 VERSION 1.0 * ************************************************************************ IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 A, V, CR, CS INTEGER*4 SW DIMENSION A(N1,N), E(N), V(N1,N), W(N1,7), LW(N) C IF (N .LE. 1 .OR. NE .EQ. 0) THEN IER = 2 WRITE(*,*) '(SUBR. EIGCH) INVALID ARGUMENT.', N, NE RETURN ENDIF NEA = IABS(NE) NVA = IABS(NV) IF (N1 .LT. N .OR. N .LT. NEA .OR. NEA .LT. NVA) THEN IER = 2 WRITE(*,*) '(SUBR. EIGCH) INVALID ARGUMENT.', N1, N, NEA, NVA RETURN ENDIF IF (EPS .LT. 0.0) EPS = 1.0D-14 IF (N .EQ. 2) THEN W(1,1) = A(1,1) W(2,1) = A(2,2) W(1,2) = cdABS(A(2,1)) A(1,1) = A(2,1) / W(1,2) T = 0.5D0 * (W(1,1) + W(2,1)) R = W(1,1) * W(2,1) - W(1,2) ** 2 D = T * T - R Q = DABS(T) + DSQRT(D) IF (T .LT. 0.0D0) Q = - Q T = T * DBLE(NE) IF (T .GE. 0.0D0) THEN E(1) = Q IF (NEA .EQ. 2) E(2) = R / Q ELSE E(1) = R / Q IF (NEA .EQ. 2) E(2) = Q ENDIF ELSE C REDUCE TO TRIDIAGONAL FORM BY HOUSEHOLDER'S METHOD. DO 60 I=1,N 60 W(I,1) = A(I,I) DO 190 K=1,N-2 S = 0.0D0 DO 70 I=K+1,N 70 S = S + DREAL(A(I,K))**2 + DIMAG(A(I,K))**2 SR = DSQRT(S) T = cdABS(A(K+1,K)) W(K,2) = - SR IF (T .EQ. 0.0) THEN A(K,K) = 1.0D0 ELSE A(K,K) = A(K+1,K) / T ENDIF IF (S .NE. 0.0D0) THEN R = 1.0D0 / (S + T * SR) A(K+1,K) = A(K+1,K) + SR * A(K,K) DO 140 I=K+1,N CS = (0.0D0, 0.0D0) IF (I .NE. K+1) THEN DO 110 J=K+1,I-1 110 CS = CS + A(I,J) * A(J,K) ENDIF CS = CS + W(I,1) * A(I,K) IF (I .NE. N) THEN DO 130 J=I+1,N 130 CS = CS + DCONJG(A(J,I)) * A(J,K) ENDIF 140 A(I,I) = CS * R CS = (0.0D0, 0.0D0) DO 150 I=K+1,N 150 CS = CS + DCONJG(A(I,K)) * A(I,I) CS = 0.5D0 * R * CS DO 160 I=K+1,N 160 A(I,I) = A(I,I) - CS * A(I,K) DO 170 I=K+1,N 170 W(I,1) = W(I,1) - 2.0D0 * DREAL(A(I,K) * DCONJG(A(I,I))) DO 180 I=K+2,N DO 180 J=K+1,I-1 180 A(I,J)=A(I,J)-A(I,K)*DCONJG(A(J,J))-A(I,I)*DCONJG(A(J,K)) ENDIF 190 CONTINUE W(N-1,2) = cdABS(A(N,N-1)) IF (W(N-1,2) .EQ. 0.0D0) THEN A(N-1,N-1) = 0.0D0 ELSE A(N-1,N-1) = A(N,N-1) / W(N-1,2) ENDIF C COMPUTE EIGENVALUES BY BISECTION METHOD. R=DMAX1((DABS(W(1,1))+DABS(W(1,2))),(DABS(W(N-1,2))+DABS(W(N,1) * ))) DO 200 I=2,N-1 T = DABS(W(I-1,2)) + DABS(W(I,1)) + DABS(W(I,2)) IF (T .GT. R) R = T 200 CONTINUE EPS1 = R * 1.0D-14 EPS2 = R * EPS DO 210 I=1,N-1 210 W(I,3) = W(I,2) ** 2 IF (NE .LT. 0) R = - R F = R DO 220 I=1,NEA 220 E(I) = - R DO 300 K=1,NEA D = E(K) 230 CONTINUE T = 0.5D0 * (D + F) IF (DABS(D-T) .GT. EPS2 .AND. DABS(F-T) .GT. EPS2) THEN J = 0 I = 1 240 Q = W(I,1) - T 250 IF (Q .GE. 0.0D0) J = J + 1 IF (Q .EQ. 0.0D0) GO TO 260 I = I + 1 IF (I .GT. N) GO TO 270 Q = W(I,1) - T - W(I-1,3) / Q GO TO 250 260 I = I + 2 IF (I .LE. N) GO TO 240 270 IF (NE .LT. 0) J = N - J IF (J .LT. K) THEN F = T ELSE D = T M = MIN0(J,NEA) DO 290 I=K,M 290 E(I) = T ENDIF GO TO 230 ENDIF 300 E(K) = T ENDIF C COMPUTE EIGENVECTORS BY INVERSE ITERATION. IF (NV .EQ. 0) RETURN MM = 584287 DO 490 I=1,NVA DO 320 J=1,N W(J,3) = W(J,1) - E(I) W(J,4) = W(J,2) 320 W(J,7) = 1.0D0 SW = 0 C REDUCE TO TRIANGULAR FORM DO 340 J=1,N-1 IF (DABS(W(J,3)) .GE. DABS(W(J,2))) THEN IF (W(J,3) .EQ. 0.0D0) W(J,3) = 1.0D-30 W(J,6) = W(J,2) / W(J,3) LW(J) = 0 W(J+1,3) = W(J+1,3) - W(J,6) * W(J,4) W(J,5) = 0.0D0 ELSE W(J,6) = W(J,3) / W(J,2) LW(J) = 1 W(J,3) = W(J,2) T = W(J,4) W(J,4) = W(J+1,3) W(J,5) = W(J+1,4) W(J+1,3) = T - W(J,6) * W(J,4) W(J+1,4) = - W(J,6) * W(J,5) ENDIF 340 CONTINUE IF (W(N,3) .EQ. 0.0D0) W(N,3) = 1.0D-30 C BEGIN BACK SUBSTITUTION IF (I .NE. 1 .AND. DABS(E(I)-E(I-1)) .LT. EPS1) THEN C GENERATE RANDOM NUMBERS DO 350 J=1,N mm=mod(1229*mm+351750, 1664501) c MM = iand(lambda*mm,mask) c MM*48828125 c Call urand1(n,x,ir) 350 W(J,7) = DFLOAT(MM) * 0.4656613D-9 ENDIF 360 T = W(N,7) R = W(N-1,7) W(N,7) = T / W(N,3) W(N-1,7) = (R - W(N-1,4) * W(N,7)) / W(N-1,3) IF (N .NE. 2) THEN K = N - 2 400 T = W(K,7) W(K,7) = (T - W(K,4)*W(K+1,7) - W(K,5)*W(K+2,7)) / W(K,3) K = K - 1 IF (K .GT. 0) GO TO 400 ENDIF 440 IF (SW .NE. 0) THEN SW = 1 DO 460 J=1,N-1 IF (LW(J) .NE. 0) THEN W(J+1,7) = W(J+1,7) - W(J,6) * W(J,7) ELSE T = W(J,7) W(J,7) = W(J+1,7) W(J+1,7) = T - W(J,6) * W(J+1,7) ENDIF 460 CONTINUE GO TO 360 ENDIF DO 480 J=1,N 480 V(J,I) = W(J,7) 490 CONTINUE C BEGIN BACK TRANSFORM (1) CR = 1.0D0 DO 500 J=2,N CR = CR * A(J-1,J-1) DO 500 I=1,NVA 500 V(J,I) = V(J,I) * CR C BEGIN BACK TRANSFORM (2) IF (N .NE. 2) THEN DO 590 I=1,NVA DO 580 K=N-2,1,-1 CR = - A(K+1,K) * DCONJG(A(K,K)) * W(K,2) IF ( DREAL(CR) .NE. 0.0D0 .OR. DIMAG(CR) .NE. 0.0D0) THEN CR = 1.0D0 / CR CS = (0.D0, 0.0D0) DO 560 J=K+1,N 560 CS = CS + DCONJG(A(J,K)) * V(J,I) CR = CR * CS DO 570 J=K+1,N 570 V(J,I) = V(J,I) - CR * A(J,K) ENDIF 580 CONTINUE 590 CONTINUE ENDIF C NORMALIZE EIGENVECTORS C NORMALIZE AS MAXIMUM ELEMENT = 1 DO 620 I=1,NVA T = DABS(DREAL(V(1,I)))+DABS(DIMAG(V(1,I))) K = 1 DO 610 J=2,N R = DABS(DREAL(V(J,I))) + DABS(DIMAG(V(J,I))) IF (T .LT. R) THEN T = R K = J ENDIF 610 CONTINUE CR = 1.0D0 / V(K,I) DO 620 J=1,N 620 V(J,I) = V(J,I) * CR IF (NV .LT. 0) RETURN C ORTHOGONALIZE AS NORM=1 DO 680 I=1,NVA IF (I.EQ.1 .OR. DABS(E(I)-E(I-1)).GT.EPS) THEN M = I ELSE DO 640 J=M,I-1 CS = (0.0D0, 0.0D0) DO 630 K=1,N 630 CS = CS + DCONJG(V(K,J)) * V(K,I) DO 640 K=1,N 640 V(K,I) = V(K,I) - CS * V(K,J) ENDIF C NORMALIZE AS NORM=1 S = 0.0D0 DO 670 J=1,N 670 S = S + DREAL(V(J,I))**2 + DIMAG(V(J,I))**2 T = DSQRT(1.0D0 / S) DO 680 J=1,N 680 V(J,I) = V(J,I) * T RETURN END