Tensor Product Variational Appriximation for 3D Ising Model
c123456789012345678901234567890123456789012345678901234567890123456789012
c ----- -----
c ----- IRF-variational approximation for 3D Ising model -----
c ----- assisted by CTMRG. -----
c ----- -----
c ----- '99/12/22, Ver.1.1 by T.Nishino -----
c ----- -----
module var_size
integer,parameter :: iq = 2, mVL = iq**8, mXX=20
integer,parameter :: mIT =10000, mXT = mVL*mXX**2
end module
c
program main
use var_size
integer mX, L, LL, m(0:1), mN(0:1), q(0:1), i, j, N, info
real(8) K, fP(0:1), G1(0:mXT,0:1)
real(8) x(0:3), fC(0:1), W(0:mVL,0:3), G2(0:mXT,0:1)
real(8) y(0:3), eP(0:1), P(0:mXT,0:1), G3(0:mXT,0:1)
real(8) u, eC(0:1), C(0:mXT,0:1), G4(0:mXT,0:1)
real(8) v, Z(0:1,0:mIT), F(0:mIT,0:1), G5(0:mXT,0:1)
real(8) eps, mg(0:1,0:mIT)
c
call getprm( L, mX, K, eps ) ! --- Parameter Imput ---
call initw( q, m, K, W, P, C ) ! --- Initialize TM ---
c
LL = 2; eP(0:1)= 0.0d0; eC(0:1)= 0.0d0
goto 20
c
10 continue
call joint1( 0,0,1, q, m, C, P, G1 ) ! --- Conserve G1 ---
call joint1( 1,0,1, q, m, P, G1, G2 )
call joint1( 0,0,2, q, m, C, G2, G3 )
call jointX( q, m, G3, G2 ) ! --- Conserve G2 ---
call obsrv1( x,y, q, W, G2 )
c
c -------------- Improve IRF weight -----------------------
c
G3(0:q(1)**4-1,1) = G2(0:q(1)**4-1,1)*W(0:q(1)**4-1,3)
c
call atorame( q(0), G3(0,1), G4(0,1) )
N = q(0)**4
G3(0:N-1,1) = 0.0d0
do i = 0, N-1
G3(0:N-1,1)= G3(0:N-1,1) + G4(i*N:i*N+N-1,1)*W(i,2)
end do
c
G3(0:N-1,1) = G3(0:N-1,1) / G2(0:N-1,0)
c
call norm( N, G3(0,1), u )
c
call vratio( q, u, G2, W )
c
W(0:N-1,2)= W(0:N-1,2) + eps * G3(0:N-1,1)
c
call norm( N, W(0,2), v )
call bigmac( q, W )
call vratio( q, v, G2, W )
c
c ----------------------------------------- Create DMat ---
call extP( q, m, W, P, G2 ) ! --- Conserve G2 ---
call extC( q, m, G1, G2, G3 ) ! --- Conserve G3 ---
call joint3( q, m, G3, G4 ) ! --- Expire G1 ---
call joint3( q, m, G4, G5 )
call obsrv2( q, m, mg(0,LL), G5 )
call credm( q, m, G5, G4 ) ! --- Dens. Mat. ---
call crepf( q, m, Z(0,LL), G4 )
c
F(LL,0:1)= eP(0:1)*8 + eC(0:1)*4 + log(Z(0:1,LL))
c
c ------------------------------------- Display Results ---
if(LL .gt. 8)
& write(6,100) LL, v/u-1.0d0,
& mg(0,LL )/Z(0,LL ), mg(1,LL )/Z(1,LL ), y(0),
& ( F(LL,0)+F(LL-4,0)-2*F(LL-2,0)
& -F(LL,1)-F(LL-4,1)+2*F(LL-2,1) )/8
100 format('#',I4,E11.4, F15.10, F15.10, F12.7, F15.10)
c
call diag( q, m, G4, G1, G5 )
c
mN(0:1)=min( mX, q(0:1)*m(0:1) )
c
call renop( q, m, mN, G2, P, G4 ) ! --- Expire G2 ---
call renoc( q, m, mN, G3, C, G4 ) ! --- Expire G3 ---
m(0:1) = mN(0:1)
c
call symchp( q(0), m(0), P(0,0), ' -a0- ' ) ! #
call symchp( q(1), m(1), P(0,1), ' -a1- ' ) ! Debug
call symchc( q(0), m(0), C(0,0), ' -a0- ' ) ! Routines
call symchc( q(1), m(1), C(0,1), ' -a1- ' ) ! #
c
20 continue
call norm( q(0)**2 * m(0)**2, P(0,0), fP(0) )
call norm( q(1)**2 * m(1)**2, P(0,1), fP(1) )
call norm( q(0) * m(0)**2, C(0,0), fC(0) )
call norm( q(1) * m(1)**2, C(0,1), fC(1) )
c
eC(0:1)= eC(0:1) + 2*eP(0:1) + log(fC(0:1))
eP(0:1)= eP(0:1) + log(fP(0:1))
LL = LL + 2
if( LL .LT. L ) goto 10
end program
ccccc
subroutine errstp( ch )
character(6) ch
write(6,*) '*** Error stop in subroutine -',ch,'- ***'
stop "Error_stop"
end
ccccc
subroutine getprm( L, mX, K, eps )
use var_size
integer L, mX; real(8) K, eps
c
write(6,*) 'Input Max. Iteration # L'; read(5,*) L
if( L .GT. mIT ) call errstp('L >mIT')
c
write(6,*) 'Input Maximum Dim. mX>= m'; read(5,*) mX
if( mX .GT. mXX ) call errstp('mX>mXX')
c
write(6,*) 'Input improvement fac.eps'; read(5,*) eps
c
write(6,*) 'Input J/T K '; read(5,*) K
write(6,*) 'L = ',L,' mX = ',mX, ' K = ',K
write(6,*) 'eps= ', eps
end
ccccc
subroutine initw( q, m, K, W, P, CC )
use var_size
integer q(0:*), m(0:*), a, b, c, d
real(8) K, f(0:1,0:1), W(0:mVL,0:*), P(0:mXT,0:*), CC(0:mXT,0:*)
c
integer :: s1(0:3)=(/0,1,0,1/), s2(0:3)=(/0,0,1,1/)
c
f(0,0)= exp( K/4 ); f(0,1)= exp( -K/4 ) ! --- Bond Factors ---
f(1,1)= exp( K/4 ); f(1,0)= exp( -K/4 )
c
q(0)= iq
do a = 0, q(0) - 1 ! d----c
do b = 0, q(0) - 1 ! | W2 | = Initial Face Factor
do c = 0, q(0) - 1 ! a----b
do d = 0, q(0) - 1
W(((a*q(0)+b)*q(0)+c)*q(0)+d,2)= f(a,b)*f(b,c)*f(c,d)*f(d,a)
end do; end do; end do; end do
c
q(1)= iq**2
do a = 0, q(1) - 1 ! d----c
do b = 0, q(1) - 1 ! | W3 | = Cube Factor
do c = 0, q(1) - 1 ! a----b
do d = 0, q(1) - 1
W((( a *q(1)+ b) *q(1)+ c) *q(1)+ d, 3)=
& W(((s1(a)*q(0)+s1(b))*q(0)+s1(c))*q(0)+s1(d),2)*
& W(((s2(a)*q(0)+s2(b))*q(0)+s2(c))*q(0)+s2(d),2)*
& f( s1(a), s2(a))* f(s1(b), s2(b) )*
& f( s1(c), s2(c))* f(s1(d), s2(d) )
end do; end do; end do; end do
c
call bigmac( q, W )
c
do a = 0, 1
m(a)=1 ! --- Fixed Boundary Condition ---
do b = 0, q(a)-1
do c = 0, q(a)-1 ! 0------c 0------c
P(b*q(a)+c,a) ! | P0/1 | = | W0/1 |
& =W(b*q(a)+c,a) ! 0------b 0------b
end do; end do ! (use rotated W)
do b = 0, q(a)-1 ! +------0 0------0
CC(b,a)= W(b,a) ! | C0/1 | = | W0/1 |
end do ! 0------b 0------b
end do ! (use rotated W)
end
ccccc
subroutine bigmac( q, W )
use var_size
integer i, j, q(0:*); real(8) W(0:mVL,0:*)
real(8),allocatable :: tp(:)
c ! d----c d----c d----c
W(0:q(0)**4-1,0) = ! | W0 | = | W2 | * | W2 |
& W(0:q(0)**4-1,2)**2 ! a----b a----b a----b
c
allocate( tp(0:q(0)**8-1) ); call atorame( q(0), W(0,3), tp )
c
do i = 0, q(0)**4-1 ! d1---c1 d12--c12 d2---c2 d----c
do j = 0, q(0)**4-1 ! | W2 | * | W3 | * | W2 | = | W1 |
tp(i*q(0)**4+j)= W(i,2)! a1---b1 a12--b12 a2---b2 a----b
& * tp(i*q(0)**4+j)* W(j,2)
end do; end do
c
call torame( q(0), tp, W(0,1) ); deallocate( tp )
end
ccccc
subroutine torame( q, W, Y )
integer a,b,c,d,e,f,g,h, q; real(8) W(0:*), Y(0:*)
c
do a = 0, q-1; do b = 0, q-1; do c = 0, q-1; do d = 0, q-1
do e = 0, q-1; do f = 0, q-1; do g = 0, q-1; do h = 0, q-1
c
Y(((((((a*q+e)*q+b)*q+f)*q+c)*q+g)*q+d)*q+h)
& = W(((((((a*q+b)*q+c)*q+d)*q+e)*q+f)*q+g)*q+h)
end do; end do; end do; end do; end do; end do; end do; end do
end
ccccc
subroutine atorame( q, W, Y )
integer a,b,c,d,e,f,g,h, q; real(8) W(0:*), Y(0:*)
c
do a = 0, q-1; do b = 0, q-1; do c = 0, q-1; do d = 0, q-1
do e = 0, q-1; do f = 0, q-1; do g = 0, q-1; do h = 0, q-1
c
Y(((((((a*q+b)*q+c)*q+d)*q+e)*q+f)*q+g)*q+h)
& = W(((((((a*q+e)*q+b)*q+f)*q+c)*q+g)*q+d)*q+h)
end do; end do; end do; end do; end do; end do; end do; end do
end
ccccc
subroutine joint1( p1, p2, p3, q, m, A, B, C )
use var_size
integer i,j,k,l,n, p1, p2, p3, q(0:*), m(0:*)
real(8) A(0:mXT,0:*), B(0:mXT,0:*), C(0:mXT,0:*)
c
do i = 0, 1
C(0:m(i)**2 *q(i)**(p1+p3+1)-1,i)= 0.0d0
do j = 0, m(i)*q(i)** p1-1
do k = 0, q(i) -1 ! +-----L-----+
do l = 0, m(i)*q(i)** p2-1 ! | 2l |
do n = 0, m(i)*q(i)** p3-1 ! | 1 | 3 |
C((j*q(i)+k)*q(i)** p3*m(i)+n,i) = ! J-j---k---n-N
& C((j*q(i)+k)*q(i)** p3*m(i)+n,i) +
& A((j*q(i)+k)*q(i)** p2*m(i)+l,i) *
& B((l*q(i)+k)*q(i)** p3*m(i)+n,i)
end do; end do; end do; end do
end do
end
ccccc
subroutine jointX( q, m, A, B )
use var_size
integer i,j,k,l,n, p1, p2, q(0:*), m(0:*)
real(8) A(0:mXT,0:*), B(0:mXT,0:*)
do i = 0, 1
B(0:q(i)**4-1,i)= 0.0d0 ! +--p1 p1-----+
do j = 0, q(i)-1 ! | | | |
do k = 0, q(i)-1 ! | l l--n |
do l = 0, q(i)-1 ! | | | |
do n = 0, q(i)-1 ! | k--j j |
do p1 = 0, m(i)-1 ! | | | |
do p2 = 0, m(i)-1 ! +-----p2 p2--+
c
B( (( j *q(i)+k)*q(i)+l)*q(i)+ n,i)
& = B( (( j *q(i)+k)*q(i)+l)*q(i)+ n,i)
& + A((((p2*q(i)+j)*q(i)+k)*q(i)+l)*m(i)+p1,i)
& * A((((p1*q(i)+l)*q(i)+n)*q(i)+j)*m(i)+p2,i)
c
end do; end do; end do; end do; end do; end do
end do
end
ccccc
subroutine joint3( q, m, A, B )
use var_size
integer q(0:*), m(0:*), i, c, d, e, f
real(8) A(0:mXT,0:*), B(0:mXT,0:*)
c
do i = 0, 1 ! +---f
B(0:q(i)**3*m(i)**2-1,i)=0 ! | A |
do c = 0, q(i)* m(i)-1 ! e---d
do d = 0, q(i) -1 ! | A |
do e = 0, q(i)* m(i)-1 ! +---c
do f = 0, q(i)* m(i)-1
B((c*q(i)+d)*q(i)*m(i)+f,i) = B((c*q(i)+d)*q(i)*m(i)+f,i)
& + A((c*q(i)+d)*q(i)*m(i)+e,i)
& * A((e*q(i)+d)*q(i)*m(i)+f,i)
end do; end do; end do; end do
end do
end
ccccc
subroutine obsrv1( x, y, q, W, G )
use var_size
integer :: i, q(0:*), dg(0:4)=(/1,2,4,8,1/)
real(8) z, x(0:*), y(0:*), W(0:mVL,0:*), G(0:mXT,0:*)
c
x(0:3) = 0.0d0; y(0:3) = 0.0d0; z = 0.0d0
c
do i = 0, q(0)**4-1
z = z + G(i,0)*W(i,0)
x(0:3)= x(0:3) + G(i,0)*W(i,0) * (mod(i/dg(0:3),2)*2-1)
y(0:3)= y(0:3) + G(i,0)*W(i,0) * (mod(i/dg(0:3),2)*2-1)
& * (mod(i/dg(1:4),2)*2-1)
end do
x(0:3) = x(0:3) / z
y(0:3) = y(0:3) / z
end
ccccc
subroutine obsrv2( q, m, mg, G )
use var_size
integer i,j,k, q(0:*), m(0:*)
real(8) mg(0:*), G(0:mXT,0:*)
c
mg(0:1)= 0.0d0
do i = 0, 1
do j = 0, q(i) -1
do k = 0, q(i)*m(i)-1
mg(i)= mg(i) + G((k*q(i)+j)*q(i)*m(i)+k,i) * (mod(j,2)*2-1)
end do; end do; end do
end
ccccc
subroutine vratio( q, u, G2, W )
use var_size
integer i, j, q(0:*)
real(8) a(0:1), u, G2(0:mXT,0:*), W(0:mVL,0:*)
c
do i = 0, 1
a(i)= 0.0d0
do j = 0, q(i)**4-1
a(i)= a(i) + W(j,i)*G2(j,i)
end do
end do
u = a(1)/a(0)
end
ccccc
subroutine extP( q, m, W, P, G )
use var_size
integer q(0:*), m(0:*), a, b, c, d, e, f, i
real(8) W(0:mVL,0:*), P(0:mXT,0:*), G(0:mXT,0:*)
c
do i = 0, 1
do a = 0, m(i)-1 ! f----e---d
do b = 0, q(i)-1 ! | P | W |
do c = 0, q(i)-1 ! a----b---c
do d = 0, q(i)-1 !
do e = 0, q(i)-1 ! careful to the a-b order
do f = 0, m(i)-1
G(((((b*m(i)+a)*q(i)+c)*q(i)+d)*q(i)+e)*m(i)+f,i) =
& P(((a*q(i)+b)*q(i) +e)*m(i)+f,i) *
& W(((b *q(i)+c)*q(i)+d)*q(i)+e, i)
end do; end do; end do; end do; end do; end do
end do
end
ccccc
subroutine extC( q, m, G1, G2, G3 )
use var_size
integer q(0:*), m(0:*), a, b, c, d, e, f, i
real(8) G1(0:mXT,0:*), G2(0:mXT,0:*), G3(0:mXT,0:*)
c
do i = 0, 1
G3(0:q(i)**3*m(i)**2-1,i) = 0
do a = 0, q(i) *m(i)-1 ! +-------f
do b = 0, q(i) -1 ! | G1 |
do c = 0, q(i) -1 ! e---d---c
do d = 0, q(i) -1 ! | G2 |
do e = 0, m(i)-1 ! A---a---b
do f = 0, m(i)-1
G3( ((a*q(i)+b)*q(i)+c) *m(i)+f,i) =
& G3( ((a*q(i)+b)*q(i)+c) *m(i)+f,i) +
& G2((((a*q(i)+b)*q(i)+c)*q(i)+d)*m(i)+e,i) *
& G1( ((e*q(i)+d)*q(i)+c) *m(i)+f,i)
end do; end do; end do; end do; end do; end do
end do
end
ccccc
subroutine credm( q, m, A, B )
use var_size
integer q(0:*), m(0:*), i, c, d, e
real(8) A(0:mXT,0:*), B(0:mXT,0:*)
c
do i = 0, 1
B(0:q(i)**2*m(i)**2-1,i)=0
do c = 0, q(i)* m(i)-1
do d = 0, q(i) -1
do e = 0, q(i)* m(i)-1
B(c*q(i)* m(i)+e,i)= B( c *q(i)*m(i)+e,i)
& + A((c*q(i)+d)*q(i)*m(i)+e,i)
end do; end do; end do
end do
end
ccccc
subroutine crepf( q, m, Z, A )
use var_size
integer q(0:*), m(0:*), i, j
real(8) Z(0:*), A(0:mXT,0:*)
do i = 0, 1
Z(i)= 0.0d0
do j = 0, q(i)*m(i)-1
Z(i)= A(j*q(i)*m(i)+j,i) + Z(i)
end do; end do
end
ccccc
subroutine diag( q, m, A, E, tp )
use var_size
integer q(0:*), m(0:*), N, info
real(8) A(0:mXT,0:*), E(0:mXT,0:*), tp(0:*)
c
do i = 0, 1
N = q(i)*m(i)
call DSYEV('V','U', N, A(0,i), N, E(0,i), tp, (N+2)*N, info )
c call dcsmaa( A(0,i), N, N, E(0,i), tp, info )
c
if( info .ne. 0 ) then
write(6,*) 'dsyev_info = ', info
call errstp(' diag ')
end if
c write(6,'(f30.10)') (E(j,i),j=N-1,N-min(N-1,8)-1,-1)
end do
end
cccccc
subroutine renop( q, m, mN, G, P, R )
use var_size
integer i, a, b, c, q(0:*), m(0:*), N, mN(0:*)
real(8) G(0:mXT,0:*), P(0:mXT,0:*), R(0:mXT,0:*)
c
do i = 0, 1
N = q(i) *m(i) ! c---c c---c
P(0:q(i)**2*N*mN(i)-1,i)= 0.0d0 ! | | | |
c ! c P | = c G b
do a = 0, mN(i)-1 ! | | | |R\
do b = 0, N -1 ! C---a C---B--a
do c = 0, q(i)**2*N-1
P(a*q(i)**2*N+c,i) = P(a*q(i)**2*N+c,i)
& + G(b*q(i)**2*N+c,i)
& * R( (N-a-1)*N+b,i)
end do; end do; end do
c ! a---a a---a
G(0:q(i)**2*mN(i)**2-1,i) = 0.0d0 ! | | | |
c ! | G | = c P |
do a = 0, q(i)**2*mN(i)-1 ! | | /R| |
do b = 0, mN(i)-1 ! b---a b--C---a
do c = 0, N -1
G(a*mN(i)+b,i) = G( a* mN(i)+b,i)
& + P( a* N +c,i)
& * R((N-b-1)*N+c,i)
end do; end do; end do
c
P(0:q(i)**2*mN(i)**2-1,i) = G(0:q(i)**2*mN(i)**2-1,i)
end do
end
ccccc
subroutine renoc( q, m, mN, G, C, R )
use var_size
integer i, a,b, d, q(0:*), m(0:*), N, mN(0:*)
real(8) G(0:mXT,0:*), C(0:mXT,0:*), R(0:mXT,0:*)
c
do i = 0, 1 ! d--d--d d--d--d
N = q(i) *m(i) ! | | | |
C(0:q(i)*N*mN(i)-1,i) = 0.0d0 ! | C | = | G b
do a = 0, mN(i)-1 ! | | | |R\
do b = 0, N -1 ! +-----a +-----b--a
do d = 0, q(i)*N -1
C(a*q(i)*N+d,i)= C( a*q(i)*N+d,i)
& + G( b*q(i)*N+d,i)
& * R((N-a-1)*N+b,i)
end do; end do; end do
c ! b
G(0:q(i)*mN(i)**2-1,i) = 0.0d0 ! |R\
do a = 0, q(i)*mN(i)-1 ! b-----a d--d--a
do b = 0, mN(i)-1 ! | G | = | C |
do d = 0, N -1 ! +-----a +-----a
G(a*mN(i)+b,i)= G( a* mN(i)+b,i)
& + C( a* N +d,i)
& * R((N-b-1)*N +d,i)
end do; end do; end do
c
C(0:q(i)*mN(i)**2-1,i) = G(0:q(i)*mN(i)**2-1,i)
end do
end
ccccc
subroutine norm( N, V, W )
integer i, N; real(8) V(0:*), W
W = 0.0d0
do i = 0, N-1
if( abs(V(i)) .GT. W ) W = abs(V(i))
end do
V(0:N-1)= V(0:N-1) / W
end
ccccc
subroutine symchp( q, m, P, ch )
integer i,j,k,l, q, m; real(8) P(0:*), S; character(6) ch
S = 0.0d0
do i = 0, m-1
do j = 0, q-1
do k = 0, q-1
do l = 0, m-1
S = S + ABS( P(((i*q+j)*q+k)*m+l)-P(((l*q+k)*q+j)*m+i) )
end do; end do; end do; end do
if( S .gt. 1.0E-8 ) write(6,*) 'sym P = ', S, ch
end
ccccc
subroutine symchc( q, m, C, ch )
integer i,j,k, q, m; real(8) C(0:*), S; character(6) ch
S = 0.0d0
do i = 0, m-1
do j = 0, q-1
do k = 0, m-1
S = S + ABS( C((i*q+j)*m+k)-C((k*q+j)*m+i) )
end do; end do; end do;
if( S .gt. 1.0E-8 ) write(6,*) 'sym C = ', S, ch
end