Lanczos Diagonalization for S=1 Heisenberg Spin Ring
@
C123456789012345678901234567890123456789012345678901234567890123456789012 C C +-----------------------------------------------------+ C | | C | Diagonalization of S = 1 Heisenberg Model | C | | C | H = S*S + a(S*S)^2, Generalized S=1 chain... | C | | C | Bisection Search Method is used. | C | Translational Invariance and Parity are used. | C | Periodic Boundary Condition is assumed. | C | | C | * This program assumes "Long Integer Mode"; some * | C | * integers should be 64-bit integer. * | C | | C +-----------------------------------------------------+ C C Site# Config. # Reduced Config. C (Sz=0) (*About*) C 3: 7 3 C 4: 19 7 C 5: 51 11 C 6: 141 23* Dump mode C 7: 393 56* C 8: 1,107 138* C 9: 3,139 348* C 10: 8,953 895* C 11: 25,653 2,332* C 12: 73,789 6,149* C 13: 212,941 16,380* C 14: 616,227 44,016* C 15: 1,787,607 119,173* Debug Mode C 16: 5,196,627 324,789* Standard C 17: 15,134,931 890,290* C 18: 44,152,809 2,452,933* (3^ 9 = 19683) C 19: 128,996,853 6,789,308* C 20: 377,379,369 18,868,968* (3^10 = 59049) C 21: 1,105,350,729 52,635,749* (3^11 = 177147) C C +---------------------------------------------------+ C | Limits: NS =< 40. | C | Patch : if NS > 20, modify 59049 ( = 3^10). | C | if NS > 20, modify 9953 in Subroutine | C | "CRESTA" (>Largest Subspace for NS=10) | C | if NS > 15, modify 122,222 and 244,444 | C | (=Vector Dimensions) | C | if NS > 32 modify Subroutine "CRENUM" | C +---------------------------------------------------+ C INTEGER LV, NS, NST, BH(0:59049,0:2), ST(0:122222) INTEGER NK, NX, TZ, BL(0:59049,0:2), EXS(0:999) REAL*8 A, CS(0:40), SN(0:40), NRM(0:999) REAL*8 VEC(0:244444,0:3), DIA(0:122222) C C +---------------------------------------------------------------+ C | Debug Level: 0(None) 1(Normal) 2(Detail) 3(Dump) 4(Core Dump) | C +---------------------------------------------------------------+ LV = 1 C +------------------------+ C | System Size "NS" | C | NS Should be even!! | C | Total Sz "Tz" | C +------------------------+ NS = 8 DO 95 J = NS-2, 0, -1 TZ = J C C *** Parity is not used yet *** C +---------------------------------+ C | Coefficient of the (S*S)^2 term | C | A=1: Permutation Model | C | A=?: AKLT Model | C +---------------------------------+ A = 0.0D0 C WRITE(6,*) WRITE(6,*) WRITE(6,*) '*** Parameters:' WRITE(6,*) '*** Ns, Total Sz, Alpha =', NS, TZ, A C C +----------------------------+ C | Total Momentum = pi*(I/NS) | C +----------------------------+ DO 90 I = 0, NS/2 C-DB- DO 90 I = 0, 0 NK = I WRITE(6,*) 'Total Momentum = ', 2*NK, '/', NS, '*pi' CALL CRENUM( NS, BH,BL ) CALL CRESTA( NS,LV,TZ, NST,ST, BH,BL ) C IF( NST .EQ. 1 ) THEN WRITE(6,*) 'You Can Get Analytic Formula' GOTO 80 END IF C CALL ORTSET( NS,LV,NK, SN,CS ) CALL NORMAL( NS,LV, NX,NST,ST,EXS, CS,NRM, VEC ) CALL CREDIA( NS,LV, NST,ST, A,DIA ) CALL LANZOS( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC ) CALL VECDMP( NS,LV,NK, NST,ST, VEC ) CALL SZDIST( NS, NK, NST,ST, BH,BL, VEC ) 80 CONTINUE 90 CONTINUE 95 CONTINUE C STOP END CCC CC C SUBROUTINE CRENUM( NS, BH, BL ) INTEGER M, NS, T(0:20), BH(0:59049,0:2), BL(0:59049,0:2) INTEGER B0(0:2), B1(0:2), B2(0:2) INTEGER A0, A1, A2, A3, A4, A5, A6, A7, A8, A9 INTEGER AA, AB, AC, AD, AE, AF C C Bh Bl C [NS-1 .... NS/2] [NS/2-1 .... 0] C M = NS / 2 DO 10 I = 0, 20 IF( I .LT. M ) THEN T(I)= 1 ELSE T(I)= 0 END IF 10 CONTINUE C DO 20 I = 0, 2 B0(I) = MOD(I+1,2) - I/2 B1(I) = MOD(I ,2) B2(I) = I/2 20 CONTINUE C DO 30 I = 0, 3**M - 1 A0 = MOD(I , 3) A1 = MOD(I/3 , 3) A2 = MOD(I/3** 2, 3) A3 = MOD(I/3** 3, 3) A4 = MOD(I/3** 4, 3) A5 = MOD(I/3** 5, 3) A6 = MOD(I/3** 6, 3) A7 = MOD(I/3** 7, 3) A8 = MOD(I/3** 8, 3) A9 = MOD(I/3** 9, 3) AA = MOD(I/3**10, 3) AB = MOD(I/3**11, 3) AC = MOD(I/3**12, 3) AD = MOD(I/3**13, 3) AE = MOD(I/3**14, 3) AF = MOD(I/3**15, 3) BH(I,0) = T( 0)*B0(A0) +T( 1)*B0(A1) +T( 2)*B0(A2) +T( 3)*B0(A3) & + T( 4)*B0(A4) +T( 5)*B0(A5) +T( 6)*B0(A6) +T( 7)*B0(A7) & + T( 8)*B0(A8) +T( 9)*B0(A9) +T(10)*B0(AA) +T(11)*B0(AB) & + T(12)*B0(AC) +T(13)*B0(AD) +T(14)*B0(AE) +T(15)*B0(AF) BH(I,1) = T( 0)*B1(A0) +T( 1)*B1(A1) +T( 2)*B1(A2) +T( 3)*B1(A3) & + T( 4)*B1(A4) +T( 5)*B1(A5) +T( 6)*B1(A6) +T( 7)*B1(A7) & + T( 8)*B1(A8) +T( 9)*B1(A9) +T(10)*B1(AA) +T(11)*B1(AB) & + T(12)*B0(AC) +T(13)*B0(AD) +T(14)*B0(AE) +T(15)*B0(AF) BH(I,2) = T( 0)*B2(A0) +T( 1)*B2(A1) +T( 2)*B2(A2) +T( 3)*B2(A3) & + T( 4)*B2(A4) +T( 5)*B2(A5) +T( 6)*B2(A6) +T( 7)*B2(A7) & + T( 8)*B2(A8) +T( 9)*B2(A9) +T(10)*B2(AA) +T(11)*B2(AB) & + T(12)*B2(AC) +T(13)*B2(AD) +T(14)*B2(AE) +T(15)*B2(AF) 30 CONTINUE C M = NS - M DO 40 I = 0, 20 IF( I .LT. M ) THEN T(I)= 1 ELSE T(I)= 0 END IF 40 CONTINUE C DO 50 I = 0, 3**M - 1 A0 = MOD(I , 3) A1 = MOD(I/3 , 3) A2 = MOD(I/3** 2, 3) A3 = MOD(I/3** 3, 3) A4 = MOD(I/3** 4, 3) A5 = MOD(I/3** 5, 3) A6 = MOD(I/3** 6, 3) A7 = MOD(I/3** 7, 3) A8 = MOD(I/3** 8, 3) A9 = MOD(I/3** 9, 3) AA = MOD(I/3**10, 3) AB = MOD(I/3**11, 3) AC = MOD(I/3**12, 3) AD = MOD(I/3**13, 3) AE = MOD(I/3**14, 3) AF = MOD(I/3**15, 3) BL(I,0) = T( 0)*B0(A0) +T( 1)*B0(A1) +T( 2)*B0(A2) +T( 3)*B0(A3) & + T( 4)*B0(A4) +T( 5)*B0(A5) +T( 6)*B0(A6) +T( 7)*B0(A7) & + T( 8)*B0(A8) +T( 9)*B0(A9) +T(10)*B0(AA) +T(11)*B0(AB) & + T(12)*B0(AC) +T(13)*B0(AD) +T(14)*B0(AE) +T(15)*B0(AF) BL(I,1) = T( 0)*B1(A0) +T( 1)*B1(A1) +T( 2)*B1(A2) +T( 3)*B1(A3) & + T( 4)*B1(A4) +T( 5)*B1(A5) +T( 6)*B1(A6) +T( 7)*B1(A7) & + T( 8)*B1(A8) +T( 9)*B1(A9) +T(10)*B1(AA) +T(11)*B1(AB) & + T(12)*B1(AC) +T(13)*B1(AD) +T(14)*B1(AE) +T(15)*B1(AF) BL(I,2) = T( 0)*B2(A0) +T( 1)*B2(A1) +T( 2)*B2(A2) +T( 3)*B2(A3) & + T( 4)*B2(A4) +T( 5)*B2(A5) +T( 6)*B2(A6) +T( 7)*B2(A7) & + T( 8)*B2(A8) +T( 9)*B2(A9) +T(10)*B2(AA) +T(11)*B2(AB) & + T(12)*B2(AC) +T(13)*B2(AD) +T(14)*B2(AE) +T(15)*B2(AF) 50 CONTINUE C RETURN END CCC CC C SUBROUTINE CRESTA( NS, LV, TZ, NST, ST, BH, BL) INTEGER NL,LZ, M, NS, LV, TZ, NST, ST(0:122222) INTEGER BH(0:59049,0:2), PT(0:131441), FG(0:131441) INTEGER BL(0:59049,0:2), TP(0:131441), NLS(-10:10) INTEGER LS(0:9953,-10:10) C M = NS / 2 C --- Make Table --- DO 20 I = -M, M NL = 0 DO 10 J = 0, 3**M -1 LZ = BL(J,2) - BL(J,0) IF( LZ .EQ. I ) THEN LS(NL,I) = J NL = NL + 1 END IF 10 CONTINUE NLS(I)= NL 20 CONTINUE C NST = 0 DO 80 I = 0, 3**(NS-M) -1 LZ = TZ - BH(I,2) + BH(I,0) C --- Select Polalization --- IF( ABS(LZ) .LE. M ) THEN DO 30 J = 0, NLS(LZ) -1 PT(J)= I*3**M + LS(J,LZ) FG(J)= 0 30 CONTINUE C --- Copy States --- DO 40 K = 0, NLS(LZ) -1 TP(K)= PT(K) 40 CONTINUE C --- Mark States --- DO 60 K = 1, NS-1 DO 50 J = 0, NLS(LZ) -1 TP(J)= MOD(TP(J)*3,3**NS) + TP(J)/3**(NS-1) IF( PT(J) .GT. TP(J) ) FG(J) = 1 50 CONTINUE 60 CONTINUE C --- Select States --- DO 70 J = 0, NLS(LZ) -1 IF( FG(J) .EQ. 0 ) THEN ST(NST) = PT(J) NST = NST + 1 END IF 70 CONTINUE END IF 80 CONTINUE C IF( LV .GE. 1 ) WRITE(6,*) 'Total # of basis=', NST IF( LV .GE. 3 ) THEN DO 90 I = 0, NST -1 CALL STDISP( NS, ST(I), 0.0D0 ) 90 CONTINUE END IF C RETURN END CCC SUBROUTINE STDISP( NS, ST, F ) INTEGER TP, NS, ST, LW(0:40) REAL*8 F TP = ST DO 10 I = NS -1, 0, -1 LW(I)= MOD(TP,3) TP = TP/3 10 CONTINUE WRITE(6,100)ST, F, (LW(I),I=0,NS-1) 100 FORMAT(1H ,I14,' ---', F20.10,' ', 40I1) RETURN END CCC CC C SUBROUTINE ORTSET( NS, LV, NK, SN , CS ) INTEGER NS, LV, NK REAL*8 SN(0:40), CS(0:40) REAL*8 PI C IF( LV .GE. 2 ) WRITE(6,*) '--- Subroutine ORTSET ---' PI = 4.0D0*ATAN(1.0D0) C DO 10 I = 0, NS-1 SN(I)= SIN( 2.0D0*I*NK*PI/NS ) CS(I)= COS( 2.0D0*I*NK*PI/NS ) IF( LV .GE. 2 ) WRITE(6,100) SN(I), CS(I) 100 FORMAT(1H ,' SIN=', F10.5, ' COS=', F10.5 ) 10 CONTINUE C RETURN END CCC CC C SUBROUTINE NORMAL( NS, LV, NX, NST, ST, EXS, CS, NRM, VEC ) INTEGER RT, NS, LV, NX, NST, ST(0:122222),EXS(0:999) INTEGER TP(0:122222) REAL*8 F, CS(0:40), NRM(0:999), VEC(0:244444,0:3) C IF( LV .GE. 2 ) WRITE(6,*) '--- Subroutine NORMAL ---' C DO 10 I = 0, NST -1 VEC(I,0)= 0.0D0 TP(I )= ST(I) 10 CONTINUE C --- Mark Self Similar States --- DO 30 I = 1, NS -1 DO 20 J = 0, NST -1 TP(J)= MOD(TP(J)*3,3**NS) + TP(J)/3**(NS-1) IF( ST(J) .EQ. TP(J) ) THEN VEC(J,0) = 1.0D0 END IF 20 CONTINUE 30 CONTINUE C --- Pick up Self Similar States --- NX = 0 DO 40 I = 0, NST -1 IF( VEC(I,0) .GT. 1.0E-8 ) THEN EXS(NX)= I NX = NX + 1 END IF 40 CONTINUE C --- Normalize --- DO 60 I = 0, NX -1 RT = ST(EXS(I)) F = 1.0D0 DO 50 J = 1, NS -1 RT = RT/3**(NS-1) + MOD(RT*3,3**NS) IF( RT .EQ. ST(EXS(I)) ) F = F + CS(J) 50 CONTINUE IF( F .GT. 1.0E-8 ) THEN NRM(I) = 1.0D0 / SQRT( F ) ELSE NRM(I) = 0.0D0 END IF 60 CONTINUE C IF( NX .GT. 999) THEN WRITE(6,*)'Degenerate Basis # > ',999 do 99 i = 0, 100 call STDISP(NS,ST(EXS(i)),0.0D0) 99 continue WRITE(6,*)'...Increase 999 in the source list' STOP END IF C IF( LV .GE. 1 ) WRITE(6,200) NX 200 FORMAT(1H ,' # of Auto Correlated Basis=', I10 ) C IF( LV .GE. 2 ) THEN DO 70 I = 0, NX -1 CALL STDISP( NS, ST(EXS(I)), NRM(I) ) 70 CONTINUE END IF C RETURN END CCC CC C SUBROUTINE CREDIA( NS, LV, NST, ST, A, DIA ) INTEGER NS, LV, NST, ST(0:122222) REAL*8 AA(0:4), A,DIA(0:122222) C IF( LV .GE. 2 ) WRITE(6,*)'SUBROUTINE CREDIA' C AA(0)= A AA(1)= A AA(2)= A*2 AA(3)= A AA(4)= A C C --- Boundary Part --- DO 10 J = 0, NST -1 DIA(J)= ( MOD(ST(J),3)-1 )*( ST(J)/3**(NS-1)-1 ) & + AA( MOD(ST(J),3) + ST(J)/3**(NS-1) ) 10 CONTINUE C DO 30 I = 0, NS -2 DO 20 J = 0, NST -1 DIA(J)= DIA(J) & + ( MOD(ST(J)/3**I,3)-1 )*( MOD(ST(J)/3**(I+1),3)-1 ) & + AA( MOD(ST(J)/3**I,3) + MOD(ST(J)/3**(I+1),3) ) 20 CONTINUE 30 CONTINUE C IF( LV .GE. 3 ) THEN DO 40 I = 0, NST-1 WRITE(6,100) I, DIA(I) 40 CONTINUE 100 FORMAT(1H , I4, F15.7 ) END IF C RETURN END CCC CC C SUBROUTINE INPRO( NV, ID, IS, PROD, VEC) INTEGER NV, ID, IS REAL*8 PROD, VEC(0:244444,0:3) PROD = 0.0D0 DO 10 I = 0, NV-1 PROD = PROD + VEC(I,ID) * VEC(I,IS) 10 CONTINUE RETURN END CCC SUBROUTINE ADDER( NV, ID, IS1, IS2, A, B, VEC) INTEGER NV, ID, IS1, IS2 REAL*8 A, B, VEC(0:244444,0:3) DO 10 I = 0, NV-1 VEC(I,ID) = A*VEC(I,IS1) + B*VEC(I,IS2) 10 CONTINUE RETURN END CCC SUBROUTINE DIVER( NV, ID, IS, D,VEC) INTEGER NV, ID, IS REAL*8 D,VEC(0:244444,0:3) DO 10 I = 0, NV-1 VEC(I,ID) = VEC(I,IS) / D 10 CONTINUE RETURN END CCC SUBROUTINE VEC0( NV, VEC) INTEGER NV REAL*8 VEC(0:244444,0:3) DO 20 I = 0, 3 DO 10 J = 0, NV -1 VEC(J,I)= 0.0D0 10 CONTINUE 20 CONTINUE RETURN END CCC SUBROUTINE VMOVE( NV, ID, IS, VEC) INTEGER NV, ID, IS REAL*8 VEC(0:244444,0:3) DO 10 I = 0, NV-1 VEC(I,ID)= VEC(I,IS) 10 CONTINUE RETURN END CCC CC C SUBROUTINE DISPLN( NS, VAL ) INTEGER NS REAL*8 VAL(0:40) WRITE(6,100) ( VAL(I),I=0,NS-1 ) 100 FORMAT(1H ,8F9.6) RETURN END CCC CC C SUBROUTINE VECINI( NS, NK, NV, NX, EXS, NRM, VEC) INTEGER NS, NK, NV, NX, EXS(0:999) REAL*8 PROD,VEC(0:244444,0:3), NRM(0:999) C IF( NK.EQ.0 .OR. NS.EQ.NK*2 ) THEN CALL VEC0( NV, VEC ) DO 10 I = 0, NV-1 VEC(I,0)= (MOD(I,13)+MOD(I,19)+00.1) / SQRT(DBLE(NV)) 10 CONTINUE C *VDIR NODEP DO 20 I = 0, NX-1 VEC(EXS(I),0)= VEC(EXS(I),0) * NRM(I) 20 CONTINUE C CALL INPRO(NV,0,0,PROD,VEC) PROD = SQRT(PROD) CALL DIVER(NV,0,0,PROD,VEC) ELSE CALL VEC0( 2*NV, VEC ) DO 30 I = 0, NV-1 VEC(I ,0)= (MOD(I,13)+MOD(I,19)+00.1) / SQRT(DBLE(NV)) VEC(I+NV,0)= (MOD(I,17)+MOD(I,23)+00.1) / SQRT(DBLE(NV)) 30 CONTINUE C *VDIR NODEP DO 40 I = 0, NX-1 VEC(EXS(I) ,0)= VEC(EXS(I) ,0) * NRM(I) VEC(EXS(I)+NV,0)= VEC(EXS(I)+NV,0) * NRM(I) 40 CONTINUE C CALL INPRO(NV*2,0,0,PROD,VEC) PROD = SQRT(PROD) CALL DIVER(NV*2,0,0,PROD,VEC) END IF C RETURN END CCC CC C SUBROUTINE SCALE( ID, IS, NS, LV, NK,NST, NX,EXS,NRM,VEC) INTEGER ID, IS, NS, LV, NK,NST, NX,EXS(0:999) REAL*8 NRM(0:999), VEC(0:244444,0:3) C IF( LV .GE. 3 ) WRITE(6,*) '--- Subroutine SCALE ---' C IF( NK.EQ.0 .OR. NS.EQ.NK*2 ) THEN DO 10 I = 0, NST -1 VEC(I,ID)= 0.0D0 10 CONTINUE *VDIR NODEP DO 20 I = 0, NX -1 VEC(EXS(I),IS)= VEC(EXS(I),IS) * NRM(I) 20 CONTINUE ELSE DO 30 I = 0, NST -1 VEC(I ,ID)= 0.0D0 VEC(I+NST,ID)= 0.0D0 30 CONTINUE *VDIR NODEP DO 40 I = 0,NX -1 VEC(EXS(I) ,IS)= VEC(EXS(I) ,IS) * NRM(I) VEC(EXS(I)+NST,IS)= VEC(EXS(I)+NST,IS) * NRM(I) 40 CONTINUE END IF C RETURN END CCC CC C SUBROUTINE RESCAL( ID, IS, NS, LV, NK,NST, NX,EXS,NRM,VEC) INTEGER ID, IS, NS, LV, NK,NST, NX,EXS(0:999) REAL*8 NRM(0:999), VEC(0:244444,0:3) C IF( LV .GE. 3 ) WRITE(6,*) '--- Subroutine RESCAL ---' C IF( NK.EQ.0 .OR. NS.EQ.NK*2 ) THEN *VDIR NODEP DO 10 I = 0, NX -1 VEC(EXS(I),ID)= VEC(EXS(I),ID) * NRM(I) IF( NRM(I) .GT. 1.0E-8 ) THEN VEC(EXS(I),IS)= VEC(EXS(I),IS) / NRM(I) END IF 10 CONTINUE ELSE *VDIR NODEP DO 20 I = 0, NX -1 VEC(EXS(I) ,ID)= VEC(EXS(I) ,ID) * NRM(I) VEC(EXS(I)+NST,ID)= VEC(EXS(I)+NST,ID) * NRM(I) IF( NRM(I) .GT. 1.0E-8 ) THEN VEC(EXS(I) ,IS)= VEC(EXS(I) ,IS) / NRM(I) VEC(EXS(I)+NST,IS)= VEC(EXS(I)+NST,IS) / NRM(I) END IF 20 CONTINUE END IF C RETURN END CCC CC C SUBROUTINE MULT(ID,IS,NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC) INTEGER IP, NW,ID,IS,NS,LV,NK,NX,NST,ST(0:122222) INTEGER VLE, IX,EXS(0:999) INTEGER V0(0:141072), V2(0:141072), V4(0:141072), V5(0:141072) REAL*8 F0(0:141072), F1(0:141072), SN(0:40), A, CS(0:40) REAL*8 DIA(0:122222), VEC(0:244444,0:3), NRM(0:999) C IF( LV .GE. 3 ) WRITE(6,*) '--- Subroutine MULT ---' C --- Scale with VEC clear --- CALL SCALE( ID, IS, NS, LV, NK,NST, NX, EXS, NRM, VEC) C DO 90 I = 0, NST-1, 141072 IX = I VLE = MIN(141072,NST) C --- State Copy --- DO 10 J = 0, VLE -1 V0(J) = ST(IX+J) 10 CONTINUE C --- Loop: (Shift)^N --- DO 80 N = 0, NS-1 C --- Hamiltonian Decomposition --- IF( LV .GE. 3 ) WRITE(6,*) 'Rotation Level =', N DO 60 K = 0, NS-1 IP = K C --- S+S- Part --- CALL NNSPSM( NS, IP, VLE, V0, V2, A, F0 ) CALL COMPRS( ST(NST-1), IX, VLE, NW, V2, F0, V4, V5, F1 ) CALL BISECT( ST,NST, NW, V2, V4 ) C IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN *VDIR NODEP DO 15 L = 0, NW -1 IF( V4(L) .EQ. ST(V2(L)) ) THEN VEC(V5(L) ,ID) = VEC(V5(L) ,ID) & + VEC(V2(L) ,IS)*F1(L)*CS(N) END IF 15 CONTINUE ELSE *VDIR NODEP DO 20 L = 0, NW -1 IF( V4(L) .EQ. ST(V2(L)) ) THEN VEC(V5(L) ,ID) = VEC(V5(L) ,ID) & + VEC(V2(L) ,IS)*F1(L)*CS(N) & - VEC(V2(L)+NST,IS)*F1(L)*SN(N) VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID) & + VEC(V2(L) ,IS)*F1(L)*SN(N) & + VEC(V2(L)+NST,IS)*F1(L)*CS(N) END IF 20 CONTINUE END IF C --- S++S-- Part --- CALL NNPPMM( NS, IP, VLE, V0, V2, A, F0 ) CALL COMPRS( ST(NST-1), IX, VLE, NW, V2, F0, V4, V5, F1 ) CALL BISECT( ST,NST, NW, V2, V4 ) C IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN *VDIR NODEP DO 25 L = 0, NW -1 IF( V4(L) .EQ. ST(V2(L)) ) THEN VEC(V5(L) ,ID) = VEC(V5(L) ,ID) & + VEC(V2(L) ,IS)*F1(L)*CS(N) END IF 25 CONTINUE ELSE *VDIR NODEP DO 30 L = 0, NW -1 IF( V4(L) .EQ. ST(V2(L)) ) THEN VEC(V5(L) ,ID) = VEC(V5(L) ,ID) & + VEC(V2(L) ,IS)*F1(L)*CS(N) & - VEC(V2(L)+NST,IS)*F1(L)*SN(N) VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID) & + VEC(V2(L) ,IS)*F1(L)*SN(N) & + VEC(V2(L)+NST,IS)*F1(L)*CS(N) END IF 30 CONTINUE END IF C --- S-S+ Part --- CALL NNSMSP( NS, IP, VLE, V0, V2, A, F0 ) CALL COMPRS( ST(NST-1), IX, VLE, NW, V2, F0, V4, V5, F1 ) CALL BISECT( ST,NST, NW, V2, V4 ) C IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN *VDIR NODEP DO 35 L = 0, NW -1 IF( V4(L) .EQ. ST(V2(L)) ) THEN VEC(V5(L) ,ID) = VEC(V5(L) ,ID) & + VEC(V2(L) ,IS)*F1(L)*CS(N) END IF 35 CONTINUE ELSE *VDIR NODEP DO 40 L = 0, NW -1 IF( V4(L) .EQ. ST(V2(L)) ) THEN VEC(V5(L) ,ID) = VEC(V5(L) ,ID) & + VEC(V2(L) ,IS)*F1(L)*CS(N) & - VEC(V2(L)+NST,IS)*F1(L)*SN(N) VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID) & + VEC(V2(L) ,IS)*F1(L)*SN(N) & + VEC(V2(L)+NST,IS)*F1(L)*CS(N) END IF 40 CONTINUE END IF C --- S--S++ Part --- CALL NNMMPP( NS, IP, VLE, V0, V2, A, F0 ) CALL COMPRS( ST(NST-1), IX, VLE, NW, V2, F0, V4, V5, F1 ) CALL BISECT( ST,NST, NW, V2, V4 ) C IF( NK .EQ. 0 .OR. NK .EQ. NS/2 ) THEN *VDIR NODEP DO 45 L = 0, NW -1 IF( V4(L) .EQ. ST(V2(L)) ) THEN VEC(V5(L) ,ID) = VEC(V5(L) ,ID) & + VEC(V2(L) ,IS)*F1(L)*CS(N) END IF 45 CONTINUE ELSE *VDIR NODEP DO 50 L = 0, NW -1 IF( V4(L) .EQ. ST(V2(L)) ) THEN VEC(V5(L) ,ID) = VEC(V5(L) ,ID) & + VEC(V2(L) ,IS)*F1(L)*CS(N) & - VEC(V2(L)+NST,IS)*F1(L)*SN(N) VEC(V5(L)+NST,ID) = VEC(V5(L)+NST,ID) & + VEC(V2(L) ,IS)*F1(L)*SN(N) & + VEC(V2(L)+NST,IS)*F1(L)*CS(N) END IF 50 CONTINUE END IF 60 CONTINUE C --- State Translation --- IF( N .EQ. NS-1 ) GOTO 80 DO 70 J = 0, VLE -1 V0(J)= MOD(V0(J)*3,3**NS) + V0(J)/3**(NS-1) 70 CONTINUE 80 CONTINUE 90 CONTINUE CALL RESCAL( ID, IS, NS, LV, NK,NST, NX,EXS,NRM,VEC) C C --- Diagonal Part --- IF( NK.EQ.0 .OR. NK.EQ.NS/2 ) THEN DO 95 J = 0, NST -1 VEC( J,ID)= VEC( J,ID) + DIA(J)*VEC( J,IS) 95 CONTINUE ELSE DO 96 J = 0, NST -1 VEC( J,ID)= VEC( J,ID) + DIA(J)*VEC( J,IS) VEC(NST+J,ID)= VEC(NST+J,ID) + DIA(J)*VEC(NST+J,IS) 96 CONTINUE END IF C IF( LV .GE. 3 ) THEN WRITE(6,*) 'Stopped by Debug Switch LV=', LV STOP END IF C RETURN END CCC SUBROUTINE NNSPSM( NS, IP, VLE, V0, V2, A, F0 ) INTEGER NS, IP, VLE, V0(0:141072), V2(0:141072) INTEGER P0, P1, T0, T1 REAL*8 A, F0(0:141072) C C --- S+[T1] S-[T0] Part --- C IF( IP .NE. NS-1 ) THEN P0 = 3** IP P1 = 3**(IP+1) DO 10 I = 0, VLE -1 T0 = MOD(V0(I)/P0,3) T1 = MOD(V0(I)/P1,3) IF( T1.NE.2 .AND. T0.NE.0 ) THEN V2(I)= V0(I) + P1 - P0 F0(I)= 1.0D0 - A*( 1 - ABS(T0+T1-2) ) ELSE F0(I)= 0.0D0 END IF 10 CONTINUE ELSE P0 = 1 P1 = 3**(NS-1) DO 20 I = 0, VLE -1 T0 = MOD(V0(I) ,3) T1 = MOD(V0(I)/P1,3) IF( T1.NE.2 .AND. T0.NE.0 ) THEN V2(I)= V0(I) + P1 - P0 F0(I)= 1.0D0 - A*( 1 - ABS(T0+T1-2) ) ELSE F0(I)= 0.0D0 END IF 20 CONTINUE END IF C RETURN END CCC SUBROUTINE NNPPMM( NS, IP, VLE, V0, V2, A, F0 ) INTEGER NS, IP, VLE, V0(0:141072), V2(0:141072) INTEGER P0, P1, T0, T1 REAL*8 A, F0(0:141072) C C --- S++[T1] S--[T0] Part --- C IF( IP .NE. NS-1 ) THEN P0 = 3** IP P1 = 3**(IP+1) DO 10 I = 0, VLE -1 T0 = MOD(V0(I)/P0,3) T1 = MOD(V0(I)/P1,3) IF( T1.EQ.0 .AND. T0.EQ.2 ) THEN V2(I)= V0(I) + P1 + P1 - P0 - P0 F0(I)= A ELSE F0(I)= 0.0D0 END IF 10 CONTINUE ELSE P0 = 1 P1 = 3**(NS-1) DO 20 I = 0, VLE -1 T0 = MOD(V0(I) ,3) T1 = MOD(V0(I)/P1,3) IF( T1.EQ.0 .AND. T0.EQ.2 ) THEN V2(I)= V0(I) + P1 + P1 - P0 - P0 F0(I)= A ELSE F0(I)= 0.0D0 END IF 20 CONTINUE END IF C RETURN END CCC SUBROUTINE NNSMSP( NS, IP, VLE, V0, V2, A, F0 ) INTEGER NS, IP, VLE, V0(0:141072), V2(0:141072) INTEGER P0, P1, T0, T1 REAL*8 A, F0(0:141072) C C --- S-[T1] S+[T0] Part --- C IF( IP .NE. NS-1 ) THEN P0 = 3** IP P1 = 3**(IP+1) DO 10 I = 0, VLE -1 T0 = MOD(V0(I)/P0,3) T1 = MOD(V0(I)/P1,3) IF( T1.NE.0 .AND. T0.NE.2 ) THEN V2(I)= V0(I) - P1 + P0 F0(I)= 1.0D0 - A*( 1 - ABS(T0+T1-2) ) ELSE F0(I)= 0.0D0 END IF 10 CONTINUE ELSE P0 = 1 P1 = 3**(NS-1) DO 20 I = 0, VLE -1 T0 = MOD(V0(I) ,3) T1 = MOD(V0(I)/P1,3) IF( T1.NE.0 .AND. T0.NE.2 ) THEN V2(I)= V0(I) - P1 + P0 F0(I)= 1.0D0 - A*( 1 - ABS(T0+T1-2) ) ELSE F0(I)= 0.0D0 END IF 20 CONTINUE END IF C RETURN END CCC SUBROUTINE NNMMPP( NS, IP, VLE, V0, V2, A, F0 ) INTEGER NS, IP, VLE, V0(0:141072), V2(0:141072) INTEGER P0, P1, T0, T1 REAL*8 A, F0(0:141072) C C --- S--[T1] S++[T0] Part --- C IF( IP .NE. NS-1 ) THEN P0 = 3** IP P1 = 3**(IP+1) DO 10 I = 0, VLE -1 T0 = MOD(V0(I)/P0,3) T1 = MOD(V0(I)/P1,3) IF( T1.EQ.2 .AND. T0.EQ.0 ) THEN V2(I)= V0(I) - P1 - P1 + P0 + P0 F0(I)= A ELSE F0(I)= 0.0D0 END IF 10 CONTINUE ELSE P0 = 1 P1 = 3**(NS-1) DO 20 I = 0, VLE -1 T0 = MOD(V0(I) ,3) T1 = MOD(V0(I)/P1,3) IF( T1.EQ.2 .AND. T0.EQ.0 ) THEN V2(I)= V0(I) - P1 - P1 + P0 + P0 F0(I)= A ELSE F0(I)= 0.0D0 END IF 20 CONTINUE END IF C RETURN END CCC SUBROUTINE COMPRS( MST, IX, VLE, NW, V2, F0, V4, V5, F1 ) INTEGER MST, IX, VLE, NW, V2(0:141072) INTEGER V4(0:141072), V5(0:141072) REAL*8 F0(0:141072), F1(0:141072) C NW = 0 DO 10 I = 0, VLE -1 IF( ABS(F0(I)) .GT. 1.0E-16 .AND. V2(I) .LE. MST ) THEN V4(NW) = V2(I) F1(NW) = F0(I) V5(NW) = IX + I NW = NW + 1 END IF 10 CONTINUE C RETURN END CCC SUBROUTINE BISECT( ST, NST, NW, V2, V4 ) INTEGER DW, ST(0:122222), NST, NW, V2(0:141072), V4(0:141072) C DW = NST - NST/2 DO 10 I = 0, NW -1 IF( V4(I) .GE. ST(DW) ) THEN V2(I) = DW ELSE V2(I) = 0 END IF 10 CONTINUE C 20 CONTINUE DW = DW - DW/2 DO 30 I = 0, NW -1 IF( V4(I) .GE. ST( MIN(NST-1,V2(I)+DW) ) ) THEN V2(I) = MIN(NST-1,V2(I)+DW) END IF 30 CONTINUE IF( DW .NE. 1 ) GOTO 20 C RETURN END CCC CC C SUBROUTINE LANZOS(NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC ) INTEGER LK,LLANCZ,NS,LV,NK,NX,NST,NV,EXS(0:999), ST(0:122222) REAL*8 E2, E2P, ENERGY, VEC(0:244444,0:3) REAL*8 SN(0:40),CS(0:40), NRM(0:999),A,DIA(0:122222) REAL*8 MER(1000,5),TER(1000), ALP(1000), BET(1000) C E2 = 0.0D0 C IF(NK.EQ.0 .OR.NK.EQ.NS/2) THEN NV = NST ELSE NV = 2*NST END IF C CALL VEC0( NV, VEC) CALL VECINI( NS, NK, NST, NX, EXS, NRM, VEC) CALL VMOVE( NV, 3, 0, VEC) C (I) CALL MULT( 1,0, NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC) C (II) C ******** LK is an Index for ALP and BET ******** LK = 1 LLANCZ = 0 10 CONTINUE C ******** Report every 5 Steps ******** E2P = E2 DO 20 J = 0, 4 CALL MACRO1( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC, & TER,ALP,BET, LK, 0 ) 20 CONTINUE CALL MALANZ( LK, LV, ALP,BET, ENERGY,E2, MER) LLANCZ = LLANCZ + 1 IF( LLANCZ.LE.40 .AND. ABS(E2P-E2).GE.1.0E-16) GOTO 10 LLANCZ = LLANCZ - 1 WRITE(6,100)ENERGY 100 FORMAT(1H , 'ENERGY= ',E30.20) C IF( ABS(ENERGY) .LT. 1.0E-6 ) THEN WRITE(6,*) 'Energy is accidentally Zero' WRITE(6,*) 'Do not believe the following data' RETURN END IF C DO 30 I = 1,1000 TER(I)= MER(I,1) 30 CONTINUE C CALL VMOVE( NV, 0, 3, VEC) CALL ADDER( NV, 3, 3, 0,TER(1) , 0.0D0, VEC) CALL MULT( 1,0, NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC) LK=1 DO 70 I = 0, LLANCZ DO 60 J = 0, 4 CALL MACRO1( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC, & TER,ALP,BET, LK, 1 ) 60 CONTINUE 70 CONTINUE CALL INVERS( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,ENERGY,VEC) WRITE(6,110) ENERGY 110 FORMAT(1H ,'Energy = ', F20.15) C RETURN END CCC SUBROUTINE MACRO1( NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC, & TER,ALP,BET, LK,ISW ) INTEGER ISW, LK,NS,LV,NK,NX,NST,NV,EXS(0:999) INTEGER ST(0:122222) REAL*8 PROD, VEC(0:244444,0:3) REAL*8 SN(0:40), CS(0:40), NRM(0:999),A,DIA(0:122222) REAL*8 TER(1000), ALP(1000), BET(1000) C IF(NK.EQ.0 .OR.NK.EQ.NS/2) THEN NV = NST ELSE NV = 2*NST END IF C CALL INPRO( NV, 0, 1, PROD,VEC) ALP(LK) = PROD PROD = -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( 1, 2,NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC) PROD = -PROD CALL ADDER( NV, 1, 1, 0,1.0D0,PROD,VEC) CALL VMOVE( NV, 0, 2, VEC) LK = LK + 1 IF(ISW.EQ.1) CALL ADDER( NV, 3, 3, 0, 1.0D0, TER(LK), VEC) C RETURN END CCC SUBROUTINE MALANZ( LK, LV, ALP, BET, ENERGY, E2, VE) INTEGER N,M, LNV, LK, LV, ISW,IERR, IW1(1000) REAL*8 ENERGY, E2,ALP(1000), BET(1000) REAL*8 TP, EPS, E(1000), W1(6000), VE(1000,5) C N = LK -1 EPS =-1.0 M = 5 LNV = 1000 ISW =-1 TP = BET(N) CALL BCSTSS(ALP,N,BET,EPS,E,M,VE,LNV,ISW,IW1,W1,IERR) BET(N)= TP IF(IERR.NE.0) THEN WRITE(6,100) IERR 100 FORMAT(1H ,' Error No. ',I10,' in RCSTSS') END IF ENERGY= E(1) E2 = E(2) IF( LV .GE. 2 ) WRITE(6,110) (E(J),J=1,3) 110 FORMAT(1H ,'*****',3F20.15,'*****') C RETURN END CCC CC C SUBROUTINE INVERS(NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,LMD,VEC) INTEGER LROT,NS,LV,NK,NX,NST,NV, EXS(0:999) INTEGER ST(0:122222) REAL*8 LMD,A,DIA(0:122222),VEC(0:244444,0:3) REAL*8 SN(0:40),CS(0:40),DTP, DTR, NRM(0:999) C IF(NK.EQ.0 .OR.NK.EQ.NS/2) THEN NV = NST ELSE NV = 2*NST END IF C CALL VMOVE( NV, 0, 3, VEC) LROT = 0 GOTO 20 C 10 CONTINUE CALL CGMETH(NS,LV,NK, NX,NST,ST,EXS,SN,CS,NRM,A,DIA,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( 2,1, NS, LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC) CALL INPRO( NV, 2, 2, DTP, VEC) DTP = SQRT( DTP ) CALL DIVER( NV, 2, 2, DTP, VEC) C 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 C CALL INPRO( NV, 3, 3, DTR, VEC) WRITE(6,100) LMD, DTP, SQRT(DTR) LROT = LROT + 1 100 FORMAT(1H , 'LMD=',E20.10, ' ? ', E20.10, ' D ',E20.10) IF( LROT .GE. 20 ) THEN WRITE(6,110) 110 FORMAT(1H ,'Did not converge in INVERS') STOP END IF IF( SQRT(DTR) .GE. 0.5E-11 ) GOTO 10 C RETURN END CCC SUBROUTINE CGMETH(NS,LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,LMD,VEC) INTEGER NV, LROT,NS,LV,NK,NX,NST,EXS(0:999), ST(0:122222) REAL*8 LMD,DB,A,ALP,BET, VEC(0:244444,0:3),DIA(0:122222) REAL*8 EPS, DTP,SN(0:40), CS(0:40), DR1,DR2, NRM(0:999) C IF(NK.EQ.0 .OR.NK.EQ.NS/2) THEN NV = NST ELSE NV = 2*NST END IF 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) EPS = 1.0E-8 LROT = 0 CALL INPRO( NV, 1, 1, DR1, VEC) 20 CALL MULT( 3,2, NS, LV,NK,NX,NST,ST,EXS,SN,CS,NRM,A,DIA,VEC) 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 C-WR- WRITE(6,100) DR2 / DB C 100 FORMAT(1H ,E30.20) IF( DR2/DB .LE. EPS ) GOTO 40 C BET = DR2 / DR1 DR1 = DR2 LROT = LROT + 1 CALL ADDER( NV, 2, 1, 2, 1.0D0, BET,VEC) IF( LROT .GE. 20 ) THEN WRITE(6,110) 110 FORMAT(1H ,'Did not converge in CG') GOTO 50 END IF GOTO 20 40 CONTINUE 50 WRITE(6,120) LROT 120 FORMAT(1H ,I10,'After Rotation in CG') C RETURN END CCC CC C SUBROUTINE VECDMP( NS, LV, NK, NST, ST, VEC ) INTEGER IP, NS, LV, NK, NST, ST(0:122222) REAL*8 FTP, VEC(0:244444,0:3) C IF( LV .EQ. 0 ) RETURN C IF( NST .GT. 500 ) THEN WRITE(6,*) ' Too Many Basis ' RETURN END IF C DO 10 I = 0, NST -1 IF( NK.EQ.0 .OR. NK.EQ.NS/2 ) THEN VEC(I,3) = VEC(I,1)*VEC(I,1) ELSE VEC(I,3) = VEC(I,1)*VEC(I,1) + VEC(NST+I,1)*VEC(NST+I,1) END IF 10 CONTINUE C DO 30 I = 0, MIN(25,NST) -1 FTP = 0.0D0 IP = 0 DO 20 J = 0, NST -1 IF( FTP .LT. VEC(J,3) ) THEN FTP = VEC(J,3) IP = J END IF 20 CONTINUE IF( FTP .GT. 1.0E-8 ) CALL STDISP( NS, ST(IP), VEC(IP,3) ) VEC(IP,3) = 0.0D0 30 CONTINUE C RETURN END CCC CC C SUBROUTINE SZDIST( NS, NK, NST, ST, BH, BL, VEC ) INTEGER TP, M, NS, NK, NST, ST(0:122222) INTEGER TH,TL,BH(0:59049,0:2), BL(0:59049,0:2) REAL*8 FC, F2, F1, F0,VEC(0:244444,0:3) C M = NS /2 TP = 3** M F0 = 0.0D0 F1 = 0.0D0 F2 = 0.0D0 C IF( NK.EQ.0 .OR. NK.EQ.NS/2 ) THEN DO 10 I = 0, NST -1 TH = ST(I)/TP TL = MOD(ST(I),TP) FC = VEC(I,1) * VEC(I,1) F0 = F0 + FC*( BH(TH,0) + BL(TL,0) ) F1 = F1 + FC*( BH(TH,1) + BL(TL,1) ) F2 = F2 + FC*( BH(TH,2) + BL(TL,2) ) 10 CONTINUE ELSE DO 20 I = 0, NST -1 TH = ST(I)/TP TL = MOD(ST(I),TP) FC = VEC(I,1) *VEC(I,1) + VEC(NST+I,1) *VEC(NST+I,1) F0 = F0 + FC*( BH(TH,0) + BL(TL,0) ) F1 = F1 + FC*( BH(TH,1) + BL(TL,1) ) F2 = F2 + FC*( BH(TH,2) + BL(TL,2) ) 20 CONTINUE END IF C WRITE(6,100) F0/NS, F1/NS, F2/NS 100 FORMAT(1H ,'v 0 ^ Ratio = ', 3F15.7) C RETURN END CCC CC C SUBROUTINE BCSTSS(A,N,B,EPS,E,M,VE,LNV,ISW,WN,W1,IERR) INTEGER MWN,NC,N,M,SP,BP,ISW,LNV,IERR,WN(0:1000) REAL*8 EPS,A(0:1000),B(0:1000),W1(0:6000),E(0:1000) REAL*8 MIN,MDL,BB(0:1000),VE(0:1000,0:2) C -DUMMY- IERR = ISW IERR = LNV IERR = 0 BP = 0 MIN = ABS(A(0)) + ABS(B(0)) DO 10 I = 0,N -3 MIN = MAX(MIN,ABS(B(I))+ABS(A(I+1))+ABS(B(I+1))) BB(I)= B(I)*B(I) 10 CONTINUE MIN = MAX(MIN,ABS(B(N-2))+ABS(A(N-1))) BB(N-2)= B(N-2)*B(N-2) C CALL COUNT(N,NC,-MIN,A,BB) MIN = -MIN MWN = NC CALL COUNT(N,NC,-MIN,A,BB) SP = 0 W1(SP)= -MIN WN(SP)= NC C 20 CONTINUE MDL =(MIN+W1(SP)) / 2 IF(W1(SP)-MIN.LT.EPS.OR.W1(SP)-MDL.LT.ABS(MDL)*1.0E-17.OR.MDL & -MIN.LT.ABS(MDL)*1.0E-17)GOTO 25 C CALL COUNT(N,NC,MDL,A,BB) C -WR- WRITE(6,198) SP,MWN-WN(SP), W1(SP)-MDL, MDL-MIN C 198 FORMAT(1H , 'SP,MWN-WN(SP), W1(SP)-MDL, MDL-MIN', 2I3, 2E20.10 ) IF(MWN-NC.EQ.0)THEN MIN = MDL MWN = NC GOTO 20 END IF C SP = SP + 1 W1(SP)= MDL WN(SP)= NC GOTO 20 C 25 CONTINUE C DO 30 I = 0,MWN - WN(SP) - 1 E(BP+I)= MDL 30 CONTINUE BP = BP + MWN - WN(SP) C IF(BP.LT.M) THEN 40 CONTINUE IF(WN(SP).EQ.WN(SP-1)) THEN SP = SP-1 GOTO 40 END IF MIN = W1(SP) MWN = WN(SP) SP = SP - 1 GOTO 20 END IF C CALL INVERT(N,E(0),A,B,W1,VE) C RETURN END CCC CC C SUBROUTINE COUNT(N,NC,LAM,A,BB) INTEGER I,N,NC REAL*8 G(0:1000),GRN,EPS,LAM,A(0:1000),BB(0:1000) C NC = 0 EPS = 1.0E-18 GRN = 1.0E+18 G(0) = LAM - A(0) DO 10 I = 1,N- 1 IF(ABS(G(I-1)).GT.EPS) THEN G(I) = LAM - A(I) - BB(I-1)/G(I-1) ELSE G(I) = GRN END IF 10 CONTINUE C DO 20 I = 0,N-1 IF(G(I).LE.0) NC = NC + 1 20 CONTINUE C RETURN END CCC CC C SUBROUTINE INVERT(N,E,A,B,W1,VE) INTEGER N REAL*8 PROD,MROD REAL*8 E,A(0:1000),B(0:1000),W1(0:1000,0:5),VE(0:1000,0:2) C VE(0,0)= 1.0D0 DO 10 I = 1,N-1 VE(I,0)= 0.0D0 10 CONTINUE C CALL CHRSKY(N,E,A,B,W1) C 20 CONTINUE DO 30 I = 0,N-1 VE(I,1)= VE(I,0) 30 CONTINUE CALL VSOLVE(N,W1,VE) C PROD = 0.0D0 DO 40 I = 0,N-1 PROD = PROD + VE(I,0)*VE(I,0) 40 CONTINUE C PROD = SQRT(PROD) DO 50 I = 0,N-1 VE(I,0)= VE(I,0) / PROD 50 CONTINUE C PROD = 0.0D0 MROD = 0.0D0 DO 60 I = 0,N-1 PROD = PROD +(VE(I,1)-VE(I,0))*(VE(I,1)-VE(I,0)) MROD = MROD +(VE(I,1)+VE(I,0))*(VE(I,1)+VE(I,0)) 60 CONTINUE C -WR- WRITE(6, * )' MROD, PROD ', MROD, PROD IF(PROD.GT.1.0E-16.AND.MROD.GT.1.0E-16) GOTO 20 C RETURN END CCC CC C SUBROUTINE CHRSKY(N,E,A,B,W1) INTEGER N REAL*8 E,A(0:1000),B(0:1000),W1(0:1000,0:5) C W1(0,0)= A(0)- E DO 10 I = 0,N - 2 W1(I,1)= B(I)/W1(I,0) W1(I+1,0)= A(I+1)- E - B(I)*W1(I,1) 10 CONTINUE C RETURN END CCC CC C SUBROUTINE VSOLVE(N,W1,VE) INTEGER N REAL*8 VE(0:1000,0:2),W1(0:1000,0:5) C DO 20 I = 1,N-1 VE(I,0)= VE(I,0) - W1(I-1,1)*VE(I-1,0) 20 CONTINUE C DO 30 I = 0,N-1 VE(I,0)= VE(I,0)/W1(I,0) 30 CONTINUE C DO 40 I = N -2,0,-1 VE(I,0)= VE(I,0) - W1(I,1)*VE(I+1,0) 40 CONTINUE C RETURN END