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