PROGRAM TryDiago C experimental matrix diagonalization IMPLICIT none integer*4 N,i,IERR,nb,ndim,ib,j,k,iter,it0,m real*8 err,eold,elim real*8,allocatable::a(:,:),C(:,:),E(:),B(:,:),AC(:,:),CB(:,:), 1EB(:),CNEW(:,:) integer*4,allocatable::db(:),bs(:) C open(8,file='A.SCR',form='unformatted') read(8)N write(6,*)'N = ',N allocate(a(N,N),C(N,N),E(N)) write(6,*)'allocated' read(8)a close(8) write(6,*)'read' CALL TRED12(N,A,C,E,2,IERR) write(6,*)'diagonalized - direct, IERR: ',IERR open(8,file='C.SCR.TXT') do 1 i=1,N 1 write(8,800)i,E(I) 800 format(i12,F12.6) close(8) write(6,600)E(1),E(2),E(N-1),E(N) 600 format(4g12.3) write(6,*)'written' c block-iterative diagonalization: write(6,*)'Number of blocks:' read(5,*)nb allocate(db(nb),bs(nb)) ndim=0 do 2 i=1,nb-1 db(i)=N/nb bs(i)=ndim+1 2 ndim=ndim+db(i) bs(nb)=ndim+1 db(nb)=N-ndim write(6,*)'block dimension start end' do 3 i=1,nb 3 write(6,*)i,db(i),bs(i),bs(i)+db(i)-1 c initial coefficients: do 4 i=1,N do 4 j=1,N 4 c(i,j)=0.0d0 do 5 j=1,N 5 c(j,j)=1.0d0 eold=1.0d12 elim=0.000001d0 iter=0 it0=100 1000 iter=iter+1 c diagonalizations in each block: do 6 ib=1,nb m=db(ib) c form matrix: allocate(B(m,m),AC(N,m),CB(m,m),EB(m),CNEW(N,m)) do 7 i=1,N do 7 j=1,m AC(i,j)=0.0d0 do 7 k=1,N 7 AC(i,j)=AC(i,j)+A(i,k)*c(k,bs(ib)+j-1) do 8 i=1,m do 8 j=1,m B(i,j)=0.0d0 do 8 k=1,N 8 B(i,j)=B(i,j)+c(k,bs(ib)+i-1)*AC(k,j) CALL TRED12(m,B,CB,EB,2,IERR) write(6,601)ib,EB(1),EB(2),EB(m-1),EB(m) c new vectors do 9 i=1,N do 9 j=1,m CNEW(i,j)=0.0d0 do 9 k=1,m 9 CNEW(i,j)=CNEW(i,j)+c(i,bs(ib)+k-1)*CB(k,j) c replace the old vectors: do 10 j=1,m do 10 i=1,N 10 c(i,bs(ib)+j-1)=CNEW(i,j) c replace energies do 11 j=1,m 11 E(bs(ib)+j-1)=EB(j) 6 deallocate(B,AC,CB,EB,CNEW) c reorder energies and eigenvectors c call reporder(N,E,C) call reporder2(N,E,C,nb,bs) err=dabs(E(1)-eold) write(6,601)iter,E(1),E(2),E(N-1),E(N),err,elim 601 format(i5,6g12.3) eold=E(1) c orthonormalize: c call ortog(N,C) if(err.gt.elim.and.iter.le.it0)goto 1000 end subroutine ortog(N,C) implicit none integer*4 N,i,j,k real*8 C(N,N),sij,norm do 2 i=1,N do 3 j=1,i-1 sij=0.0d0 do 4 k=1,N 4 sij=sij+C(k,i)*C(k,j) do 3 k=1,N 3 C(k,i)=C(k,i)-sij*C(k,j) norm=0.0d0 do 5 k=1,N 5 norm=norm+C(k,i)**2 norm=1.0d0/dsqrt(norm) do 2 k=1,N 2 C(k,i)=C(k,i)*norm return end subroutine reporder2(N,E,C,nb,ns) implicit none integer*4 N,i,ib,k,nb,ns(*) real*8 E(*),C(N,N),el real*8 ,allocatable::cl(:) allocate(cl(N)) c from each block one down and rest up do 1 ib=2,nb c lowest eigenvalue and eigenvector in this block: el=E(ns(ib)) do 10 k=1,N 10 cl(k)=C(k,ns(ib)) do 11 i=ns(ib),ib+1,-1 E(i)=E(i-1) do 11 k=1,N 11 C(k,i)=C(k,i-1) E(ib)=el do 1 k=1,N 1 C(k,ib)=cl(k) return end subroutine reporder(N,E,C) implicit none integer*4 N,i,j,ji,k real*8 E(*),C(N,N),t do 2 i=1,N ji=i do 3 j=i+1,N 3 if(E(j).lt.E(ji))ji=j if(ji.ne.i)then t=E(i) E(i)=E(ji) E(ji)=t do 4 k=1,N t=C(k,i) C(k,i)=C(k,ji) 4 C(k,ji)=t endif 2 continue return end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(*) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(N,N),D(*),E(*) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 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 IF (ICODE .NE. 2) GO TO 300 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 END