DMRG for 2D Ising model represented as q**4-vertex model
@
This program runs with NEC's ASL library
(Please contact me if you wish to know about ASL.)

c123456789012345678901234567890123456789012345678901234567890123456789012
c     -----                                                      -----
c     -----                                                      -----
c     -----     DMRG for q-state lattice polymer model           -----
c     -----                                                      -----
c     -----                      1999/3/18, Ver.0.9 by T.Nishino -----
c     -----                                                      -----
c     -----                                                      -----
c     -----                                                      -----
      parameter(iq=2,NX=200, mXX=100)
c
      integer  i, ii, j, k, L, N, m(0:NX), mX, LV, Bc, NV
      real*8   Z,Tp
      real*8   T(0:iq*mXX**2,0:NX), w(0:iq**4),s(0:iq**4), Mg(0:NX)
      real*8   A(0:iq*mXX**2,0:NX), v(0:iq**2*mXX**2,0:3), Fc(0:NX)
      real*8                        d(0:iq**3*mXX**2,0:2)
c
c     ----- Initialization -----
      call getprm( LV, Bc, mX, N, L, Tp )
      call initpt(      i,  j, N, m, Fc )
      call inittm( LV, Bc, Tp, s, w, T(0,i), T(0,j) )
c
c     <<<<<<<<<<<<<<<<<<<<<< infinite process >>>>>>>>>>>>>>>>>>>
      do 10 k = 1, N/2-2
      call psiini(     NV, m(i), m(j),                       V )
      call lanzos( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      call invers( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      if( LV.ge.1 )  call pdisp( i, j, N, Fc, Z )
      call exchg3(       m(i),iq,iq*m(j), V(0,1), V(0,2)         )
      call ljoint(    iq*m(i),   iq*m(j),         V(0,2), A(0,i) )
      call   diag( LV,iq*m(i),            V(0,2), V(0,3), A(0,i) )
      call hjoint(       m(i)*iq,iq*m(j), V(0,1),         A(0,j) )
      call   diag( LV,           iq*m(j), V(0,2), V(0,3), A(0,j) )
c
      call   extl( mX,m(i),m(i+1),T(0,i),T(0,i+1),w,A(0,i),Fc(i+1),d )
                        i  = i+1
      call   extr( mX,m(j),m(j-1),T(0,j),T(0,j-1),w,A(0,j),Fc(j-1),d )
                        j  = j-1
   10 continue
c 
      do 90 ii= 1, L
c                             finite process >>>>>>>>>>>>>>>>>>>>
      do 30 k = 1, N/2-2
      call psiini(     NV, m(i), m(j),                       V )
      call lanzos( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      call invers( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      if( LV.ge.1 )  call pdisp( i, j, N, Fc, Z )
      call exchg3(       m(i),iq,iq*m(j), V(0,1), V(0,2)         )
      call ljoint(    iq*m(i),   iq*m(j),         V(0,2), A(0,i) )
      call   diag( LV,iq*m(i),            V(0,2), V(0,3), A(0,i) )
c
      call   extl( mX,m(i),m(i+1),T(0,i),T(0,i+1),w,A(0,i),Fc(i+1),d )
                        i  = i+1
                        j  = j+1
   30 continue
c     <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      do 40 k = 1, N-4
      call psiini(     NV, m(i), m(j),                       V )
      call lanzos( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      call invers( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      if( LV.ge.1 )  call pdisp( i, j, N, Fc, Z )
      call hjoint(       m(i)*iq,iq*m(j), V(0,1),         A(0,j) )
      call   diag( LV,           iq*m(j), V(0,2), V(0,3), A(0,j) )
c
      call   extr( mX,m(j),m(j-1),T(0,j),T(0,j-1),w,A(0,j),Fc(j-1),d )
                        i  = i-1
                        j  = j-1
   40 continue
c     >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      do 50 k = 1, N/2-2
      call psiini(     NV, m(i), m(j),                       V )
      call lanzos( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      call invers( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      if( LV.ge.1 )  call pdisp( i, j, N, Fc, Z )
      call exchg3(       m(i),iq,iq*m(j), V(0,1), V(0,2)         )
      call ljoint(    iq*m(i),   iq*m(j),         V(0,2), A(0,i) )
      call   diag( LV,iq*m(i),            V(0,2), V(0,3), A(0,i) )
c
      call   extl( mX,m(i),m(i+1),T(0,i),T(0,i+1),w,A(0,i),Fc(i+1),d )
                        i  = i+1
                        j  = j+1
   50 continue
   90 continue
c
      call psiini(     NV, m(i), m(j),                       V )
      call lanzos( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      call invers( LV, NV, m(i), m(j), W, T(0,i), T(0,j), Z, V )
      if( LV.ge.1 )  call pdisp( i, j, N, Fc, Z )
      write(6,*) 'program end'
c
      end
ccc
      subroutine errstp( ch )
      character*6        ch
      write(6,*) '*** Error stop in subroutine -',ch,'- ***'
      stop
      end
ccc
      subroutine dump( ch, N, a  )
      character*4      ch
      integer           i, N
      real*8                  a(0:*)
      write(6,*) ch
      write(6,*) (a(i),i=0,N-1)
      end
ccc
      subroutine getprm( LV, Bc, mX, N, L, Tp )
      parameter(NX=200, mXX=100)
      integer            LV, Bc, mX, N, L
      real*8                               Tp
      write(6,*)    'Debuf Level 0:none, 1:norm, 2:monitor 3:dump'
       read(5,*) LV
c     write(6,*)    'Boundary Condition: 0:Free 1:Ferro 2:Twist '
c      read(5,*) Bc
                 Bc = 0
      write(6,*)    'Maximum m, less than ', mXX
       read(5,*) mX
      write(6,*)    'Iteration Number L'
       read(5,*) L
      write(6,*)    'System size N, less than ', NX
       read(5,*) N
      write(6,*)    'Temperature T'
       read(5,*) Tp
      write(6,*)    'LV, mX, N, L, Bc, Tp = ',LV,mX,N,L,Bc,Tp
      end
ccc
      subroutine initpt( i, j, N, m,     Fc )
      parameter(iq=2)
      integer            i, j, N, m(0:*)
      real*8                             Fc(0:*)
              i = 1
              j = N
            m(i)= iq
            m(j)= iq
           Fc(i)= 1.0d0
           Fc(j)= 1.0d0
      end
ccc
      subroutine inittm(LV,Bc,T,s,     w,     TL,     TR    )
      parameter(iq=2)
      integer   i,j,k,l,LV,Bc
      real*8          a, x,   T,s(0:*),w(0:*),TL(0:*),TR(0:*)
      real*8                                       Y(0:1,0:1)
c
      if(  T .lt. 100.0d0 ) then
            x = exp(2.0d0/T) +  sqrt( exp(4.0d0/T)-1.0d0 )
            a =   ( 2.0d0*(x+1.0d0/x) )**(1.0d0/4.0d0)
        else
            write(6,*) 'T = ', T,' is replaced by T = inf.'
            x =      1.0d0
            a = sqrt(2.0d0)
        end if
c
        Y(0,0)=      sqrt(x) / a
        Y(1,1)=      sqrt(x) / a
        Y(0,1)=  1 /(sqrt(x) * a)
        Y(1,0)=  1 /(sqrt(x) * a)
      do 40 i = 0, 1
      do 30 j = 0, 1
      do 20 k = 0, 1
      if( Bc.eq.0) then
          TL(4*i+2*j+k)= Y(i,0)*Y(j,0)*Y(k,0) + Y(i,1)*Y(j,1)*Y(k,1)
          TR(4*i+2*j+k)= Y(i,0)*Y(j,0)*Y(k,0) + Y(i,1)*Y(j,1)*Y(k,1)
        end if
      if( Bc.eq.1) then
          TL(4*i+2*j+k)= Y(i,0)*Y(j,0)*Y(k,0)
          TR(4*i+2*j+k)= Y(i,0)*Y(j,0)*Y(k,0)
        end if
      if( Bc.eq.2) then
          TL(4*i+2*j+k)= Y(i,0)*Y(j,0)*Y(k,0)
          TR(4*i+2*j+k)=                        Y(i,1)*Y(j,1)*Y(k,1)
        end if
      do 10 l = 0, 1
       s(i*8+j*4+k*2+l)= Y(i,0)*Y(j,0)*Y(k,0)*Y(l,0)
     &                 - Y(i,1)*Y(j,1)*Y(k,1)*Y(l,1)
       w(i*8+j*4+k*2+l)= Y(i,0)*Y(j,0)*Y(k,0)*Y(l,0)
     &                 + Y(i,1)*Y(j,1)*Y(k,1)*Y(l,1)
   10 continue
   20 continue
   30 continue
   40 continue
      if( LV .ge. 2 ) call dump( '-ww-', iq**4, w  )
      if( LV .ge. 2 ) call dump( '-TL-', iq**3, TL )
      if( LV .ge. 2 ) call dump( '-TR-', iq**3, TR )
      end
ccc
      subroutine pdisp( i, j, N, Fc,      Zp )
      integer           i, j, N, k
      real*8                  F, Fc(0:*), Zp
            F =   -log( -Zp )
      do 10 k = 1, i
            F = F -log(Fc(k))
   10 continue
      do 20 k = j, N
            F = F -log(Fc(k))
   20 continue
       write(6,*) i, j, F
      end
ccc
      subroutine diag( LV, N, E,      Wk,      A    )
      integer    ierr, LV, N, i, j
      real*8                  E(0:*), Wk(0:*), A(0:*)
      call dcsmaa( A, N, N, E, Wk, ierr )
      if( ierr .NE. 0 ) then
          write(6,*) '--- Error Number ',IERR,' in dcsmaa ---'
          stop
        end if
      if( LV .ge. 2 ) call dump( '-dm-', N, E )
      do 20 i = 0, N-1
      do 10 j = 0, N-1
      Wk(i*N+j)= A((N-i-1)*N+j)
   10 continue
   20 continue
      call matcpy( N*N, Wk, A )
      end
ccc
      subroutine extl( mX,m,mN, T,TN,w,A,Fc,d )
      parameter(iq=2,mXX=100)
      integer          mX,m,mN
      real*8    T(0:*),TN(0:*),w(0:*),A(0:*),Fc,d(0:iq**3*mXX**2,0:*)
c
      call  joint( iq,iq*iq*iq,   m*m, w, T, d(0,0)               )
      call exchg3(    iq,iq*iq,   m*m,       d(0,0),d(0,1)        )
      call exchg3(       iq,iq*iq*m,m,              d(0,1),d(0,0) )
c
      mN =  min( iq*m, mX )
      call reno( m, mN, A, TN, Fc, d(0,0), d(0,1), d(0,2) )
      end
ccc
      subroutine extr( mX,m,mN, T,TN, w, A, Fc, d )
      parameter(iq=2,mXX=100)
      integer          mX,m,mN
      real*8    T(0:*),TN(0:*),w(0:*),A(0:*),Fc,d(0:iq**3*mXX**2,0:*)
c
      call exchg2(    iq*iq*iq,iq,     w, d(0,0)                 )
      call  joint( iq,iq*iq*iq,   m*m,    d(0,0),T,d(0,1)        )
      call exchg4(    iq*iq,iq,   m,m,             d(0,1),d(0,0) )
c
      mN =  min( iq*m, mX )
      call reno( m, mN, A, TN, Fc, d(0,0), d(0,1), d(0,2) )
      end
ccc
      subroutine reno( m,mN,A, TN, Fc,  d0,      d1,      d2    )
      parameter(iq=2)
      integer          m,mN
      real*8       Fc, A(0:*), TN(0:*), d0(0:*), d1(0:*), d2(0:*)
c
      call exchg2(         mN,iq*m,       A, d2         )
      call exchg2( iq,iq*m   *iq*m,          d0, d1     )
      call  joint(    iq*m,mN,iq*m   *iq,    d2, d1, d0 )
      call exchg2(         mN,iq*m   *iq,    d0, d1     )
      call  joint(            iq*m,mN,iq*mN, d2, d1, d0 )
      call exchg2(                 mN,iq*mN, d0, TN     )
      call   norm(                    iq*mN*mN,  TN, Fc )
      end
ccc
      subroutine norm( L, V,      W )
      integer       i, L
      real*8              V(0:*), W
      call fmax( L, V, W )
      do 10 i = 0, L-1
          V(i)= V(i) / W
   10 continue
      end
ccc
      subroutine fmax( L, V,      W )
      integer       i, 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
      end
ccc
      subroutine mult( LV, id, is, mL, mR, W, TL, TR, V )
      parameter(iq=2,mXX=100)
      integer          LV, id, is, mL, mR
      real*8       W(0:*),TL(0:*),TR(0:*),  V(0:iq**2*mXX**2,0:*)
      real*8           A0(0:iq**3*mXX**2), A1(0:iq**3*mXX**2)
c
c     ***** Caution: In this program, (-1)*TM is applied to *****
c     *****          the state vectors.                     *****
c
      call exchg2( iq,mL*mL,                   TL, A0             )
      call  joint(    mL,mL*iq,iq*iq*mR,           A0,V(0,is), A1 )
      call exchg2(       mL,iq*iq*iq*mR,       A1, A0             )
      call  joint(    iq*iq,iq*iq,iq*mR   *mL,  W, A0, A1         )
      call exchg2(          iq,iq*iq*mR   *mL,         A1, A0     )
      call  joint(    iq*iq,   iq*iq,mR   *mL*iq,       W, A0, A1 )
      call exchg2(             iq,iq*mR   *mL*iq,      A1, A0     )
      call  joint(                iq*mR,mR,mL*iq*iq,   TR, A0, A1 )
      call exchg2(                      mR,mL*iq*iq,   A1,V(0,id) )
      call chsign(                         mL*iq*iq*mR,   V(0,id) )
      end
ccc
      subroutine exchg2( MH, ML, A,      B    )
      integer      i, j, MH, ML
      real*8                     A(0:*), B(0:*)
      do 20 i = 0, MH-1
      do 10 j = 0, ML-1
      B(j*MH+i) = A(i*ML+j)
   10 continue
   20 continue
      end
ccc
      subroutine exchg3(  MG, MH, ML, A,      B    )
      integer    i, j, k, MG, MH, ML
      real*8                          A(0:*), B(0:*)
      do 30 i = 0, MG-1
      do 20 j = 0, MH-1
      do 10 k = 0, ML-1
      B((j*MG+i)*ML+k) = A((i*MH+j)*ML+k)
   10 continue
   20 continue
   30 continue
      end
ccc
      subroutine exchg4(  MG, MH, MM, ML, A,      B    )
      integer i, j, k, L, MG, MH, MM, ML
      real*8                              A(0:*), B(0:*)
      do 40 i = 0, MG-1
      do 30 j = 0, MH-1
      do 20 k = 0, MM-1
      do 10 L = 0, ML-1
      B(((i*MM+k)*MH+j)*ML+L) = A(((i*MH+j)*MM+k)*ML+L)
   10 continue
   20 continue
   30 continue
   40 continue
      end
ccc
      subroutine joint( MG, MH, ML, A,      B,      C    )
      integer  i, j, k, MG, MH, ML
      real*8                        A(0:*), B(0:*), C(0:*)
      call   vec0( MH*ML, C )
      do 30 i = 0, MG-1
      do 20 j = 0, MH-1
      do 10 k = 0, ML-1
      C(j*ML+k)= C(j*ML+K) + A(i*MH+j)*B(i*ML+k)
   10 continue
   20 continue
   30 continue
      end
ccc
      subroutine ljoint( MH, ML, A,      B   )
      integer   i, j, k, MH, ML
      real*8                     A(0:*), B(0:*)
      call   vec0( MH*MH, B )
      do 30 i = 0, MH-1
      do 20 j = 0, MH-1
      do 10 k = 0, ML-1
      B(i*MH+j)= B(i*MH+j) + A(i*ML+k)*A(j*ML+k)
   10 continue
   20 continue
   30 continue
      end
ccc
      subroutine hjoint( MH, ML, A,      B   )
      integer   i, j, k, MH, ML
      real*8                     A(0:*), B(0:*)
      call   vec0( ML*ML, B )
      do 30 i = 0, MH-1
      do 20 j = 0, ML-1
      do 10 k = 0, ML-1
      B(j*ML+k)= B(j*ML+k) + A(i*ML+j)*A(i*ML+k)
   10 continue
   20 continue
   30 continue
      end
ccc
      subroutine chsign( N, V )
      integer            N
      real*8                V(0:*)
      do 10 i = 0, N-1
          V(i)=-V(i)
   10 continue
      end
ccc
      subroutine matcpy( N, A,      B    )
      integer         i, N
      real*8                A(0:*), B(0:*)
      do 10 i = 0, N-1
          B(i)= A(i)
   10 continue
      end
ccc
      subroutine psiini( NV, mL, mR, V )
      parameter(iq=2)
      integer            NV, mL, mR
      real*8                         V(0:*)
           NV = mL*iq*iq*mR
      do 10 I = 0,NV-1
          V(I)=( dble(irand()) + 00.1 ) / sqrt( dble(NV) )
   10 continue
      end
ccc
cc-------------------------------------------------
c Vector Processors For the Lanczos Diagonalization
c -------------------------------------------------
      subroutine inpro( NV, ID, IS, PRD, V)
      parameter(iq=2,mXX=100)
      integer           NV, ID, IS
      real*8                        PRD, V(0:iq**2*mXX**2,0:*)
         PRD = 0.0d0
      do 1 I = 0, NV-1
         PRD = PRD + V(I,ID)*V(I,IS)
    1 continue
      end
ccc
      subroutine adder( NV, ID, IS1, IS2, A, B, V)
      parameter(iq=2,mXX=100)
      integer           NV, ID, IS1, IS2
      real*8                              A, B, V(0:iq**2*mXX**2,0:*)
      do 1 I  = 0, NV-1
       V(I,ID)= A*V(I,IS1) + B*V(I,IS2)
    1 continue
      end
ccc
      subroutine diver( NV, ID, IS, D, V)
      parameter(iq=2,mXX=100)
      integer           NV, ID, IS
      real*8                        D, V(0:iq**2*mXX**2,0:*)
      do 1 I  = 0, NV-1
       V(I,ID)= V(I,IS) / D
    1 continue
      end
ccc
      subroutine vec0( NV, V)
      integer          NV
      real*8               V(0:*)
      do 1 J = 0, NV-1
         V(J)= 0.0d0
    1 continue
      end
ccc
      subroutine ivec0( NV,IV)
      integer           NV
      integer              IV(0:*)
      do 1 J = 0, NV-1
        IV(J)= 0
    1 continue
      end
ccc
      subroutine vmove( NV, ID, IS, V )
      parameter(iq=2,mXX=100)
      integer           NV, ID, IS
      real*8                        V(0:iq**2*mXX**2,0:*)
      do 1 I  = 0, NV-1
       V(I,ID)= V(I,IS)
    1 continue
      end
ccc            --- copy  to  A from  B ---
      subroutine veccpy( NV, A,      B )
      integer            NV
      real*8                 A(0:*), B(0:*)
      do 1 I = 0, NV-1
         A(I)= B(I)
    1 continue
      end
ccc
      subroutine vecini( NV, V)
      parameter(iq=2,mXX=100)
      integer            NV
      real*8            PRD, V(0:iq**2*mXX**2,0:*)
      call  vec0( NV, V(0,1) )
      call  vec0( NV, V(0,2) )
      call  vec0( NV, V(0,3) )
      call inpro( NV, 0, 0, PRD, V )
          PRD = sqrt(PRD)
      call diver( NV, 0, 0, PRD, V )
      end
ccc
      subroutine lanzos( LV, NV, mL, mR, W, TL, TR, ENG, V )
      parameter( mil = 999 )
      integer            LV, NV, mL, mR
      integer      MLS,  LS, LK
      real*8                 W(0:*), TL(0:*),   TR(0:*)
      real*8   EPS, ENG, ALP(0:mil), TE(0:mil),  V(0:*)
      real*8    E2, E2P, BET(0:mil), ME(0:mil,0:4)
c
          MLS = 190
          EPS = 1.0E-10
           E2 = 0.0d0
      call vecini( NV,          V )
      call  vmove( NV,  3, 0,   V )
      call   mult( LV,  1, 0, mL, mR, W, TL, TR, V )
           LK = 0
           LS = 0
   10 continue
          E2P = E2
      do 20 J = 0, 4
      call macro1( LV,NV,mL,mR,W,TL,TR, V,TE,ALP,BET,LK,0 )
   20 continue
      call malanz( LV, LK, ALP, BET, ENG, E2, ME )
      LS = LS + 1
      IF(  LS .LE. MLS  .AND.  ABS(E2P-E2) .GE. EPS ) goto 10
      LS = LS - 1
      if( LV .ge. 2 ) write(6,100) ENG
  100 format(1H ,' Energy Upper Bound = ',E30.20)
c
      call veccpy( mil+1, TE, ME )
c
      call vmove( NV, 0, 3,                        V )
      call adder( NV, 3, 3,  0,  TE(0), 0.0d0,     V )
      call  mult( LV, 1, 0, mL, mR, W, TL, TR, V )
            LK = 0
      do  50 I = 0, LS
      do  40 J = 0,  4
      call macro1( LV,NV,mL,mR,W,TL,TR, V,TE,ALP,BET,LK,1 )
   40 continue
   50 continue
      end
ccc
      subroutine macro1(LV,NV,mL,mR,W,TL,TR, V,TE,ALP,BET,LK,ISW)
      integer           LV,NV,mL,mR,                      LK,ISW
      real*8                         W(0:*), TL(0:*),  TR(0:*)
      real*8       PROD, BET(0:*), ALP(0:*), TE(0:*),   V(0:*)
      call inpro( NV, 0, 1,           PROD,   V )
      ALP(LK) = PROD
         PROD = -1.0 * PROD
      call adder( NV, 2, 1, 0, 1.0d0, PROD,   V )
      call inpro( NV, 2, 2,           PROD,   V )
         PROD = sqrt(PROD)
      BET(LK) =      PROD
      call diver( NV, 2, 2,           PROD,   V )
      call  mult( LV, 1, 2, mL, mR, W, TL, TR, V )
         PROD = -1.0 * PROD
      call adder( NV, 1, 1, 0, 1.0d0, PROD,   V )
      call vmove( NV, 0, 2,                   V )
           LK = LK + 1
      if( ISW .eq. 1 ) call adder( NV, 3, 3, 0, 1.0d0, TE(LK), V)
      end
ccc
      subroutine malanz( LV, LK, ALP, BET, ENG,  E2,  VE )
      parameter( mil = 999 )
      integer LV,  N,    LK,   LNV, IW1(0:mil), IERR
      real*8  ENG,E2, ALP(0:*), VE(0:mil,0:4)
      real*8  EPS,    BET(0:*),  E(0:mil),   W1(0:6*mil)
        N = LK
      EPS = -1.0d0
      LNV = MIL+1
      call dcstss( ALP,N,BET,EPS,E,5,VE,LNV,-1,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)
      IF( LV .ge. 2 ) write(6,200) (E(J),J = 0,2), LK
  200 format(1H ,'*',3F22.15,'* ',I3)
      end
ccc                      
      subroutine INVERS( LV, NV, mL, mR, W, TL, TR, LMD, V )
      integer    MLR,LR, LV, NV, mL, mR
      real*8  DTP,EPS,DTR,LMD, W(0:*), TL(0:*), TR(0:*), V(0:*)
       LR =  0
      MLR = 40
      EPS = 0.5E-11
      call vmove( NV,  0, 3,                    V )
      GOTO 20
   10 continue
      call cgmeth( LV, NV, mL, mR, W, TL, TR, LMD,V )
   20 call  vmove( NV, 1, 0,                    V )
      call  inpro( NV, 1, 1, DTP,               V )
      DTP =  sqrt(DTP)
      call  diver( NV, 1, 1, DTP,               V )
      call   mult( LV, 2, 1, mL, mR, W, TL, TR, V )
      call  inpro( NV, 2, 2, DTP,               V )
      DTP =  sqrt(DTP)
      call  diver( NV, 2, 2, DTP,               V )
      IF( LMD .lt. 0 ) THEN
          call adder( NV, 3, 2, 1, 1.0d0,1.0d0, V )
        ELSE
          call adder( NV, 3, 2, 1,-1.0d0,1.0d0, V )
        END IF
      call  inpro( NV, 3, 3, DTR,               V )
      if( LV .ge. 2 ) write(6,100) LMD, DTP, sqrt(DTR)
  100 format(1H ,3E25.16 )
c
c     ***** Obtain Minimum Eigenvalue (Hamiltonian/(-1)*TM) *****
c
      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.'
        end if
      LR = LR + 1
      if( LR .ge. MLR ) then
          write(6,*) 'Eigen vector is not obtained in', MLR,' Steps.'
          call errstp( 'invers' )
        end if
      if( sqrt(DTR) .ge. EPS ) goto 10
      end
CCC
      subroutine cgmeth( LV, NV, mL, mR, W, TL, TR, LMD,V )
      parameter(iq=2,mXX=100)
      integer    MLR,LR, LV, NV, mL, mR
      real*8        LMD,  W(0:*), TL(0:*),TR(0:*),V(0:iq**2*mXX**2,0:*)
      real*8  DB,DR1,DR2,ALP,BET,EPS,DTP
       LR = 0
      MLR = 40
      EPS =1.0E-8
C (I)
      call inpro( NV, 1, 1, DB,            V )
      call vmove( NV, 2, 1,                V )
      call  vec0( NV,                      V )
C (II)
      call inpro( NV, 1, 1, DR1,           V )
   20 call  mult( LV, 3, 2, mL, mR, W, TL, TR, V )
      do 30 I = 0, NV-1
      V(I,3)= V(I,3) - LMD*V(I,2)
   30 continue
      call inpro( NV, 3, 2,DTP,            V )
      ALP=DR1/DTP
      call adder( NV, 0, 0, 2, 1.0d0, ALP, V )
      call adder( NV, 1, 1, 3, 1.0d0,-ALP, V )
      call inpro( NV, 1, 1,      DR2,      V )
C
C     write(6,100) DR2 / DB, LR
C 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, V)
      if( LR .ge. MLR ) then
          goto 50
        end if
      goto 20
   40 continue
   50 continue
C     write(6,*) LR,' Rotations in CG'
      end