This is the source code of DMRG for Ising Model
This program runs with NEC's ASL library
(Please contact me if you wish to know about ASL.)
::Infinite System Method::
C123456789012345678901234567890123456789012345678901234567890123456789012
C ----- -----
C ----- Density Matrix Renormalization Group for a 2-state -----
C ----- symmetric IRF model. -----
C ----- '95/6/11, Ver.1.0 by T.Nishino -----
C ----- -----
C ----- Transfer Matrix: -----
C ----- @----O---O----@ O:2-state spin -----
C ----- TM = | TL | W | TR | -----
C ----- @----O---O----@ @:m-state spin -----
C ----- -----
C ----- C---D -----
C ----- | W | = Boltsmann Weight for the spin configuration -----
C ----- A---B A,B,C,D = +(-) 1. -----
C ----- -----
C
C +---- Parameters: MXDM ( Max. dim. of "V" ) -----------------+
C | MXTM ( Max. dim. of "TM" ) |
C +------------ MXRT ( Max. num. of the iteration ) -------+
PARAMETER( MXDM=1000000, MXTM=1000, MXRT=2400 )
INTEGER MX, N, M, L, LX(0:2,0:20)
REAL*8 R, TL(0:MXTM,0:MXTM), CV(0:MXRT,0:5)
REAL*8 T, DT, TR(0:MXTM,0:MXTM), LNZ(0:MXRT,0:2)
REAL*8 T0, T1, DM(0:MXTM,0:MXTM), ENG(0:MXRT,0:2)
REAL*8 T2, ER, LM(0:MXTM), PS(0:MXDM,0:3)
REAL*8 FL, W(0:3,0:3)
C
CALL GETPRM( L, MX, T1, T2, DT )
C
DO 90 K = 0, 4
T0 = T1+(T2-T1)*K/ 5.0D0
DO 20 J = 0, 1
WRITE(6,*)
T = T0 + DT*(J*2.0D0-1.0D0)/2.0D0
R = 0.0D0
CALL INITTM( M, W, T, TL, TR )
CALL VECINI( 4*M*M, PS(0,1))
I =-1
10 CONTINUE
I = I + 1
CALL NORMTM( M, R, TL, TR )
CALL NORMAL( 4*M*M, PS(0,1))
CALL LANZOS( M, W, LM, TL, TR, PS )
CALL GETWWW( M, W,ENG(I,J),PS(0,1))
LNZ(I,J)= R + LOG(-LM(0) )
CALL GETDML( M, DM, PS(0,1))
CALL DIAGAL( M, DM, LM )
CALL SELECT( M, N, MX, LM, ER )
CALL ROTATL( M, N, DM, TL, W )
M = N
CALL TCOPY( M, TL, TR )
C
IF( I.LT.12 ) GOTO 10
IF( J.NE.0 .AND. I.LT.LX(J-1,K) ) GOTO 10
FL = ABS( LNZ(I,J) -2*LNZ(I-1,J) +LNZ(I-2,J) )
IF( I.LT.L .AND. FL.GT.DT/10**11 ) GOTO 10
LX(J,K)= I
20 CONTINUE
CALL GETCV( L,T0, DT,ENG, ER, LX(0,K), CV(0,K) )
90 CONTINUE
CALL SUMMRY( CV, LX )
C
STOP
END
CCC
CC
C
SUBROUTINE GETPRM( L, MX, T1, T2, DT )
INTEGER L, MX
REAL*8 T1, T2, DT
C Write(6,*) 'Input System Size L'
C Read(5,*) L
L = 256
C Write(6,*) 'Input Min. Temperature T1'
C Read(5,*) T1
T1 = 3.0D0
C Write(6,*) 'Input Max. Temperature T2'
C Read(5,*) T2
T2 = 4.0D0
C Write(6,*) 'Input Max. Dim. 128 > mx >= m, '
C Read(5,*) MX
MX = 40
DT = 0.00010D0
Write(6,100) L, MX, T1, T2, DT
100 FORMAT(1H , 2I4, 3F15.7 )
RETURN
END
CCC
CC
C
SUBROUTINE INITTM( M, W, T, TL, TR )
PARAMETER( MXTM=1000 )
INTEGER M, A, B, C, D
REAL*8 W(0:3,0:3), T, TL(0:MXTM,0:MXTM), TR(0:MXTM,0:MXTM)
C +---------------------------------------------------------+
C | Boltsmann Weight for the Triangular Lattice Ising Model |
C +---------------------------------------------------------+
M = 2
DO 10 I = 0, 3
A =( I/2 )*2 -1
B = MOD(I,2) *2 -1
DO 10 J = 0, 3
C =( J/2 )*2 -1
D = MOD(J,2) *2 -1
W(J,I)= EXP( ((A*B+C*D+A*C+B*D)/2.0D0) / T )
TL(J,I)= EXP( ((A*B+C*D+A*C+B*D)/2.0D0) / T )
TR(J,I)= EXP( ((A*B+C*D+A*C+B*D)/2.0D0) / T )
10 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE NORMTM( M, R, TL, TR )
PARAMETER( MXTM=1000 )
INTEGER M
REAL*8 Q, R, TL(0:MXTM,0:MXTM), TR(0:MXTM,0:MXTM)
C
Q = ABS( TL(0,0) )
DO 10 J = 0, 2*M-1
DO 10 I = 0, 2*M-1
IF( Q .LT. ABS( TL(I,J) ) ) Q = ABS( TL(I,J) )
10 CONTINUE
DO 20 J = 0, 2*M-1
DO 20 I = 0, 2*M-1
TL(I,J)= TL(I,J) / Q
20 CONTINUE
R = R + LOG(Q)
C
Q = ABS( TR(0,0) )
DO 30 J = 0, 2*M-1
DO 30 I = 0, 2*M-1
IF( Q .LT. ABS( TR(I,J) ) ) Q = ABS( TR(I,J) )
30 CONTINUE
DO 40 J = 0, 2*M-1
DO 40 I = 0, 2*M-1
TR(I,J)= TR(I,J) / Q
40 CONTINUE
R = R + LOG(Q)
C
RETURN
END
CCC
CC
C
C +---------------------------------------------------+
C | A Very Rapid Multiplcation Algorithm by S.R.White |
C +---------------------------------------------------+
SUBROUTINE MULT( M, W, TL, TR, V1, V0 )
PARAMETER(MXTM=1000, MXDM=1000000, MXQM=2000000 )
INTEGER M, A, B, C, D, E, F, G, H
REAL*8 TL(0:MXTM,0:MXTM), V1(0:MXDM), X(0:MXQM)
REAL*8 W(0:3,0:3), TR(0:MXTM,0:MXTM), V0(0:MXDM), Y(0:MXQM)
REAL*8 WW(0:MXTM*2,0:3)
C
CALL VEC0( 8*M*M, X )
DO 10 A = 0, 2*M-1
DO 10 C = 0, 1
DO 10 D = 0, M-1
DO 10 H = 0, 2*M-1
X((A*2+C)*2*M+H)= X((A*2+C)*2*M+H)
& + TR(H,C*M+D)*V0((A*2+C)*M+D)
10 CONTINUE
C
DO 15 J = 0, 3
DO 15 I = 0, M-1
WW( I,J)= W(0,J)
WW( M+I,J)= W(1,J)
WW(2*M+1,J)= W(2,J)
WW(3*M+1,J)= W(3,J)
15 CONTINUE
C
CALL VEC0( 8*M*M, Y )
DO 20 A = 0, M-1
DO 20 B = 0, 1
DO 20 C = 0, 1
DO 20 F = 0, 1
DO 20 G = 0, 1
DO 20 H = 0, M-1
Y((((A*2+B)*2+F)*2+G)*M+H)= Y((((A*2+B)*2+F)*2+G)*M+H)
& + W(F*2+G,B*2+C)* X((((A*2+B)*2+C)*2+G)*M+H)
20 CONTINUE
C
CALL VEC0( 4*M*M, V1 )
DO 30 A = 0, M-1
DO 30 B = 0, 1
DO 30 E = 0, M-1
DO 30 F = 0, 1
DO 30 H = 0, 2*M-1
V1((E*2+F)*2*M+H)= V1((E*2+F)*2*M+H)
& - TL(E*2+F,A*2+B)*Y(((A*2+B)*2+F)*2*M+H)
30 CONTINUE
C
RETURN
END
CCC
CC
C
SUBROUTINE GETWWW( M, W, ENG, PS )
INTEGER M
REAL*8 W(0:3,0:3), ENG, PS(0:*)
C
ENG = 0.0D0
DO 10 I = 0, 2*M-1
DO 10 J = 0, 2*M-1
ENG = ENG + PS(I*2*M+J)*PS(I*2*M+J)
& *(2*MOD(I,2)-1)*(2*(J/M)-1)
10 CONTINUE
C
RETURN
END
CCC
CC
C
SUBROUTINE GETDML( M, DM, PS )
PARAMETER( MXTM=1000 )
INTEGER M
REAL*8 DM(0:MXTM,0:MXTM), PS(0:*)
DO 10 J = 0, 2*M-1
DO 10 I = 0, 2*M-1
DM(I,J)= 0.0D0
DO 10 K = 0, 2*M-1
DM(I,J)= DM(I,J) + PS(I*2*M+K) * PS(J*2*M+K)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE GETDMR( M, DM, PS )
PARAMETER( MXTM=1000 )
INTEGER M
REAL*8 DM(0:MXTM,0:MXTM), PS(0:*)
DO 10 J = 0, 2*M-1
DO 10 I = 0, 2*M-1
DM(I,J)= 0.0D0
10 CONTINUE
DO 20 K = 0, 2*M-1
DO 20 J = 0, 2*M-1
DO 20 I = 0, 2*M-1
DM(I,J)= DM(I,J) + PS(K*2*M+I) * PS(K*2*M+J)
20 CONTINUE
RETURN
END
CCC
SUBROUTINE DIAGAL( M, DM, LM )
PARAMETER( MXTM=1000 )
INTEGER M, LNA, IERR
REAL*8 DM(0:MXTM,0:MXTM), LM(0:*), W1(0:MXTM)
LNA = MXTM+1
CALL DCSMAA( DM, LNA, 2*M, LM, W1, IERR )
IF( IERR .NE. 0 ) THEN
WRITE(6,*) 'IERR = ', IERR,' in DIAGAL'
STOP
END IF
RETURN
END
CCC
CC
C
SUBROUTINE SELECT( M, N, MX, LM, ER )
INTEGER M, N, MX
REAL*8 R, LM(0:*), ER
C
C Write(6,*)
C Write(6,100) (LM(I),I=2*M-1,0,-1)
100 FORMAT(1H ,' Singular Values ',E15.7)
C
R = 1.1
N = 2*M
IF( 2*M .LE. MX ) THEN
ER = 0.0D0
RETURN
END IF
C ----------- Precise Cut-Off -----------
10 CONTINUE
N = N - 1
IF( N .EQ. 1 ) RETURN
IF( N .LE. MX .AND. LM(2*M-N)/LM(2*M-N-1) .GT. R ) THEN
C Write(6,*) N,' States are kept from ', 2*M,': ', LM(2*M-N)
ER = LM(2*M-N)
RETURN
END IF
GOTO 10
C
END
CCC
CC
C
SUBROUTINE ROTATL( M, N, U, TL, W )
PARAMETER( MXTM=1000, MXTP=2000 )
INTEGER M, N, A, B, C, D, X, Y, Z, Q, BX, BZ
REAL*8 U(0:MXTM,0:MXTM), TL(0:MXTM,0:MXTM)
REAL*8 TP(0:MXTP,0:MXTP), W(0:3,0:3)
C
DO 10 I = 0, 2*N-1
DO 10 J = 0, 4*M-1
TP(J,I)= 0.0D0
10 CONTINUE
C
DO 20 X = 0, N-1
BX = 2*M-1-X
DO 20 Y = 0, 1
DO 20 A = 0, M-1
DO 20 B = 0, 1
DO 20 D = 0, 1
DO 20 Q = 0, 1
DO 20 C = 0, M-1
TP(C*4+D*2+Q,X*2+Y) = TP(C*4+D*2+Q,X*2+Y)
& + TL(C*2+D,A*2+B) * W(D*2+Q,B*2+Y) * U(A*2+B,BX)
20 CONTINUE
C
DO 30 I = 0, 2*N-1
DO 30 J = 0, 2*N-1
TL(J,I)= 0.0D0
30 CONTINUE
C
DO 40 X = 0, N-1
DO 40 Y = 0, 1
DO 40 Z = 0, N-1
BZ = 2*M-1-Z
DO 40 D = 0, 1
DO 40 Q = 0, 1
DO 40 C = 0, M-1
TL(Z*2+Q,X*2+Y) = TL(Z*2+Q,X*2+Y)
& + TP(C*4+D*2+Q,X*2+Y) * U(C*2+D,BZ)
40 CONTINUE
C
RETURN
END
CCC
SUBROUTINE ROTATR( M, N, U, TR, W )
PARAMETER( MXTM=1000, MXTP=2000 )
INTEGER M, N, A, B, C, D, X, Y, Z, Q, BY, BQ
REAL*8 U(0:MXTM,0:MXTM), W(0:3,0:3)
REAL*8 TR(0:MXTM,0:MXTM), TP(0:MXTP,0:MXTP)
C
DO 10 I = 0, 2*N-1
DO 10 J = 0, 4*M-1
TP(J,I)= 0.0D0
10 CONTINUE
C
DO 20 X = 0, 1
DO 20 Y = 0, N-1
BY = 2*M-1-Y
DO 20 A = 0, 1
DO 20 B = 0, M-1
DO 20 Z = 0, 1
DO 20 C = 0, 1
DO 20 D = 0, M-1
TP(Z*2*M+C*M+D,X*N+Y) = TP(Z*2*M+C*M+D,X*N+Y)
& + W(Z*2+C,X*2+A) * TR(C*M+D,A*M+B) * U(A*M+B,BY)
20 CONTINUE
C
DO 30 I = 0, 2*N-1
DO 30 J = 0, 2*N-1
TR(J,I)= 0.0D0
30 CONTINUE
C
DO 40 X = 0, 1
DO 40 Y = 0, N-1
DO 40 Z = 0, 1
DO 40 Q = 0, N-1
BQ = 2*M-1-Q
DO 40 C = 0, 1
DO 40 D = 0, M-1
TR(Z*N+Q,X*N+Y) = TR(Z*N+Q,X*N+Y)
& + TP(Z*2*M+C*M+D,X*N+Y) * U(C*M+D,BQ)
40 CONTINUE
C
RETURN
END
CCC
SUBROUTINE TCOPY( M, TL, TR )
PARAMETER( MXTM=1000 )
INTEGER M
REAL*8 TL(0:MXTM,0:MXTM), TR(0:MXTM,0:MXTM)
C
DO 10 I = 0, 1
DO 10 J = 0, M-1
DO 10 K = 0, 1
DO 10 L = 0, M-1
TR(K*M+L,I*M+J)=TL(2*L+K,2*J+I)
10 CONTINUE
C
RETURN
END
CCC
CC
C
SUBROUTINE DUMP( N, T )
PARAMETER( MXTM=1000 )
INTEGER N
REAL*8 T(0:MXTM,0:MXTM)
Write(6,*)
DO 10 I = 0, 2*N-1
WRITE(6,100) (T(I,J),J=0,2*N-1)
100 FORMAT(1H ,8F9.5)
10 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE GETCV( L, T0, DT, ENG, ER, LX, CV )
PARAMETER( MXRT=2400 )
INTEGER L,LX(0:2), LXM
REAL*8 T0, DT, ENG(0:MXRT,0:2), CV(0:MXRT)
LXM = MIN(LX(0),LX(1))
DO 10 I = 0, L
IF( I .LE. LXM ) THEN
CV(I)=(ENG(I,0)-ENG(I,1)) *2.0D0 / DT
ELSE
CV(I)= 0.0D0
END IF
10 CONTINUE
WRITE(6,100) T0, ER, CV(8), CV(16), CV(32), CV(64), CV(LXM)
LX(0)= LXM
100 FORMAT(1H ,F5.3, E10.2, 7F9.5)
RETURN
END
CCC
SUBROUTINE SUMMRY( CV, LX )
PARAMETER( MXRT=2400 )
INTEGER LX(0:2,0:20)
REAL*8 CV(0:MXRT,0:5), DP(0:25,0:5)
C
DO 10 I = 0, 19
IF( LX(0,0)-I .GE. 0 ) DP(I,0)= CV(LX(0,0)-I,0)
IF( LX(0,1)-I .GE. 0 ) DP(I,1)= CV(LX(0,1)-I,1)
IF( LX(0,2)-I .GE. 0 ) DP(I,2)= CV(LX(0,2)-I,2)
IF( LX(0,3)-I .GE. 0 ) DP(I,3)= CV(LX(0,3)-I,3)
IF( LX(0,4)-I .GE. 0 ) DP(I,4)= CV(LX(0,4)-I,4)
C IF( LX(0,5)-I .GE. 0 ) DP(I,5)= CV(LX(0,5)-I,5)
10 CONTINUE
C
Write(6,*)
DO 20 I = 0, 19
Write(6,100) I, (DP(I,J),J=0,4)
100 FORMAT(1H , I4, 5F14.10 )
20 CONTINUE
C
Write(6,*)
Write(6,*) LX(0,0),LX(0,1),LX(0,2),LX(0,3),LX(0,4)
C
RETURN
END
CCC
CC-------------------------------------------------
C Vector Processors For the Lanczos Diagonalization
C -------------------------------------------------
SUBROUTINE INPRO( NV, ID, IS, PRD, V)
PARAMETER( MXDM=1000000 )
INTEGER NV, ID, IS
REAL*8 PRD, V(0:MXDM,0:3)
PRD = 0.0D0
DO 10 I = 0, NV-1
PRD = PRD + V(I,ID)*V(I,IS)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE ADDER( NV, ID, IS1, IS2, A, B, V)
PARAMETER( MXDM=1000000 )
INTEGER NV, ID, IS1, IS2
REAL*8 A, B, V(0:MXDM,0:3)
DO 10 I = 0, NV-1
V(I,ID)= A*V(I,IS1) + B*V(I,IS2)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE DIVER( NV, ID, IS, D, V)
PARAMETER( MXDM=1000000 )
INTEGER NV, ID, IS
REAL*8 D, V(0:MXDM,0:3)
DO 10 I = 0, NV-1
V(I,ID)= V(I,IS) / D
10 CONTINUE
RETURN
END
CCC
SUBROUTINE VEC0( NV, V)
INTEGER NV
REAL*8 V(0:*)
DO 10 J = 0, NV-1
V(J)= 0.0D0
10 CONTINUE
RETURN
END
CCC
SUBROUTINE IVEC0( NV,IV)
INTEGER NV
INTEGER IV(0:*)
DO 10 J = 0, NV-1
IV(J)= 0.0D0
10 CONTINUE
RETURN
END
CCC
SUBROUTINE VMOVE( NV, ID, IS, V )
PARAMETER( MXDM=1000000 )
INTEGER NV, ID, IS
REAL*8 V(0:MXDM,0:3)
DO 10 I = 0, NV-1
V(I,ID)= V(I,IS)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE VECCPY( NV, A, B )
INTEGER NV
REAL*8 A(0:*), B(0:*)
DO 10 I = 0, NV-1
B(I)= A(I)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE VECINI( NV, V )
INTEGER NV
REAL*8 V(0:*)
DO 10 I = 0,NV-1
V(I)=( DBLE(IRAND()) + 00.1 )/SQRT( DBLE(NV) )
10 CONTINUE
CALL NORMAL( NV, V )
RETURN
END
CCC
SUBROUTINE NORMAL( NV, V )
INTEGER NV
REAL*8 PRD, V(0:*)
CALL INPRO( NV, 0, 0, PRD, V )
PRD = SQRT(PRD)
CALL DIVER( NV, 0, 0, PRD, V )
RETURN
END
CCC
CC
C
SUBROUTINE LANZOS( M, W, ENG, TL, TR, VEC )
PARAMETER( MXDM=1000000, MXTM=1000 )
INTEGER MLS,LS,NV,LK,M
REAL*8 TL(0:MXTM,0:MXTM), TR(0:MXTM,0:MXTM), W(0:3,0:3)
REAL*8 EPS, ENG, ALP(0:999), TE(0:999),VEC(0:MXDM,0:3)
REAL*8 E2, E2P, BET(0:999), ME(0:999,0:4)
C
IF( LV.GE.3 ) WRITE(6,*) '-vvvv- Sub. LANcZOS -vvvvv-'
C +----------------------------------------+
C | MLS - Maximum Lanczos Step 999 |
C | EPS - Eigenvalue Convergence Parameter |
C +----------------------------------------+
LV = 0
NV = 4*M*M
MLS = 100
EPS = 1.0E-8
C ++++++ Initialization Steps ++++++
E2 = 0.0D0
CALL VMOVE( NV, 0, 1, VEC )
CALL VMOVE( NV, 3, 1, VEC )
CALL VEC0( NV, VEC(0,1) )
CALL VEC0( NV, VEC(0,2) )
C
CALL MULT( M, W, TL, TR, VEC(0,1), VEC(0,0) )
IF( LV .GE. 2 ) WRITE(6,*)
IF( LV .GE. 1 ) WRITE(6,*) '-+- Eigenvalues -+-'
C ++++++ LS = Lanczos Step, LK = 5*(LS/5) ++++++
LK = 0
LS = 0
10 CONTINUE
E2P = E2
C ++++++ Five Laanczos Steps ++++++
DO 20 J = 0, 4
CALL MACRO1( M, W, TL, TR, VEC,TE,ALP,BET,LK,0)
20 CONTINUE
C
C ++++++ Diagonalization of the Tri-Diagonal matrix ++++++
CALL MALANZ( LK, LV, ALP,BET,ENG, E2, ME)
LS = LS + 1
CC
CC ++++++ Safety Trap ++++++
IF( LS .EQ. 95 ) STOP
CC +++++++++++++++++++++++++
CC
C ++++++ Convergence Check ++++++
IF( LS .LE. MLS .AND. ABS(E2P-E2) .GE. EPS ) GOTO 10
LS = LS - 1
IF( LV .GE. 1 ) WRITE(6,100) ENG
100 FORMAT(1H ,' Energy Upper Bound = ',E30.20)
C ++++++ Lanczos Vector Creation ++++++
DO 30 I = 0,999
TE(I) = ME(I,0)
30 CONTINUE
C ++++++ Two-pass Method ++++++
CALL VMOVE( NV, 0, 3, VEC )
CALL ADDER( NV, 3, 3, 0, TE(1), 0.0D0, VEC )
CALL MULT( M, W, TL, TR, VEC(0,1), VEC(0,0) )
LK = 0
DO 50 I = 0, LS
DO 40 J = 0, 4
CALL MACRO1( M, W, TL, TR, VEC,TE,ALP,BET,LK,1)
40 CONTINUE
IF( LV .GE. 3 ) WRITE(6,*) LK, ' -th Second Pass'
50 CONTINUE
C ++++++ Inprove the Lanczos Vector via ++++++
C ++++++ the Inverse Iteration Method ++++++
CALL INVERS( LV, M, W, TL, TR, ENG, VEC )
IF( LV .GE. 1 ) WRITE(6,200) ENG
200 FORMAT(1H ,' Ground State Energy =', F20.15)
C
RETURN
END
CCC
CC
C
SUBROUTINE MACRO1( M,W,TL,TR,VEC,TE,ALP,BET,LK,ISW)
PARAMETER( MXDM=1000000, MXTM=1000)
INTEGER NV, LK, ISW,M
REAL*8 TL(0:MXTM,0:MXTM), TR(0:MXTM,0:MXTM), W(0:3,0:3)
REAL*8 PROD, BET(0:999), ALP(0:999), TE(0:999)
REAL*8 VEC(0:MXDM,0:3)
C
NV = 4*M*M
C ++++++ A Lanczos Step ++++++
CALL INPRO( NV, 0, 1, PROD, VEC )
ALP(LK) = PROD
PROD = -1.0 * PROD
CALL ADDER( NV, 2, 1, 0, 1.0D0, PROD, VEC )
CALL INPRO( NV, 2, 2, PROD, VEC )
PROD = SQRT(PROD)
BET(LK) = PROD
CALL DIVER( NV, 2, 2, PROD, VEC )
CALL MULT( M, W, TL, TR, VEC(0,1), VEC(0,2) )
PROD = -1.0 * PROD
CALL ADDER( NV, 1, 1, 0, 1.0D0, PROD, VEC )
CALL VMOVE( NV, 0, 2, VEC )
C
C ++++++ Second Path: Lanczos Vector Creation ++++++
IF( ISW .EQ. 1 ) CALL ADDER( NV, 3, 3, 0, 1.0D0, TE(LK), VEC)
LK = LK + 1
C
RETURN
END
CCC
CC
C
SUBROUTINE MALANZ( LK, LV, ALP, BET, ENG, E2, VE )
INTEGER N, M, LK, LV, LNV, ISW, IW1(0:999),IERR
REAL*8 ENG,E2, ALP(0:999), VE(0:999,0:4)
REAL*8 EPS, BET(0:999), E(0:999), W1(0:5999)
C
N = LK
EPS = -1.0
M = 5
LNV = 1000
ISW = -1
CALL DCSTSS( ALP,N,BET,EPS,E,M,VE,LNV,ISW,IW1,W1,IERR )
IF( IERR .NE. 0 ) THEN
WRITE(6,100) IERR
100 FORMAT(1H ,' ERROR NO. in DCSTSS(ASL)',I10)
END IF
ENG = E(0)
E2 = E(1)
C
IF( LV .GE. 1 ) WRITE(6,200) (E(J),J = 0,2), LK
200 FORMAT(1H ,'***',3F20.15,'*** ',I3,'- Steps')
C
RETURN
END
CCC
CC
C
SUBROUTINE INVERS( LV, M, W, TL, TR, LMD, VEC)
PARAMETER( MXDM=1000000, MXTM=1000)
INTEGER NV,MLR, LR,LV, M
REAL*8 TL(0:MXTM,0:MXTM), TR(0:MXTM,0:MXTM), W(0:3,0:3)
REAL*8 DTP, EPS, DTR, LMD, VEC(0:MXDM,0:3)
C
IF( LV.GE.3 ) WRITE(6,*) '-vvvv- Sub. INVERS -vvvvv-'
C +-----------------------------------------+
C | MLR - Maximum Inverse Iteration # |
C | EPS - Eigenvector Convergence Parameter |
C +-----------------------------------------+
LR = 0
MLR = 20
EPS = 0.5E-11
NV = 4*M*M
C
CALL VMOVE( NV, 0, 3, VEC )
GOTO 20
10 CONTINUE
CALL CGMETH( LV, M, W, TL, TR, LMD,VEC )
20 CALL VMOVE( NV, 1, 0, VEC )
CALL INPRO( NV, 1, 1, DTP, VEC )
DTP = SQRT(DTP)
CALL DIVER( NV, 1, 1, DTP, VEC )
CALL MULT( M, W, TL, TR, VEC(0,2), VEC(0,1) )
CALL INPRO( NV, 2, 2, DTP, VEC )
DTP = SQRT(DTP)
CALL DIVER( NV, 2, 2, DTP, VEC )
IF( LMD .LT. 0 ) THEN
CALL ADDER( NV, 3, 2, 1, 1.0D0,1.0D0, VEC )
ELSE
CALL ADDER( NV, 3, 2, 1,-1.0D0,1.0D0, VEC )
END IF
CALL INPRO( NV, 3, 3, DTR, VEC )
IF( LV .GE. 1 ) WRITE(6,100) LMD, DTP, SQRT(DTR)
100 FORMAT(1H ,'Lmd=',E17.10,' Lmd2=',E17.10,' Difference=',E17.10)
C ++++++ Ground State Energy Inprovement ++++++
IF( LMD .GT. 0 .AND. DTP .LT. LMD ) LMD = DTP
IF( LMD .LT. 0 .AND. DTP .GT.-LMD ) LMD =-DTP
C
IF( ABS(LMD) .LE. EPS ) THEN
WRITE(6,*)'Eigen value ',LMD,' might be Accidentaly Zero.'
WRITE(6,*)'Please change chemical potential to get normal',
& ' result.'
IF( LV .GE. 3 ) THEN
WRITE(6,*)'.........Software Termination........'
STOP
END IF
END IF
C
LR = LR + 1
IF( LR .GE. MLR ) THEN
C WRITE(6,*) 'Eigen vector is not obtained in', MLR,
C & ' Inverse Iterations.'
C WRITE(6,*) 'Software Termination'
C STOP
RETURN
END IF
C ++++++ Convergence Check ++++++
IF( SQRT(DTR) .GE. EPS ) GOTO 10
C
RETURN
END
CCC
CC
C
SUBROUTINE CGMETH( LV, M, W, TL, TR, LMD,VEC )
PARAMETER( MXDM=1000000, MXTM=1000)
INTEGER NV, LR, LV, MLR, M
REAL*8 TL(0:MXTM,0:MXTM), TR(0:MXTM,0:MXTM), W(0:3,0:3)
REAL*8 DB, DR1, DR2, ALP, BET, EPS, DTP, LMD,VEC(0:MXDM,0:3)
C
C +-----------------------------------------+
C | MLR - Maximum Inverse Iteration # |
C | EPS - Eigenvector Convergence Parameter |
C +-----------------------------------------+
LR = 0
MLR = 20
EPS =1.0E-8
NV = 4*M*M
C (I)
CALL INPRO( NV, 1, 1, DB, VEC)
CALL VMOVE( NV, 2, 1, VEC)
DO 10 I = 0,NV-1
VEC(I,0)= 0.0D0
10 CONTINUE
C (II)
CALL INPRO( NV, 1, 1, DR1, VEC)
20 CALL MULT( M, W, TL, TR, VEC(0,3), VEC(0,2) )
DO 30 I = 0, NV-1
VEC(I,3)= VEC(I,3) - LMD*VEC(I,2)
30 CONTINUE
CALL INPRO( NV, 3, 2,DTP, VEC)
ALP=DR1/DTP
CALL ADDER( NV, 0, 0, 2, 1.0D0, ALP, VEC)
CALL ADDER( NV, 1, 1, 3, 1.0D0,-ALP, VEC)
CALL INPRO( NV, 1, 1, DR2, VEC)
C
IF( LV .GE. 3 ) WRITE(6,100) DR2 / DB, LR
100 FORMAT(1H ,' CG Conv. Parameter = ', E18.10,' -',I3,'-Step')
IF( DR2/DB .LE. EPS ) GOTO 40
BET = DR2 / DR1
DR1 = DR2
LR = LR + 1
CALL ADDER( NV, 2, 1, 2,1.0D0, BET, VEC)
IF( LR .GE. MLR ) THEN
IF(LV.GE.4) WRITE(6,*) 'Eigen vector is not obtained in ',
& MLR,' CG steps.'
GOTO 50
END IF
GOTO 20
40 CONTINUE
50 IF( LV .GE. 3 ) WRITE(6,*) LR,' Rotations in CG'
C
RETURN
END