DMRG for 3D Ising model (Old version)
This program runs with NEC's ASL library
(Please contact me if you wish to know about ASL.)
C123456789012345678901234567890123456789012345678901234567890123456789012
C ----- -----
C ----- Corner Tensor Renormalization Group Method -----
C ----- for the Cubic Lattice Ising Model Tc -----
C ----- -----
C ----- '97/ 9/24, Ver.2.0 by T.Nishino -----
C ----- Conditions:
C ----- q = 2 Ising
C ----- M1 = m = 2 Line Block Spin
C ----- M2 = M = 2 Surface Block Spin
C ----- Tensor Size:
C ----- Type A: q^6 = 2^6 = 64
C ----- Type B: q x (qm)^4 = 2^9 = 512
C ----- Type C: (qm)^2 x (qmmM)^2 = 2^12 = 4096
C ----- Type D: (qmmM)^3 = 2^12 = 4096
C ----- Type E: ((qmmM)^2)^2 = 2^16 = 65536
C ----- Type M: (qmmM)^2 = 2^8 = 256
C
C PARAMETER( MA=64, MB=512, MC=4096, MD=4096, ME=65536, MM=256 )
C
PARAMETER( IQ=2, MS=3, ML=3 )
PARAMETER( MA=IQ**6 )
PARAMETER( MB=IQ*(IQ*MS)**4 )
PARAMETER( MC=(IQ*MS)**2*(IQ*MS*MS*ML)**2 )
PARAMETER( MD=(IQ*MS*MS*ML)**3 )
PARAMETER( ME=(IQ*MS*MS*ML)**4 )
PARAMETER( MM=(IQ*MS*MS*ML)**2 )
C
INTEGER M1, M1X, M1N, L
INTEGER M2, M2X, M2N, ISW
REAL*8 EA, EB, EC, ED, KP1, Y
REAL*8 T, FB, FC, FD, KP2, Z, F
REAL*8 A1(0:MA), A2(0:MA)
REAL*8 B1(0:MB), B2(0:MB), B3(0:MB)
REAL*8 C1(0:MC), C2(0:MC), C3(0:MC)
REAL*8 D1(0:MD), D2(0:MD), D3(0:MD)
REAL*8 E1(0:ME), E2(0:ME), E3(0:ME), E4(0:ME)
REAL*8 DM1(MM), DM2(MM), W1(0:MM), W2(0:MM), PM(0:MM)
C
C ***** Initial Steps *****
CALL GETPRM( L,M1X,M2X, T )
C
C Write(6,*) ' 0:Normal 1:T=0 2:T=inf.'
C Read(5,*) ISW
CALL INITTT( A1, A2, B1, C1, D1, T )
C IF( ISW .EQ. 1 ) CALL INITTZ( A1, A2, B1, C1, D1 )
C IF( ISW .EQ. 2 ) CALL INITTI( A1, A2, B1, C1, D1 )
C
M1 = 2
M2 = 2
K = 1
PM(0)= 1.0D0
PM(1)= 1.0D0
CALL NORM( 2**6, A1, EA )
CALL NRA2( 2**6, A2, EA )
EA = LOG(EA)
EB = 0.0D0
EC = 0.0D0
ED = 0.0D0
Write(6,*) 'KP2, Kp1, F, F**(1/3), F/N, MG'
10 CONTINUE
C
C Write(6,*) '--- Step Number ', K,' ---
C ***** Stack Steps *****
CALL PROC1( M1,M2, C1, D1,D2 )
CALL PROC2( M1,M2, B1, C1,C2 )
CALL PROC3( M1,M2, C2, D2,D1 )
CALL PROC4( M1, A1,B1,B2 )
CALL PROC4( M1, A2,B1,B3 )
CALL PROC5( M1,M2, B2,C2,C3 )
CALL PROC5( M1,M2, B3,C2,C1 )
CALL PROC6( M1,M2, C3,D1,D2 )
CALL PROC6( M1,M2, C1,D1,D3 )
CALL PROC7( M1,M2, D2,E1 )
CALL PROC8( M1,M2, E1,E2 )
CALL PROC8( M1,M2, E2,E3 )
C
C ***** Trace Steps *****
CALL CRDM2( M1,M2, E3,DM2 )
CALL SYMCDM( M1,M2, PM,DM2 )
CALL CRDM1( M1,M2, DM2,DM1 )
CALL CREZ( M1,M2, Z, E3 )
C
C ***** Observation Steps *****
CALL STEP7( M1,M2, D2,D3,E4 )
CALL STEP8( M1,M2, E1,E4,E3 )
CALL STEP8( M1,M2, E2,E3,E4 )
CALL CREZ( M1,M2, Y, E4 )
C
C ***** Cut-off Steps *****
CALL DIAG( 2*M1*M1*M2, DM2, E1, W2 )
M2N = MIN( 2*M1*M1*M2, M2X )
CALL KPTWT( 2*M1*M1*M2, M2N, KP2, W2 )
CALL PARITY( M1, M2, DM2, PM, W2 )
CALL DIAG( 2*M1, DM1, E1, W1 )
M1N = MIN( 2*M1, M1X )
CALL KPTWT( 2*M1, M1N, KP1, W1 )
C
C ***** Display Steps *****
F = ( EA*8 +EB*24 +EC*24 +ED*8 +LOG(Z) )/LOG(2.0D0)
Write(6,100) KP2, KP1, F/(K*2+2)**3, Y/Z
100 FORMAT(1H ,'KW:',2F11.7,' f= ',F20.15,' m= ',F20.15 )
C
C ***** Renormalization Steps *****
CALL RENOB( M1,M1N, DM1, B2, B1 )
CALL TRPOS( M1, M2, PM, DM2, W2 )
CALL RENOC( M1,M1N,M2,M2N, DM1, DM2, W2, C3, C1 )
CALL RENOD( M1, M2,M2N, DM2, W2, D2, D1 )
M1 = M1N
M2 = M2N
C
C ***** Symmetry Check Steps *****
C Write(6,*) ' - After Renorm -'
CALL SYMCHA( A1 )
CALL SYMCHB( M1, B1 )
CALL SYMCHC( M1,M2, C1 )
CALL SYMCHD( M2, D1 )
C
20 CONTINUE
C ***** Normalization Steps *****
CALL NORM( 2*M1*M1*M1*M1, B1, FB )
CALL NORM( M1*M1*M2*M2, C1, FC )
CALL NORM( M2*M2*M2, D1, FD )
ED = ED + LOG(FD) + EA + 3*EB + 3*EC
EC = EC + LOG(FC) + EA + 2*EB
EB = EB + LOG(FB) + EA
K = K + 1
IF( K .LT. L+1 ) GOTO 10
STOP
END
CCC
CC
C
SUBROUTINE NZDUMP( N, A )
INTEGER N
REAL*8 A(0:*)
Write(6,*) 'dim.', N
DO 10 I = 0, N-1
IF( ABS(A(I)) .GT. 1.0E-8 ) Write(6,*) I, A(I)
10 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE GETPRM( L, M1X, M2X, T )
INTEGER L, M1X, M2X
REAL*8 T
C
Write(6,*) 'Input Max. States m1x'
Read(5,*) M1X
C M1X = 2
Write(6,*) 'Input Max. States m2x'
Read(5,*) M2X
C M2X = 3
Write(6,*) 'Input Max. Iteration # L'
Read(5,*) L
C L = 100
Write(6,*) 'Input Min. Temperature T'
Read(5,*) T
C T = 2.0
Write(6,*) ' m1x = ',M1X,' m2x = ',M2X
Write(6,*) ' L = ', L,' T = ', T
RETURN
END
CCC
CC
C
SUBROUTINE INITTT( A1, A2, B, C, D, T )
INTEGER P1, P2, P3, P4
REAL*8 T, XI, A1(0:*), A2(0:*), B(0:*), C(0:*), D(0:*)
C
C A1: Vertex Weight
C i j i j i j i j A2: Magnetization
C |/ |/ |/ |/
C D---k C---k m---B---k m---A---k B,C,D
C / / /| --- Fixed B.C.
C L L L n
C
C note: All indices are equivalent;
C order of indices (ijklmn) is irrelevant.
C
XI = EXP(2.0D0/T) + SQRT( EXP(2.0D0/T)**2 - 1.0D0 )
DO 10 I = 0, 1
DO 10 J = 0, 1
DO 10 K = 0, 1
P1 = I*4 + J*2 + K
D(P1)= XI**( I+J+K -1.5D0) / 2.0D0**1.5D0
C & + XI**(-I-J-K +1.5D0) / 2.0D0**1.5D0
DO 10 L = 0, 1
P2 = P1*2 + L
C(P2)= XI**( I+J+K+L -2.0D0) / 2.0D0**2.0D0
C & + XI**(-I-J-K-L +2.0D0) / 2.0D0**2.0D0
DO 10 M = 0, 1
P3 = P2*2 + M
B(P3)= XI**( I+J+K+L+M -2.5D0) / 2.0D0**2.5D0
C & + XI**(-I-J-K-L-M +2.5D0) / 2.0D0**2.5D0
DO 10 N = 0, 1
P4 = P3*2 + N
A1(P4)= XI**( I+J+K+L+M+N-3.0D0) / 2.0D0**3.0D0
& + XI**(-I-J-K-L-M-N+3.0D0) / 2.0D0**3.0D0
A2(P4)= XI**( I+J+K+L+M+N-3.0D0) / 2.0D0**3.0D0
& - XI**(-I-J-K-L-M-N+3.0D0) / 2.0D0**3.0D0
10 CONTINUE
RETURN
END
CCC
SUBROUTINE INITTZ( A1, A2, B, C, D )
INTEGER P1, P2, P3, P4
REAL*8 A1(0:*), A2(0:*), B(0:*), C(0:*), D(0:*)
C
C A1: Vertex Weight
C i j i j i j i j A2: Magnetization
C |/ |/ |/ |/
C D---k C---k m---B---k m---A---k B,C,D
C / / /| --- Fixed B.C.
C L L L n
C
C note: Zero Temperature
C
Write(6,*) '--- Debug Mode T = 0'
C
DO 10 I = 0, 1
DO 10 J = 0, 1
DO 10 K = 0, 1
P1 = I*4 + J*2 + K
IF( I+J+K .EQ. 3 ) THEN
D(P1)= 1.0D0
ELSE
D(P1)= 0.0D0
END IF
DO 10 L = 0, 1
P2 = P1*2 + L
IF( I+J+K+L .EQ. 4 ) THEN
C(P2)= 1.0D0
ELSE
C(P2)= 0.0D0
END IF
DO 10 M = 0, 1
P3 = P2*2 + M
IF( I+J+K+L+M .EQ. 5 ) THEN
B(P3)= 1.0D0
ELSE
B(P3)= 0.0D0
END IF
DO 10 N = 0, 1
P4 = P3*2 + N
A1(P4)= 0.0D0
A2(P4)= 0.0D0
IF( I+J+K+L+M+N .EQ. 6 ) THEN
A1(P4)= 1.0D0
A2(P4)= 1.0D0
END IF
IF( I+J+K+L+M+N .EQ. 0 ) THEN
A1(P4)= 1.0D0
A2(P4)=-1.0D0
END IF
10 CONTINUE
RETURN
END
CCC
SUBROUTINE INITTI( A1, A2, B, C, D )
INTEGER P1, P2, P3, P4
REAL*8 A1(0:*), A2(0:*), B(0:*), C(0:*), D(0:*)
C
C A1: Vertex Weight
C i j i j i j i j A2: Magnetization
C |/ |/ |/ |/
C D---k C---k m---B---k m---A---k B,C,D
C / / /| --- Free B.C.
C L L L n
C
C note: Infinite Temperature
C
Write(6,*) '--- Debug Mode T = inf.'
C
DO 10 I = 0, 1
DO 10 J = 0, 1
DO 10 K = 0, 1
P1 = I*4 + J*2 + K
D(P1)= 2.0D0 / 2.0D0**1.5D0
DO 10 L = 0, 1
P2 = P1*2 + L
C(P2)= 2.0D0 / 2.0D0**2.0D0
DO 10 M = 0, 1
P3 = P2*2 + M
B(P3)= 2.0D0 / 2.0D0**2.5D0
DO 10 N = 0, 1
P4 = P3*2 + N
A1(P4)= 2.0D0 / 2.0D0**3.0D0
A2(P4)= 0.0D0
10 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE PROC1( M1, M2, C1, D1, D2 )
INTEGER M1, M2, AC, AD1, AD2
REAL*8 C1(0:*), D1(0:*), D2(0:*)
C [n]
C |
C D--[m]
C /
C j[L]
C |/
C C---i
C /
C [k]
C
CALL MATCLR( M2 * M1*M2 * M1*M2, D2 )
DO 10 I = 0, M1 -1
DO 10 J = 0, M1 -1
DO 10 K = 0, M2-1
DO 10 L = 0, M2-1
AC = ((I*M1+J)*M2+K)*M2+L
DO 10 M = 0, M2-1
DO 10 N = 0, M2-1
AD1 = (L*M2+M)*M2+N
AD2 =(((K*M1+I)*M2+M)*M1+J)*M2+N
D2(AD2)= D2(AD2) + C1(AC)*D1(AD1)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE PROC2( M1, M2, B, C1, C2 )
INTEGER S,T, M1, M2, AB, AC1, AC2
REAL*8 B(0:*), C1(0:*), C2(0:*)
C n
C |
C [s]--C--[t]
C /
C (i)m
C |/
C j---B---L
C /
C k
C
CALL MATCLR( M1 * 2*M1 * M1*M2 * M1*M2, C2 )
DO 10 I = 0, 2-1
DO 10 J = 0, M1-1
DO 10 K = 0, M1-1
DO 10 L = 0, M1-1
DO 10 M = 0, M1-1
AB = (((I*M1+J)*M1+K)*M1+L)*M1+M
DO 10 N = 0, M1-1
DO 10 S = 0, M2-1
DO 10 T = 0, M2-1
AC1 = ((M*M1+N)*M2+S)*M2+T
AC2 =(((((K*2+I)*M1+N)*M1+J)*M2+S)*M1+L)*M2+T
C2(AC2)= C2(AC2) + B(AB)*C1(AC1)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE PROC3( M1, M2, C2, D2, D1 )
INTEGER M1, M2, AC, AD2, AD1
REAL*8 C2(0:*), D2(0:*), D1(0:*)
C [n] j
C | | j:2*M1
C D--[L]--C--[m] L,m,n: M1*M2
C / /
C [k] i
C
CALL MATCLR( M1*M2 * M1*M2 * 2*M1*M1*M2, D1 )
DO 10 I = 0, M1 -1
DO 10 J = 0, 2*M1 -1
DO 10 K = 0, M2-1
DO 10 L = 0, M1*M2-1
DO 10 M = 0, M1*M2-1
AC = ((I*2*M1+J)*M1*M2+L)*M1*M2+M
DO 10 N = 0, M1*M2-1
AD2 = (K*M1*M2+L)*M1*M2+N
AD1 =(((I*M2+K)*M1*M2+M)*2*M1+J)*M1*M2+N
D1(AD1)= D1(AD1) + C2(AC)*D2(AD2)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE PROC4( M1, A, B1, B2 )
INTEGER Q,R,S,T, M1, AA, AB1, AB2
REAL*8 A(0:*), B1(0:*), B2(0:*)
C
C (i L)
C |/
C (m)--A--(k)
C /|
C (j n)s
C |/
C t---B---r
C /
C q
C
CALL MATCLR( 2 * 2*M1 * 2*M1 * 2*M1 * 2*M1, B2 )
DO 10 I = 0, 2 -1
DO 10 J = 0, 2 -1
DO 10 K = 0, 2 -1
DO 10 L = 0, 2 -1
DO 10 M = 0, 2 -1
DO 10 N = 0, 2 -1
AA = ((((I*2+J)*2+K)*2+L)*2+M)*2+N
DO 10 Q = 0, M1-1
DO 10 R = 0, M1-1
DO 10 S = 0, M1-1
DO 10 T = 0, M1-1
AB1 = (((N*M1+Q)*M1+R)*M1+S)*M1+T
AB2 =(((((((I*2+J)*M1+Q)*2+K)*M1+R)*2+L)*M1+S)*2+M)*M1+T
B2(AB2)= B2(AB2) + A(AA)*B1(AB1)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE PROC5( M1, M2, B2, C2, C3 )
INTEGER S, T, M1, M2, AB, AC2, AC3
REAL*8 B2(0:*), C2(0:*), C3(0:*)
C
C k Using Mirror Symmetry of B
C |
C j---B---L j,k,L,n:2*M1
C /| s,t: M1*M2
C (i)n
C |
C [s]--C--[t]
C /
C m
C
CALL MATCLR( 2*M1 * 2*M1 * 2*M1*M1*M2 * 2*M1*M1*M2, C3 )
DO 10 I = 0, 2 -1
DO 10 J = 0, 2*M1 -1
DO 10 K = 0, 2*M1 -1
DO 10 L = 0, 2*M1 -1
DO 10 M = 0, M1 -1
DO 10 N = 0, 2*M1 -1
AB = (((I*2*M1+J)*2*M1+K)*2*M1+L)*2*M1+N
DO 10 S = 0, M1*M2-1
DO 10 T = 0, M1*M2-1
AC2 = ((M*2*M1+N)*M1*M2+S)*M1*M2+T
AC3 =(((((I*M1+M)*2*M1+K)*2*M1+J)*M1*M2+S)*2*M1+L)*M1*M2+T
C3(AC3)= C3(AC3) + B2(AB)*C2(AC2)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE PROC6( M1, M2, C3, D1, D2 )
INTEGER M1, M2, AC, AD1, AD2
REAL*8 C3(0:*), D1(0:*), D2(0:*)
C [k]
C |
C C---j
C /|
C i[n] i,j:2*M1
C | L,m: M1*M2
C D--[m] k,n:2*M1*M1*M2 (oriented)
C /
C [L]
C
CALL MATCLR( 2*M1*M1*M2 * 2*M1*M1*M2 * 2*M1*M1*M2, D2 )
DO 10 I = 0, 2*M1 -1
DO 10 J = 0, 2*M1 -1
DO 10 K = 0, 2*M1*M1*M2-1
DO 10 L = 0, M1*M2-1
DO 10 M = 0, M1*M2-1
DO 10 N = 0, 2*M1*M1*M2-1
AD1 = (L*M1*M2+M)*2*M1*M1*M2+N
AC = ((I*2*M1+J)*2*M1*M1*M2+K)*2*M1*M1*M2+N
AD2 =(((I*M1*M2+L)*2*M1+J)*M1*M2+M)*2*M1*M1*M2+K
D2(AD2)= D2(AD2) + C3(AC)*D1(AD1)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE PROC7( M1, M2, D, E )
INTEGER MF, M1, M2, AD, DD, AE
REAL*8 D(0:*), E(0:*)
C
C [M] [k]
C | |
C *D*-[i]--D Indices are oriented
C / /
C [L] [j] *D* = Mirror of D
C
MF = 2*M1*M1*M2
CALL MATCLR( MF**4, E )
DO 10 I = 0, MF-1
DO 10 J = 0, MF-1
DO 10 K = 0, MF-1
AD = (I*MF+J)*MF+K
DO 10 L = 0, MF-1
DO 10 M = 0, MF-1
DD = (I*MF+L)*MF+M
AE =((J*MF+L)*MF+K)*MF+M
E(AE)= E(AE) + D(AD)*D(DD)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE STEP7( M1, M2, D2, D3, E4 )
INTEGER MF, M1, M2, AD2, AD3, AE
REAL*8 D2(0:*), D3(0:*), E4(0:*)
C
C [M] [k]
C | |
C D3-[i]--D2 Indices are oriented
C / /
C [L] [j]
C
MF = 2*M1*M1*M2
CALL MATCLR( MF**4, E4 )
DO 10 I = 0, MF-1
DO 10 J = 0, MF-1
DO 10 K = 0, MF-1
AD2 = (I*MF+J)*MF+K
DO 10 L = 0, MF-1
DO 10 M = 0, MF-1
AD3 = (I*MF+L)*MF+M
AE =((J*MF+L)*MF+K)*MF+M
E4(AE)= E4(AE) + D2(AD2)*D3(AD3)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE PROC8( M1, M2, EA, EB )
INTEGER MFF, M1, M2
REAL*8 EA(0:*), EB(0:*)
C EA *
C /| /|
C (jXi) EA|
C | / (i)
C EA * |
C / | |
C (k) (j) * Orientation Free Step
C /
C (k)EA
C Altnatively, |/
C *
C
MFF = (2*M1*M1*M2)**2
CALL MATCLR( MFF*MFF, EB )
DO 10 I = 0, MFF-1
DO 10 J = 0, MFF-1
DO 10 K = 0, MFF-1
EB(J*MFF+K)= EB(J*MFF+K) + EA(I*MFF+J)*EA(I*MFF+K)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE STEP8( M1, M2, EA, EB, EC )
INTEGER MFF, M1, M2
REAL*8 EA(0:*), EB(0:*), EC(0:*)
C
C ***** Extension of Proc8 (Orientation Free Step) *****
C
MFF = (2*M1*M1*M2)**2
CALL MATCLR( MFF*MFF, EC )
DO 10 I = 0, MFF-1
DO 10 J = 0, MFF-1
DO 10 K = 0, MFF-1
EC(J*MFF+K)= EC(J*MFF+K) + EA(I*MFF+J)*EB(I*MFF+K)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE CREZ( M1,M2, Z, E )
INTEGER MFF, M1,M2
REAL*8 Z, E(0:*)
MFF = (2*M1*M1*M2)**2
Z = 0
DO 10 I = 0, MFF-1
Z = Z + E(I*MFF+I)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE CRDM2( M1, M2, E3, DM2 )
INTEGER MF, M1, M2
REAL*8 E3(0:*), DM2(0:*)
MF = 2*M1*M1*M2
CALL MATCLR( MF*MF, DM2 )
DO 10 I = 0, MF-1
DO 10 J = 0, MF-1
DO 10 K = 0, MF-1
DM2(J*MF+K)= DM2(J*MF+K) + E3(((I*MF+J)*MF+I)*MF+K)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE CRDM1( M1, M2, DM2, DM1 )
INTEGER M1, M2
REAL*8 DM2(0:*), DM1(0:*)
CALL MATCLR( 2*M1 * 2*M1, DM1 )
DO 10 I = 0, M1*M2-1
DO 10 J = 0, 2*M1 -1
DO 10 K = 0, 2*M1 -1
DM1(J*2*M1+K)= DM1(J*2*M1+K) + DM2(((J*M1*M2+I)*2*M1+K)*M1*M2+I)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE DIAG( N, A, E, W )
INTEGER N, IERR
REAL*8 A(0:*), E(0:*), W(0:*), WK(0:1000)
C
IF( N .GT. 1000 ) THEN
Write(6,*) '--- Matrix Size Over in Diag ---'
STOP
END IF
C
CALL DCSMAA( A, N, N, W, WK, IERR )
IF( IERR .NE. 0 ) THEN
Write(6,*) '--- Error Number ',IERR,' in DCSMAA ---'
STOP
END IF
C
CALL REORD( N, A, E, W )
C Write(6,*)
C DO 10 I = 0, 3
C Write(6,*) W(I)/W(0)
C 10 CONTINUE
RETURN
END
CCC
SUBROUTINE KPTWT( N1, N2, KP, W )
INTEGER N1, N2
REAL*8 S1, S2, KP, W(0:*)
S1 = 0.0D0
DO 10 I = 0, N1-1
S1 = S1 + W(I)
10 CONTINUE
S2 = 0.0D0
DO 20 I = 0, N2-1
S2 = S2 + W(I)
20 CONTINUE
KP = S2 / S1
C Write(6,*) '- Kept Weight -', KP
RETURN
END
CCC
SUBROUTINE REORD( N, A, E, W )
INTEGER N
REAL*8 A(0:*), E(0:*), W(0:*)
C
DO 10 I = 0, N-1
E(I)= W(N-I-1)
10 CONTINUE
CALL MATCPY( N, E, W )
C
DO 20 I = 0, N-1
DO 20 J = 0, N-1
E(I*N+J)= A((N-I-1)*N+J)
20 CONTINUE
CALL MATCPY( N*N, E, A )
RETURN
END
CCC
SUBROUTINE PARITY( M1, M2, A, PM, W )
INTEGER MF, M1, M2
REAL*8 T(0:1024), A(0:*), PM(0:*), W(0:*)
MF = 2*M1*M1*M2
DO 10 I = 0, MF-1
T(I)= 0.0D0
DO 10 J = 0, 2-1
DO 10 K = 0, M1-1
DO 10 L = 0, M1-1
DO 10 M = 0, M2-1
T(I)= T(I) + A( I*MF + ((J*M1+K)*M1+L)*M2+M )
& * A( I*MF + ((J*M1+L)*M1+K)*M2+M ) * PM(M)
10 CONTINUE
C
DO 20 I = 0, 3
Write(6,100) W(I)/W(0), T(I)
100 Format(1H , 2F15.7 )
20 CONTINUE
C
CALL MATCPY( MF, T, PM )
C
RETURN
END
CCC
SUBROUTINE TRPOS( M1, M2, PM, DM, RM )
INTEGER MF, M1, M2
REAL*8 PM(0:*), DM(0:*), RM(0:*)
MF = 2*M1*M1*M2
DO 10 I = 0, MF-1
DO 10 J = 0, 2-1
DO 10 K = 0, M1-1
DO 10 L = 0, M1-1
DO 10 M = 0, M2-1
RM( I*MF + ((J*M1+K)*M1+L)*M2+M ) =
& DM( I*MF + ((J*M1+L)*M1+K)*M2+M )
C & DM( I*MF + ((J*M1+L)*M1+K)*M2+M ) * PM(M)
10 CONTINUE
RETURN
END
CCC
SUBROUTINE RENOB( M1, M1N, DM, B2, B1 )
INTEGER M1, M1N, R, AD2, AD1
REAL*8 DM(0:*), B2(0:*), B1(0:*)
C
CALL MATCLR( 2* 2*M1 * 2*M1 * 2*M1 * M1N, B1 )
DO 10 I = 0, 2 -1
DO 10 J = 0, 2*M1-1
DO 10 K = 0, 2*M1-1
DO 10 L = 0, 2*M1-1
DO 10 M = 0, 2*M1-1
DO 10 R = 0, M1N-1
AD2 =(((I*2*M1+J)*2*M1+K)*2*M1+L)*2*M1+M
AD1 =(((I*2*M1+J)*2*M1+K)*2*M1+L)*M1N +R
B1(AD1)= B1(AD1) + B2(AD2) * DM(R*2*M1+M)
10 CONTINUE
C
CALL MATCLR( 2* 2*M1 * 2*M1 * M1N * M1N, B2 )
DO 20 I = 0, 2 -1
DO 20 J = 0, 2*M1-1
DO 20 K = 0, 2*M1-1
DO 20 L = 0, 2*M1-1
DO 20 R = 0, M1N-1
DO 20 M = 0, M1N-1
AD1 =(((I*2*M1+J)*2*M1+K)*2*M1+L)*M1N +M
AD2 =(((I*2*M1+J)*2*M1+K)*M1N +R)*M1N +M
B2(AD2)= B2(AD2) + B1(AD1) * DM(R*2*M1+L)
20 CONTINUE
C
CALL MATCLR( 2* 2*M1 * M1N * M1N * M1N, B1 )
DO 30 I = 0, 2 -1
DO 30 J = 0, 2*M1-1
DO 30 K = 0, 2*M1-1
DO 30 R = 0, M1N-1
DO 30 L = 0, M1N-1
DO 30 M = 0, M1N-1
AD2 =(((I*2*M1+J)*2*M1+K)*M1N +L)*M1N +M
AD1 =(((I*2*M1+J)*M1N +R)*M1N +L)*M1N +M
B1(AD1)= B1(AD1) + B2(AD2) * DM(R*2*M1+K)
30 CONTINUE
C
CALL MATCLR( 2* M1N * M1N * M1N * M1N, B2 )
DO 40 I = 0, 2 -1
DO 40 J = 0, 2*M1-1
DO 40 R = 0, M1N-1
DO 40 K = 0, M1N-1
DO 40 L = 0, M1N-1
DO 40 M = 0, M1N-1
AD1 =(((I*2*M1+J)*M1N +K)*M1N +L)*M1N +M
AD2 =(((I*M1N +R)*M1N +K)*M1N +L)*M1N +M
B2(AD2)= B2(AD2) + B1(AD1) * DM(R*2*M1+J)
40 CONTINUE
C
CALL MATCPY( 2 * M1N * M1N * M1N * M1N, B2, B1 )
RETURN
END
CCC
SUBROUTINE RENOC( M1,M1N, M2,M2N, DM1, DM2, RM2, C3, C1 )
INTEGER R, MF, M1,M1N, M2,M2N, A3, A1
REAL*8 DM1(0:*), DM2(0:*), C3(0:*),RM2(0:*), C1(0:*)
C
MF = 2*M1*M1*M2
CALL MATCLR( 2*M1 * 2*M1 * MF * M2N, C1 )
DO 10 I = 0, 2*M1-1
DO 10 J = 0, 2*M1-1
DO 10 K = 0, MF-1
DO 10 L = 0, MF-1
DO 10 R = 0, M2N-1
A3 =((I*2*M1+J)*MF +K)*MF +L
A1 =((I*2*M1+J)*MF +K)*M2N+R
C1(A1)= C1(A1) + C3(A3) * DM2(R*MF+L)
C *Normal*
10 CONTINUE
C
CALL MATCLR( 2*M1 * 2*M1 * M2N * M2N, C3 )
DO 20 I = 0, 2*M1-1
DO 20 J = 0, 2*M1-1
DO 20 K = 0, MF-1
DO 20 R = 0, M2N-1
DO 20 L = 0, M2N-1
A1 =((I*2*M1+J)*MF +K)*M2N+L
A3 =((I*2*M1+J)*M2N+R)*M2N+L
C3(A3)= C3(A3) + C1(A1) * RM2(R*MF+K)
C *Reflct*
20 CONTINUE
C
CALL MATCLR( 2*M1 * M1N * M2N * M2N, C1 )
DO 30 I = 0, 2*M1-1
DO 30 J = 0, 2*M1-1
DO 30 R = 0, M1N-1
DO 30 K = 0, M2N-1
DO 30 L = 0, M2N-1
A3 =((I*2*M1+J)*M2N+K)*M2N+L
A1 =((I*M1N +R)*M2N+K)*M2N+L
C1(A1)= C1(A1) + C3(A3) * DM1(R*2*M1+J)
30 CONTINUE
C
CALL MATCLR( M1N * M1N * M2N * M2N, C3 )
DO 40 I = 0, 2*M1-1
DO 40 R = 0, M1N-1
DO 40 J = 0, M1N-1
DO 40 K = 0, M2N-1
DO 40 L = 0, M2N-1
A1 =((I*M1N +J)*M2N+K)*M2N+L
A3 =((R*M1N +J)*M2N+K)*M2N+L
C3(A3)= C3(A3) + C1(A1) * DM1(R*2*M1+I)
40 CONTINUE
C
CALL MATCPY( M1N * M1N * M2N * M2N, C3, C1 )
RETURN
END
CCC
SUBROUTINE RENOD( M1, M2, M2N, DM, RM, D2, D1 )
INTEGER MF, M1, M2, M2N, R, AD2, AD1
REAL*8 RM(0:*), DM(0:*), D2(0:*), D1(0:*)
C
MF = 2*M1*M1*M2
CALL MATCLR( MF*MF*M2N, D1 )
DO 10 I = 0, MF -1
DO 10 J = 0, MF -1
DO 10 K = 0, MF -1
DO 10 R = 0, M2N-1
AD2 = (I*MF +J)*MF +K
AD1 = (I*MF +J)*M2N+R
D1(AD1)= D1(AD1) + D2(AD2) * RM(R*MF+K)
C *Rflct*
10 CONTINUE
C
CALL MATCLR( MF*M2N*M2N, D2 )
DO 20 I = 0, MF -1
DO 20 J = 0, MF -1
DO 20 R = 0, M2N-1
DO 20 K = 0, M2N-1
AD1 = (I*MF +J)*M2N+K
AD2 = (I*M2N+R)*M2N+K
D2(AD2)= D2(AD2) + D1(AD1) * DM(R*MF+J)
C *Norml*
20 CONTINUE
C
CALL MATCLR( M2N*M2N*M2N, D1 )
DO 30 I = 0, MF -1
DO 30 R = 0, M2N-1
DO 30 J = 0, M2N-1
DO 30 K = 0, M2N-1
AD2 = (I*M2N+J)*M2N+K
AD1 = (R*M2N+J)*M2N+K
D1(AD1)= D1(AD1) + D2(AD2) * RM(R*MF+I)
C *Rflct*
30 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE NORM( L, V, W )
INTEGER L
REAL*8 V(0:*), W
C
CALL FMAX( L, V, W )
DO 20 I = 0, L-1
V(I)= V(I) / W
20 CONTINUE
RETURN
END
CCC
SUBROUTINE NRA2( N, A, EA )
INTEGER N
REAL*8 A(0:*), EA
DO 10 I = 0, N-1
A(I)= A(I) / EA
10 CONTINUE
RETURN
END
CCC
SUBROUTINE FMAX( L, V, W )
INTEGER L
REAL*8 V(0:*), W
W = 0.0D0
DO 10 I = 0, L-1
IF( ABS(V(I)) .GT. W ) W = ABS(V(I))
10 CONTINUE
RETURN
END
CCC
SUBROUTINE SYMCHA( A )
REAL*8 W, A(0:*)
DO 10 I = 0, 1
DO 10 J = 0, 1
DO 10 K = 0, 1
DO 10 L = 0, 1
DO 10 M = 0, 1
DO 10 N = 0, 1
W = ABS( A( ((((I*2+J)*2+K)*2+L)*2+M)*2+N )
& -A( ((((J*2+K)*2+L)*2+M)*2+N)*2+I ) )
IF( W .GT. 1.0E-10 ) THEN
Write(6,*)'-A-',I,J,K,L,M,N,A(((((I*2+J)*2+K)*2+L)*2+M)*2+N)
Write(6,*)'-A-',J,K,L,M,N,I,A(((((J*2+K)*2+L)*2+M)*2+N)*2+I)
STOP
END IF
10 CONTINUE
RETURN
END
CCC
SUBROUTINE SYMCHB( M1, B )
INTEGER EC, M1
REAL*8 FX, W, B(0:*)
C
CALL FMAX( 2*M1*M1*M1*M1, B, FX )
C Write(6,*) 'Maximum Element in B = ',FX
EC = 0
DO 10 I = 0, 2-1
DO 10 J = 0, M1-1
DO 10 K = 0, M1-1
DO 10 L = 0, M1-1
DO 10 M = 0, M1-1
W = ABS( B( (((I*M1+J)*M1+K)*M1+L)*M1+M )
& - B( (((I*M1+K)*M1+L)*M1+M)*M1+J ) )
& + ABS( B( (((I*M1+K)*M1+L)*M1+M)*M1+J )
& - B( (((I*M1+L)*M1+M)*M1+J)*M1+K ) )
& + ABS( B( (((I*M1+L)*M1+M)*M1+J)*M1+K )
& - B( (((I*M1+M)*M1+J)*M1+K)*M1+L ) )
& + ABS( B( (((I*M1+M)*M1+J)*M1+K)*M1+L )
& - B( (((I*M1+J)*M1+K)*M1+L)*M1+M ) )
IF( W/FX .GT. 1.0E-10 .AND. EC .LE. 10 ) THEN
C Write(6,*) '-B- ',I,J,K,L,M,' = ', W
Write(6,*) ' ',I,J,K,L,M,' ',
& B( (((I*M1+J)*M1+K)*M1+L)*M1+M )
Write(6,*) ' ',I,K,L,M,J,' ',
& B( (((I*M1+K)*M1+L)*M1+M)*M1+J )
Write(6,*) ' ',I,L,M,J,K,' ',
& B( (((I*M1+L)*M1+M)*M1+J)*M1+K )
Write(6,*) ' ',I,M,J,K,L,' ',
& B( (((I*M1+M)*M1+J)*M1+K)*M1+L )
EC = EC + 1
C IF( EC .GT. 10 ) STOP
END IF
10 CONTINUE
RETURN
END
CCC
SUBROUTINE SYMCHC( M1, M2, C )
INTEGER EC, M1, M2
REAL*8 FX, W, C(0:*)
C
CALL FMAX( M1*M1*M2*M2, C, FX )
C Write(6,*) 'Maximum Element in C = ',FX
EC = 0
DO 10 I = 0, M1-1
DO 10 J = 0, M1-1
DO 10 K = 0, M2-1
DO 10 L = 0, M2-1
W = ABS( ABS(C(((I*M1+J)*M2+K)*M2+L))
& -ABS(C(((J*M1+I)*M2+K)*M2+L)) )
& + ABS( ABS(C(((I*M1+J)*M2+K)*M2+L))
& -ABS(C(((I*M1+J)*M2+L)*M2+K)) )
IF( W/FX .GT. 1.0E-10 .AND. EC .LE. 10 ) THEN
C Write(6,*) '-C- ',I,J,K,L,' = ', W
Write(6,*) ' ',I,J,K,L,' ', C(((I*M1+J)*M2+K)*M2+L)
Write(6,*) ' ',J,I,K,L,' ', C(((J*M1+I)*M2+K)*M2+L)
Write(6,*) ' ',I,J,L,K,' ', C(((I*M1+J)*M2+L)*M2+K)
EC = EC + 1
C IF( EC .GT. 10 ) STOP
END IF
10 CONTINUE
RETURN
END
CCC
SUBROUTINE SYMCHD( M2, D )
INTEGER EC, M2
REAL*8 FX, W, D(0:*)
C
CALL FMAX( M2*M2*M2, D, FX )
C Write(6,*) 'Maximum Element in D = ',FX
EC = 0
DO 10 I = 0, M2-1
DO 10 J = 0, M2-1
DO 10 K = 0, M2-1
W = ABS( ABS(D((I*M2+J)*M2+K)) - ABS(D((J*M2+K)*M2+I)) )
& + ABS( ABS(D((J*M2+K)*M2+I)) - ABS(D((K*M2+I)*M2+J)) )
& + ABS( ABS(D((K*M2+I)*M2+J)) - ABS(D((I*M2+J)*M2+K)) )
IF( W/FX .GT. 1.0E-10 .AND. EC .LE. 10 ) THEN
C Write(6,*) '-D- ',I,J,K,' = ', W
Write(6,*) ' ',I,J,K,' ', D((I*M2+J)*M2+K)
Write(6,*) ' ',J,K,I,' ', D((J*M2+K)*M2+I)
Write(6,*) ' ',K,I,J,' ', D((K*M2+I)*M2+J)
EC = EC + 1
C IF( EC .GT. 10 ) STOP
END IF
10 CONTINUE
RETURN
END
CCC
SUBROUTINE SYMCDM( M1, M2, PM, DM )
INTEGER EC, M1, M2, Q, R, S, T
REAL*8 FX, W, PM(0:*), DM(0:*)
CALL FMAX((2*M1*M1*M2)**2, DM, FX )
C Write(6,*) 'Maximum Element in DM = ',FX
EC = 0
DO 10 I = 0, 2-1
DO 10 J = 0, M1-1
DO 10 K = 0, M1-1
DO 10 L = 0, M2-1
DO 10 Q = 0, 2-1
DO 10 R = 0, M1-1
DO 10 S = 0, M1-1
DO 10 T = 0, M2-1
W = DM(((((((I*M1+J)*M1+K)*M2+L)*2+Q)*M1+R)*M1+S)*M2+T)
& - DM(((((((I*M1+K)*M1+J)*M2+L)*2+Q)*M1+S)*M1+R)*M2+T)
& *PM(L) *PM(T)
IF( ABS(W)/FX .GT. 1.0E-10 .AND. EC .LE. 10 ) THEN
Write(6,*) '-DM ',I,J,K,L,Q,R,S,T,' = ', W
EC = EC + 1
C IF( EC .GT. 10 ) STOP
END IF
10 CONTINUE
RETURN
END
CCC
CC
C
SUBROUTINE MATCLR( N, A )
INTEGER N
REAL*8 A(0:*)
DO 10 I = 0, N-1
A(I)= 0.0D0
10 CONTINUE
RETURN
END
CCC
SUBROUTINE MATCPY( N, A, B )
INTEGER N
REAL*8 A(0:*), B(0:*)
DO 10 I = 0, N-1
B(I)= A(I)
10 CONTINUE
RETURN
END