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