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