PROGRAM A22cgGS3 C CCC for at_22 potential N=66 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C W.Quapp 30.09.2004 for RGF=Newton Trajectory C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886 c test for script version - modular structure c GS cycle for RP: NewPoint calculation by projected gradient c goes along reduced gradient to chain point c c corrector step by CG+ Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Integer istatus,ITALL,itGS,Natoms(22),N3 Integer LA,NDIM,NDIM1,L PARAMETER(NDIM=66,L=25,NDIM1=67) DOUBLE PRECISION XX(NDIM1),G(NDIM1),DDOT,XTOL double precision d(ndim1),gold(ndim1),w(ndim1) double precision Fun,EPS,TLEV,X(NDIM),Proj(NDIM) c double precision time1,time2,tottime REAL*8 Y(L+1,NDIM),VV(NDIM),D1,Yold(NDIM) REAL*8 Rs(NDIM),energy,Cx,gnorm,Pena cc REAL*8 RotX(NDIM),RotY(NDIM),RotZ(NDIM) cc REAL*8 PP(NDIM,NDIM),RR(NDIM,NDIM), REAL*8 ProMat(NDIM,NDIM) logical finish integer iprint(2),iflag,icall,n,method,iter,nfun,IREST Integer JJ,I,J,Np,Np1,IterMax Character*1 CharAtoms(22) cc COMMON /RUNINF/ITER,NFUN Data CharAtoms/'h','c','h','c','h','o','n','c','h','c','h','c', 1 'o','n','c','h','h','h','h','h','h','h'/ c C the actual point, to start the correction, is Y(LA,..) c OPEN(9,FILE='at22chain.dat') OPEN(10,FILE='at22Rs.dat') OPEN(11,FILE='at22param.dat') OPEN(12,FILE='at22proj.dat') OPEN(13,FILE='at22point.dat') OPEN(14,FILE='at22ener.dat') OPEN(15,FILE='at22grad.dat') OPEN(16,FILE='at22engs.weg') OPEN(17,FILE='at22lamb.dat') OPEN(19,FILE='at22grno.dat') OPEN(21,FILE='iflag.dat') OPEN(22,FILE='cgplFi.dat') OPEN(23,FILE='cggsW.dat') OPEN(25,FILE='at22dddd.dat') OPEN(26,FILE='at22gold.dat') OPEN(28,FILE='at22yold.dat') OPEN(38,FILE='at22chain.mol',status='unknown',access='append') OPEN(44,FILE='protocol.txt',status='unknown',access='append') c C the actual point, to start the correction, is Y(LA,..) c Read problem input information c Pena= 1.00001d+0 c N=66 Np = N Np1 = Np+1 N3 = N/3 C Np1 for Lagrange formulation of one additional constraint IPRINT(1)= 0 IPRINT(2)= 0 XTOL= 11.11111D-11 ITERMAX= N + Np1 c Note: EPS1= 1.0D-5 for pure CG+ is very to small here !!! C CG+ method = 3 IREST = 1 rewind 21 read(21,*) IFLAG c write(44,*) ' start cgGS - iflag is ', IFLAG write( 6,*) ' iflag is ', IFLAG if(IFLAG.ne.0) then rewind 22 read(22,*) FINISH,ICALL,NFUN rewind 23 read(23,50) (W(J),J=1,N) rewind 25 read(25,50) (D(J),J=1,Np1) rewind 26 read(26,50) (Gold(J),J=1,Np1) rewind 28 read(28,50) (Yold(J),J=1,N) rewind 17 read(17,50) XX(Np1) endif istatus=0 rewind 11 read(11,*) JJ, LA, eps ,ITALL, itGS cc itGS=itGS+1 -->> other count: by NFUN cc JJ=L !! rewind 11 write(11,*) L, LA, eps ,ITALL, itGS write(44,*)' L=',L,' LA= ',LA,' ITall = ',ITALL,'+',itGS c c Print parameters if (iprint(1).ge.0 .AND. LA.lt.3 .AND. IFLAG.eq.0)then write (*,820) write (*,840) Np1, method, irest write (44,820) write (44,840) Np1, method, irest end if cc 20 CONTINUE If(itGS.gt. N+N+1) then write(6,*) ' Cancel Loop: ERROR , no Convergence at ' write(6,*) ' LA= ' ,LA, ' itGS = ', itGS write(44,*)' Cancel Loop: ERROR , no Convergence at ' write(44,*)' LA= ' ,LA, ' itGS = ', itGS istatus=1 iflag=0 goto 111 Endif c Call the function and gradient values here rewind 10 READ(10,50)(Rs(I),I=1,N) rewind 12 Do 3 J=1,N 3 READ(12,50)(ProMat(J,I),I=1,N) rewind 13 READ(13,50)(X(I),I=1,N) rewind 14 READ(14,*) energy WRITE(44,*)'energy=',energy rewind 9 Do 40, J=1,L+1 40 Read(9,50) (Y(J,I),I=1,N) rewind 15 READ(15,50) (VV(I),I=1,N) c c Do 1, J=1,N c 1 READ(15,*,END=111) J1, JJ, (-VV( 3*(J-1)+I ),I=1,3) c call PERROR("I/O processing") c write(6,*) ' for gradient ' c IF(J1 .ne. N3 ) WRITE(44,*) ' ERROR in gradient reading !!! ' c using the special file with grad output c from Gaussian03 !! note: - sign from gaussian gradient for VV !! c if(IFLAG.eq.0) then Do 211 I=1,N 211 Yold(I)= Y(LA,I) rewind 28 write(28,50) (Yold(I),I=1,N) endif c norm of gradient gnorm= DSQRT(DDOT(N,VV,1,VV,1)) WRITE(6,*) ' PES full gradient norm ',gnorm WRITE(44,*) ' PES full gradient norm ',gnorm c cccccccccccccccccccccccccccccccccccccccc C Newton trajectory with RGF Projector c and rotations out C Projektion of Gradient cccccccccccccccccccccccccccccccccccccccc c Do 13 J=1,N Proj(J)=VV(J) 13 XX(J)=X(J) c cc if(iflag.eq.0) : no cc enforce before every step the "true" lambda XX(Np1)= gnorm/Pena c estimation for lambda of Lagrange ansatz c c turn direction of RS into direction of gradient: If(DDOT(N,Rs,1,VV,1) .lt. 0.0d0) then Do 14 J=1,N 14 Rs(J)= -Rs(J) endif c c Projection of Gradient orthogonal to Rotations and RS call PROJECTION(N,Proj,ProMat) D1= DSQRT(DDOT(N,Proj,1,Proj,1)) WRITE(6,*) ' projected grad norm ', D1 WRITE(44,*)' projected grad norm ', D1 WRITE(44,*)' lambda guess ', XX(Np1) cc WRITE(44,*)' Penalty factor ', Pena If(D1 .lt. eps) THEN WRITE(44,*)'at IteGS=',itGS,' ok - proj norm less then', 1 ' tolerance ', eps WRITE(6, *)' IteGS=',itGS,' == ok - proj norm ',D1 WRITE(44,*)' successful convergence (no ERROR) ' itGS=itGS+NFUN NFUN=0 IFLAG=0 istatus=1 GOTO 111 ENDIF c c CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ c program for running the conjugate gradient methods described in c Gilbert, J.C. and Nocedal, J. (1992). c "Global Convergence Properties of Conjugate Gradient Methods" c SIAM Journal on Optimization, Vol.2, pp. 21-42. c Written by G. Liu, J. Nocedal and R. Waltz October 1998 c IFLAG=0 indicates an initial entry to program c wq: C one constraint: C(x)=0= Sum_i^n rs(i) * (x(i)-Yold(i) ) C Lagrange function: F(x,Lambda)=fcn(x)-Lambda*c(x) C where Lambda = xx(n+1) is the Lagrange parameter c Pena is an additional penalty factor cc Cx=0.0d0 Do 23 J=1,Np 23 Cx= Cx + (X(J)-Yold(J))*Rs(J)*Pena write(44,*) ' C(x) = ' , Cx Fun= energy - XX(Np1)*Cx c c WRITE(44,*)'Fun(x)=',Fun c Do 18 J=1,Np 18 G(J)= VV(J) -XX(Np1)*Rs(J)*Pena G(Np1)= -Cx c c TESTs c write(44,*) ' G(J) = ' , (G(J),J=1,Np1) c D1=0.0d0 c Do 1388 J=1,Np c1388 D1= D1+ G(J)*Rs(J) c write(44,*) 'sum G(J)*Rs(J) = ' ,D1 c D1=0.0d0 c Do 1381 J=1,Np c1381 D1= D1+ VV(J)*Rs(J) c write(44,*) 'sum VV(J)*Rs(J)= ' ,D1 c 30 CONTINUE c Call the main optimization code ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc CALL CGFAM(Np1,XX,Fun,G,D,GOLD,IPRINT,W, * IFLAG,IREST,METHOD,FINISH,ProMat,Rs) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c IFLAG= c 0 : successful termination c 1 : return to evaluate F and G c 2 : return with a new iterate, try termination test cccc -i : error write(44,*) ' NFUN= ', NFUN rewind 11 write(11,*) L, LA, eps ,ITALL, itGS write(44,*) ' at IFLAG',IFLAG,': lambda is after CG+ ',XX(Np1) Do 201,J=1,Np 201 X(J) = XX(J) write(44,*) c actual point approximation: rewind 13 WRITE(13,50) (XX(I),I=1,N) rewind 17 WRITE(17,50) XX(Np1) ccc for gaussian input file c rewind 32 c READ(32,*) (Natoms(I),I=1,N3) c rewind 33 c Do 7, J=1,N3 c 7 WRITE(33,*) Natoms(J), (X(3*(J-1)+I ),I=1,3) c WRITE(33,*) ' ' ccc for gaussian input file rewind 23 write(23,50) (W(J),J=1,Np1) rewind 25 write(25,50) (D(J),J=1,Np1) rewind 26 write(26,50) (Gold(J),J=1,Np1) rewind 21 write(21,*) IFLAG cc write(44,*) ' iflag nach CG+ ', IFLAG write( 6,*) ' iflag nach CG+ ', IFLAG rewind 22 write(22,*) FINISH,ICALL,NFUN IF(IFLAG.LE.0.OR.ICALL.GT.ITERMAX) GO TO 501 IF(IFLAG.EQ.1) THEN ICALL=ICALL + 1 WRITE(6,*) ' Line search LOOP, ICALL= ',ICALL rewind 22 write(22,*) FINISH,ICALL,NFUN rewind 11 write(11,*) L, LA, eps ,ITALL, itGS istatus=0 call exit(istatus) STOP cc GO TO 20 ENDIF IF(IFLAG.EQ.2) THEN c Termination Test. The user may replace it by some other test; c the parameter 'FINISH' must be set to 'TRUE' when the test is satisfied. TLEV= EPS cc wq *(1.0d0 + DABS(F)) I=0 41 I=I+1 IF(I.GT.N) THEN cc IF(I.GT.Np1) THEN cc wq IF(DABS(G(Np1))/Pena.GT.TLEV) go to 30 FINISH = .TRUE. write(44,*) ' finish = true ' itGS=itGS+NFUN NFUN=0 GO TO 30 ENDIF cc wq IF(DABS(G(I)).GT.TLEV) THEN Do 2020,J=1,Np 2020 Proj(J) = G(J) call PROJECTION(Np,Proj,ProMat) IF(DABS(Proj(I)).GT.TLEV) THEN GO TO 30 ELSE GO TO 41 ENDIF ENDIF c 501 continue IF(ICALL.eq.ITERMAX) then istatus=1 write(44,*)'Stop Corr: ICALL=ITERMAX under IFLAG=',IFLAG endif IF(IFLAG.eq.0) then write(44,*)'ICALL= ',ICALL,' under succsessful IFLAG=',IFLAG else write(44,*)'ERROR ? ICALL= ',ICALL,' under IFLAG=',IFLAG endif rewind 11 write(11,*) L, LA, eps ,ITALL, itGS write(44,*)' L=', L,' LA= ', LA ,' ITall = ',ITALL, ' + ',itGS c if (iprint(1).ge.0.and.iflag.ge.0) then write (*,890) Fun write (*,*) ' at LA ', LA end if write(44,890) energy write(44,50) (x(i),i=1,Np) write(44,*) ' C(x)= ', Cx/Pena c write(44,*) 'lambda=?gnorm? = ', XX(Np1)*Pena c * ,' at LA ', LA write(44,*) 111 continue IF(LA.lt.L) write(44,*) ' new step: ' write(38,*) ' 22 ' write(38,*) write(38,505) (x(i),i=1,Np) 505 FORMAT(1X,4H ar ,3F20.10) c cc If(IFLAG.ne.1) then c rewind 32 c READ(32,*) (Natoms(I),I=1,N3) cc Do 17, J=1,N3 c 17 WRITE(44,*) CharAtoms(J),' ',(X(3*(J-1)+I ),I=1,3) cc 17 WRITE(44,50) (X(3*(J-1)+I ),I=1,3) cc endif ccc IF(IFLAG.eq.0) then istatus=1 endif IF(IFLAG.lt.0) then write(6,*)' Cancel Loop: no Convergence in CG+ ' c write(44,*)'ERROR : Corr no Convergence in CG+ ' rewind 11 write(11,*) L, LA, eps ,ITALL, itGS write(44,*)' L=', L,' LA= ', LA ,' ITall = ',ITALL, ' + ',itGS istatus=1 call exit(istatus) STOP endif c rewind 21 write(21,*) iflag c Istatus controls the script loops If(LA.gt.L) istatus=10 c new chain point becomes: IF(finish) then Do 117 I=1,Np 117 Y(LA,I)= XX(I) rewind 9 Do 118, J=1,L+1 118 write(9,50) (Y(J,I),I=1,N) WRITE(6,*)' energy= ',energy WRITE(16,*) '{ ' IF(LA.le.L) WRITE(16,*) ' ', energy,' , ' IF(LA.eq.L+1)WRITE(16,*) ' ', energy,' } ' endif rewind 11 write(11,*) L, LA, eps ,ITALL, itGS write(44,*)' at the end of frame of CGFAM ' write(44,*)' L=', L,' LA= ', LA ,' ITall = ',ITALL, ' + ',itGS call exit(istatus) STOP c 50 FORMAT(1X,3F20.10) 820 format ( ' Conjugate Gradient Minimization Routine') 840 format ( ' n =', i6, * ' method =', i6, * ' irest =', i6) 850 format (/,' Error: negative N value') 860 format (/,' Error: N too large, increase parameter ndim') 890 format (/,' f(x*) =', 1pd16.8) c END cccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE PROJECTION(N,VEC,ProMat) c projection of vector VEC of dim N c by Projection Matrix ProMat of dim(N,N) PARAMETER(Ndim=66) REAL*8 VEC(N),ProMat(N,N),Proj(Ndim) INTEGER N,J,JJ c Do 1 J=1,N Proj(J)=0.0d0 Do 2 JJ=1,N 2 Proj(J)=Proj(J) + Vec(JJ)*ProMat(J,JJ) 1 Continue Do 3 J=1,N 3 Vec(J) = Proj(J) RETURN END C Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C C SUBROUTINE CVSMOD(N,X,F,G,S,STP,FTOL,GTOL,XTOL, * STPMIN,STPMAX,MAXFEV,INFO,NFEV,WA,DGINIT,DGOUT) INTEGER N,MAXFEV,INFO,NFEV DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX DOUBLE PRECISION X(N),G(N),S(N),WA(N) C C *** This is a modification of More's line search routine ** C * * * * * * C THE PURPOSE OF CVSMOD IS TO FIND A STEP WHICH SATISFIES C A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION. C THE USER MUST PROVIDE A SUBROUTINE WHICH CALCULATES THE C FUNCTION AND THE GRADIENT. C C AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF C UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF C UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A C MINIMIZER OF THE MODIFIED FUNCTION C C F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S). C C IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION C HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE, C THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT C CONTAINS A MINIMIZER OF F(X+STP*S). C C THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES C THE SUFFICIENT DECREASE CONDITION C C F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S), C C AND THE CURVATURE CONDITION C C ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S). C C IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION C IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES C BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH C CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING C ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY C SATISFIES THE SUFFICIENT DECREASE CONDITION. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE CVSMOD(N,X,F,G,S,STP,FTOL,GTOL,XTOL, C STPMIN,STPMAX,MAXFEV,INFO,NFEV,WA,DG,DGOUT) C WHERE C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. C X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE C BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS C X + STP*S. C F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F C AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S. C G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE C GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT C OF F AT X + STP*S. C S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE C SEARCH DIRECTION. C STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN C INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT C STP CONTAINS THE FINAL ESTIMATE. C FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION C OCCURS WHEN THE SUFFICIENT DECREASE CONDITION AND THE C DIRECTIONAL DERIVATIVE CONDITION ARE SATISFIED. C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS C WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY C IS AT MOST XTOL. C STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH C SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION C OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST C MAXFEV BY THE END OF AN ITERATION. C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: C C INFO = 0 IMPROPER INPUT PARAMETERS. C INFO =-1 A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT. C INFO = 1 THE SUFFICIENT DECREASE CONDITION AND THE C DIRECTIONAL DERIVATIVE CONDITION HOLD. C INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY C IS AT MOST XTOL. C INFO = 3 NUMBER OF CALLS TO FCN HAS REACHED MAXFEV. C INFO = 4 THE STEP IS AT THE LOWER BOUND STPMIN. C INFO = 5 THE STEP IS AT THE UPPER BOUND STPMAX. C INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS. C THERE MAY NOT BE A STEP WHICH SATISFIES THE C SUFFICIENT DECREASE AND CURVATURE CONDITIONS. C TOLERANCES MAY BE TOO SMALL. C C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN. C WA IS A WORK ARRAY OF LENGTH N. C C *** The following two parameters are a modification to the code C DG IS THE INITIAL DIRECTIONAL DERIVATIVE (IN THE ORIGINAL CODE C IT WAS COMPUTED IN THIS ROUTINE0 C DGOUT IS THE VALUE OF THE DIRECTIONAL DERIVATIVE WHEN THE WOLFE C CONDITIONS HOLD, AND AN EXIT IS MADE TO CHECK DESCENT. C C SUBPROGRAMS CALLED C CSTEPM C FORTRAN-SUPPLIED...ABS,MAX,MIN C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983 C JORGE J. MORE', DAVID J. THUENTE C INTEGER INFOC,J LOGICAL BRACKT,STAGE1 DOUBLE PRECISION DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM, * FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY, * STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,ZERO,DGOUT DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/ save OPEN(44,FILE='protocol.txt',status='unknown',access='append') OPEN(31,FILE='cvsRest.dat') C wq for scratch form IF(INFO.ne.0) then rewind31 read(31,*)DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM, * FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY, * STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,DGOUT, * INFOC, BRACKT,STAGE1 endif C end wq IF(INFO.EQ.-1) GOTO 45 IF(INFO.EQ.1) GOTO 321 INFOC = 1 C wq if(STP .LE. ZERO) then write(44,*) 'ERROR - warning: STP < 0 ??? ' goto 1 endif C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. STP .LE. ZERO .OR. FTOL .LT. ZERO .OR. * GTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. STPMIN .LT. ZERO * .OR. STPMAX .LT. STPMIN .OR. MAXFEV .LE. 0) Then write(44,*) ' return out off CVSMOD by ERROR ' RETURN endif 1 continue C C COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION C AND CHECK THAT S IS A DESCENT DIRECTION. C IF (DGINIT .GE. ZERO) then c rewind31 c write(31,*)DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM, c * FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY, c * STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,DGOUT, c * INFOC, BRACKT,STAGE1 c RETURN write(44,*) 'ERROR - warning: DGINIT.GE.ZERO 0 ??? ' endif C C INITIALIZE LOCAL VARIABLES. C BRACKT = .FALSE. STAGE1 = .TRUE. NFEV = 0 FINIT = F DGTEST = FTOL*DGINIT WIDTH = STPMAX - STPMIN WIDTH1 = WIDTH/P5 DO 20 J = 1, N WA(J) = X(J) 20 CONTINUE C C THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP, C FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP. C THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP, C FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF C THE INTERVAL OF UNCERTAINTY. C THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP, C FUNCTION, AND DERIVATIVE AT THE CURRENT STEP. C STX = ZERO FX = FINIT DGX = DGINIT STY = ZERO FY = FINIT DGY = DGINIT C C START OF ITERATION. C 30 CONTINUE C C SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND C TO THE PRESENT INTERVAL OF UNCERTAINTY. C IF (BRACKT) THEN STMIN = MIN(STX,STY) STMAX = MAX(STX,STY) ELSE STMIN = STX STMAX = STP + XTRAPF*(STP - STX) END IF C C FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN. C STP = MAX(STP,STPMIN) STP = MIN(STP,STPMAX) C C IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET C STP BE THE LOWEST POINT OBTAINED SO FAR. C IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX)) * .OR. NFEV .GE. MAXFEV-1 .OR. INFOC .EQ. 0 * .OR. (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX)) STP = STX C C EVALUATE THE FUNCTION AND GRADIENT AT STP C AND COMPUTE THE DIRECTIONAL DERIVATIVE. C DO 40 J = 1, N X(J) = WA(J) + STP*S(J) 40 CONTINUE C Return to compute function value INFO=-1 rewind31 write(31,*)DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM, * FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY, * STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,DGOUT, * INFOC, BRACKT,STAGE1 RETURN C 45 INFO=0 NFEV = NFEV + 1 DG = ZERO DO 50 J = 1, N DG = DG + G(J)*S(J) 50 CONTINUE FTEST1 = FINIT + STP*DGTEST C C TEST FOR CONVERGENCE. C IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX)) * .OR. INFOC .EQ. 0) INFO = 6 IF (STP .EQ. STPMAX .AND. * F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5 IF (STP .EQ. STPMIN .AND. * (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4 IF (NFEV .GE. MAXFEV) INFO = 3 IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2 C More's code has been modified so that at least one new C function value is computed during the line search (enforcing C at least one interpolation is not easy, since the code may C override an interpolation) IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT) * .AND. NFEV.GT.1) INFO = 1 C C CHECK FOR TERMINATION. C IF (INFO .NE. 0)THEN DGOUT=DG rewind31 write(31,*)DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM, * FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY, * STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,DGOUT, * INFOC, BRACKT,STAGE1 RETURN ENDIF 321 continue C C IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED C FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE. C IF (STAGE1 .AND. F .LE. FTEST1 .AND. * DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE. C C A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF C WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED C FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE C DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN C OBTAINED BUT THE DECREASE IS NOT SUFFICIENT. C IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN C C DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES. C FM = F - STP*DGTEST FXM = FX - STX*DGTEST FYM = FY - STY*DGTEST DGM = DG - DGTEST DGXM = DGX - DGTEST DGYM = DGY - DGTEST C C CALL CSTEPM TO UPDATE THE INTERVAL OF UNCERTAINTY C AND TO COMPUTE THE NEW STEP. C CALL CSTEPM(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM, * BRACKT,STMIN,STMAX,INFOC) C C RESET THE FUNCTION AND GRADIENT VALUES FOR F. C FX = FXM + STX*DGTEST FY = FYM + STY*DGTEST DGX = DGXM + DGTEST DGY = DGYM + DGTEST ELSE C C CALL CSTEPM TO UPDATE THE INTERVAL OF UNCERTAINTY C AND TO COMPUTE THE NEW STEP. C CALL CSTEPM(STX,FX,DGX,STY,FY,DGY,STP,F,DG, * BRACKT,STMIN,STMAX,INFOC) END IF C C FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE C INTERVAL OF UNCERTAINTY. C IF (BRACKT) THEN IF (ABS(STY-STX) .GE. P66*WIDTH1) * STP = STX + P5*(STY - STX) WIDTH1 = WIDTH WIDTH = ABS(STY-STX) END IF C C END OF ITERATION. C GO TO 30 C C LAST CARD OF SUBROUTINE CVSMOD. C END C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE CSTEPM(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT, * STPMIN,STPMAX,INFO) INTEGER INFO DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX LOGICAL BRACKT,BOUND C ********** C C SUBROUTINE CSTEPM C C THE PURPOSE OF CSTEPM IS TO COMPUTE A SAFEGUARDED STEP FOR C A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR C A MINIMIZER OF THE FUNCTION. C C THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION C VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS C ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE C DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A C MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY C WITH ENDPOINTS STX AND STY. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE CSTEPM(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT, C STPMIN,STPMAX,INFO) C C WHERE C C STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP, C THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED C SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION C OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE C SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY. C C STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP, C THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF C THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE C UPDATED APPROPRIATELY. C C STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP, C THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP. C IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE C BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP. C C BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER C HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED C THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER C IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE. C C STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER C AND UPPER BOUNDS FOR THE STEP. C C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: C IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED C ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE C INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS. C C FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983 C JORGE J. MORE', DAVID J. THUENTE C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA INFO = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF ((BRACKT .AND. (STP .LE. MIN(STX,STY) .OR. * STP .GE. MAX(STX,STY))) .OR. * DX*(STP-STX) .GE. 0.0 .OR. STPMAX .LT. STPMIN) RETURN C C DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN. C SGND = DP*(DX/ABS(DX)) C C FIRST CASE. A HIGHER FUNCTION VALUE. C THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER C TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN, C ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN. C IF (FP .GT. FX) THEN INFO = 1 BOUND = .TRUE. THETA = 3*(FX - FP)/(STP - STX) + DX + DP S = MAX(ABS(THETA),ABS(DX),ABS(DP)) GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S)) IF (STP .LT. STX) GAMMA = -GAMMA P = (GAMMA - DX) + THETA Q = ((GAMMA - DX) + GAMMA) + DP R = P/Q STPC = STX + R*(STP - STX) STPQ = STX + ((DX/((FX-FP)/(STP-STX)+DX))/2)*(STP - STX) IF (ABS(STPC-STX) .LT. ABS(STPQ-STX)) THEN STPF = STPC ELSE STPF = STPC + (STPQ - STPC)/2 END IF BRACKT = .TRUE. C C SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF C OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC C STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP, C THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN. C ELSE IF (SGND .LT. 0.0) THEN INFO = 2 BOUND = .FALSE. THETA = 3*(FX - FP)/(STP - STX) + DX + DP S = MAX(ABS(THETA),ABS(DX),ABS(DP)) GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S)) IF (STP .GT. STX) GAMMA = -GAMMA P = (GAMMA - DP) + THETA Q = ((GAMMA - DP) + GAMMA) + DX R = P/Q STPC = STP + R*(STX - STP) STPQ = STP + (DP/(DP-DX))*(STX - STP) IF (ABS(STPC-STP) .GT. ABS(STPQ-STP)) THEN STPF = STPC ELSE STPF = STPQ END IF BRACKT = .TRUE. C C THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE C SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES. C THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY C IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC C IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE C EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO C COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP C CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN. C ELSE IF (ABS(DP) .LT. ABS(DX)) THEN INFO = 3 BOUND = .TRUE. THETA = 3*(FX - FP)/(STP - STX) + DX + DP S = MAX(ABS(THETA),ABS(DX),ABS(DP)) C C THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND C TO INFINITY IN THE DIRECTION OF THE STEP. C GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S))) IF (STP .GT. STX) GAMMA = -GAMMA P = (GAMMA - DP) + THETA Q = (GAMMA + (DX - DP)) + GAMMA R = P/Q IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN STPC = STP + R*(STX - STP) ELSE IF (STP .GT. STX) THEN STPC = STPMAX ELSE STPC = STPMIN END IF STPQ = STP + (DP/(DP-DX))*(STX - STP) IF (BRACKT) THEN IF (ABS(STP-STPC) .LT. ABS(STP-STPQ)) THEN STPF = STPC ELSE STPF = STPQ END IF ELSE IF (ABS(STP-STPC) .GT. ABS(STP-STPQ)) THEN STPF = STPC ELSE STPF = STPQ END IF END IF C C FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE C SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES C NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP C IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN. C ELSE INFO = 4 BOUND = .FALSE. IF (BRACKT) THEN THETA = 3*(FP - FY)/(STY - STP) + DY + DP S = MAX(ABS(THETA),ABS(DY),ABS(DP)) GAMMA = S*SQRT((THETA/S)**2 - (DY/S)*(DP/S)) IF (STP .GT. STY) GAMMA = -GAMMA P = (GAMMA - DP) + THETA Q = ((GAMMA - DP) + GAMMA) + DY R = P/Q STPC = STP + R*(STY - STP) STPF = STPC ELSE IF (STP .GT. STX) THEN STPF = STPMAX ELSE STPF = STPMIN END IF END IF C C UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT C DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE. C IF (FP .GT. FX) THEN STY = STP FY = FP DY = DP ELSE IF (SGND .LT. 0.0) THEN STY = STX FY = FX DY = DX END IF STX = STP FX = FP DX = DP END IF C C COMPUTE THE NEW STEP AND SAFEGUARD IT. C STPF = MIN(STPMAX,STPF) STPF = MAX(STPMIN,STPF) STP = STPF IF (BRACKT .AND. BOUND) THEN IF (STY .GT. STX) THEN STP = MIN(STX+0.66*(STY-STX),STP) ELSE STP = MAX(STX+0.66*(STY-STX),STP) END IF END IF RETURN C C LAST CARD OF SUBROUTINE CSTEPM. C END C C -------------------------------------------------------------------- C Conjugate Gradient methods for solving unconstrained nonlinear C optimization problems, as described in the paper: C Gilbert, J.C. and Nocedal, J. (1992). "Global Convergence Properties C of Conjugate Gradient Methods", SIAM Journal on Optimization, Vol. 2, C pp. 21-42. C A web-based Server which solves unconstrained nonlinear optimization C problems using this Conjugate Gradient code can be found at: C http://www-neos.mcs.anl.gov/neos/solvers/UCO:CGPLUS/ C -------------------------------------------------------------------- C SUBROUTINE CGFAM(N,X,F,G,D,GOLD,IPRINT,W, * IFLAG,IREST,METHOD,FINISH,ProMat,Rs) C C Subroutine parameters C INTEGER N,oldN PARAMETER (oldN=66) DOUBLE PRECISION X(N),G(N),D(N),GOLD(N),W(N),F 1 ,Rs(oldN),ProMat(oldN,oldN) INTEGER IPRINT(2),IFLAG,IREST,METHOD,IM,NDES C C N = NUMBER OF VARIABLES C X = ITERATE C F = FUNCTION VALUE C G = GRADIENT VALUE C GOLD = PREVIOUS GRADIENT VALUE C IPRINT = FREQUENCY AND TYPE OF PRINTING C IPRINT(1) < 0 : NO OUTPUT IS GENERATED C IPRINT(1) = 0 : OUTPUT ONLY AT FIRST AND LAST ITERATION C IPRINT(1) > 0 : OUTPUT EVERY IPRINT(1) ITERATIONS C IPRINT(2) : SPECIFIES THE TYPE OF OUTPUT GENERATED; C THE LARGER THE VALUE (BETWEEN 0 AND 3), C THE MORE INFORMATION C IPRINT(2) = 0 : NO ADDITIONAL INFORMATION PRINTED C IPRINT(2) = 1 : INITIAL X AND GRADIENT VECTORS PRINTED C IPRINT(2) = 2 : X VECTOR PRINTED EVERY ITERATION C IPRINT(2) = 3 : X VECTOR AND GRADIENT VECTOR PRINTED C EVERY ITERATION C EPS = CONVERGENCE CONSTANT !! not used ?? here -> deleted C W = WORKING ARRAY C IFLAG = CONTROLS TERMINATION OF CODE, AND RETURN TO MAIN C PROGRAM TO EVALUATE FUNCTION AND GRADIENT C IFLAG = -3 : IMPROPER INPUT PARAMETERS C IFLAG = -2 : DESCENT WAS NOT OBTAINED C IFLAG = -1 : LINE SEARCH FAILURE C IFLAG = 0 : INITIAL ENTRY OR C SUCCESSFUL TERMINATION WITHOUT ERROR C IFLAG = 1 : INDICATES A RE-ENTRY WITH NEW FUNCTION VALUES C IFLAG = 2 : INDICATES A RE-ENTRY WITH A NEW ITERATE C IREST = 0 (NO RESTARTS); 1 (RESTART EVERY N STEPS) C METHOD = 1 : FLETCHER-REEVES C 2 : POLAK-RIBIERE C 3 : POSITIVE POLAK-RIBIERE: BETA=MAX{BETA,0} C Local variables C DOUBLE PRECISION GTOL,ONE,ZERO,GNORM,STP1,FTOL,XTOL,STPMIN, . STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD, . DGOUT,DG,DG1,d1,ddot INTEGER MP,LP,ITER,NFUN,MAXFEV,INFO,I,NFEV,NRST,IDES LOGICAL NEW,FINISH C C THE FOLLOWING PARAMETERS ARE PLACED IN COMMON BLOCKS SO THEY C CAN BE EASILY ACCESSED ANYWHERE IN THE CODE C C MP = UNIT NUMBER WHICH DETERMINES WHERE TO WRITE REGULAR OUTPUT C LP = UNIT NUMBER WHICH DETERMINES WHERE TO WRITE ERROR OUPUT COMMON /CGDD/MP,LP C C ITER: KEEPS TRACK OF THE NUMBER OF ITERATIONS C NFUN: KEEPS TRACK OF THE NUMBER OF FUNCTION/GRADIENT EVALUATIONS COMMON /RUNINF/ITER,NFUN SAVE DATA ONE/1.0D+0/,ZERO/0.0D+0/ OPEN(24,FILE='cgRest.dat') cc OPEN(32,FILE='alaatoms.dat',status='old') cc OPEN(33,FILE='alapoig.dat') IF(IFLAG.eq.1) then rewind 24 read(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN, . STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD, . DGOUT,DG,DG1,d1, . MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES, . NEW,FINISH ,IM,NDES endif C C IFLAG = 1 INDICATES A RE-ENTRY WITH NEW FUNCTION VALUES IF(IFLAG.EQ.1) GOTO 72 C C IFLAG = 2 INDICATES A RE-ENTRY WITH A NEW ITERATE IF(IFLAG.EQ.2) GOTO 80 C C INITIALIZE C ---------- C C IM = NUMBER OF TIMES BETAPR WAS NEGATIVE FOR METHOD 2 OR C NUMBER OF TIMES BETAPR WAS 0 FOR METHOD 3 C C NDES = NUMBER OF LINE SEARCH ITERATIONS AFTER WOLFE CONDITIONS C WERE SATISFIED C ITER= 0 IF(N.LE.0) GO TO 96 NFUN= 1 NEW=.TRUE. NRST= 0 IM=0 NDES=0 C DO 5 I=1,N 5 D(I)= -G(I) cc cc wq: here in projection, use old N c oldN=N-1 call PROJECTION(oldN,D,ProMat) c cc wq c D1= DDOT(66,D,1,Rs,1) c write(44,*) 'projection of D(I) * Rs ', D1 c write(44,*) c wq: N-1 ?? GNORM= DSQRT(DDOT(N,D,1,D,1)) cc write(44,*) 'D-norm-Np1 ', GNORM cc here already shortened length of projected D used !! STP1= ONE/GNORM C C PARAMETERS FOR LINE SEARCH ROUTINE C ---------------------------------- C C FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION C OCCURS WHEN THE SUFFICIENT DECREASE CONDITION AND THE C DIRECTIONAL DERIVATIVE CONDITION ARE SATISFIED. C C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS C WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY C IS AT MOST XTOL. C C STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH C SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. C C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION C OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST C MAXFEV BY THE END OF AN ITERATION. FTOL= 1.0D-4 GTOL= 1.0D-1 IF(GTOL.LE.1.D-04) THEN IF(LP.GT.0) WRITE(LP,145) GTOL=1.D-02 ENDIF XTOL= 1.0D-17 STPMIN= 1.0D-20 STPMAX= 1.0D+20 MAXFEV= 20 cc wq new: stpmin= 0.23081947D-5 stpmax= 0.1200023081947d+0 MAXFEV= N + N cc wq C IF(IPRINT(1).GE.0) CALL CGBD(IPRINT,ITER,NFUN, * GNORM,N,X,F,G,STP,FINISH,NDES,IM,BETAFR,BETAPR,BETA) C C MAIN ITERATION LOOP C --------------------- C cc 88888 cc 8 ITER= ITER+1 C WHEN NRST>N AND IREST=1 THEN RESTART NRST= NRST+1 INFO=0 C C CALL THE LINE SEARCH ROUTINE OF MOR'E AND THUENTE C (modified for our CG method) C ------------------------------------------------- C JJ Mor'e and D Thuente, "Linesearch Algorithms with Guaranteed C Sufficient Decrease". ACM Transactions on Mathematical C Software 20 (1994), pp 286-307. C NFEV=0 DO 70 I=1,N cc W(I)=G(I) 70 GOLD(I)= G(I) cc wq ????? c here, 66 is old N c call PROJECTION(oldN,W,ProMat) c DO 71 I=1,N c 71 GOLD(I)= W(I) cc wq DG= DDOT(N,D,1,G,1) C wq ????? Notbremse !!! c if(DG.gt.0.0d0) DG=-DG c write(44,*) ' ERROR ? Notbremse !!! DG=-DG ' DGOLD=DG STP=ONE C C Shanno-Phua's Formula For Trial Step C IF(.NOT.NEW) STP= DG0/DG IF (ITER.EQ.1) STP=STP1 IDES=0 new=.false. 72 CONTINUE C wq ????? Notbremse !!! if(stp.lt.0.0d0) then stp= 1.0D-3 write(44,*) ' ERROR ? stp is - : changed to 10^-3 ' endif write(6,*) 'step = ', stp write(44,*) 'step = ', stp C wq C C Call to the line search subroutine C CALL CVSMOD(N,X,F,G,D,STP,FTOL,GTOL, * XTOL,STPMIN,STPMAX,MAXFEV,INFO,NFEV,W,DG,DGOUT) C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: C INFO = 0 IMPROPER INPUT PARAMETERS. C INFO =-1 A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT. C INFO = 1 THE SUFFICIENT DECREASE CONDITION AND THE C DIRECTIONAL DERIVATIVE CONDITION HOLD. C INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY C IS AT MOST XTOL. C INFO = 3 NUMBER OF CALLS TO FCN HAS REACHED MAXFEV. C INFO = 4 THE STEP IS AT THE LOWER BOUND STPMIN. C INFO = 5 THE STEP IS AT THE UPPER BOUND STPMAX. C INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS. C THERE MAY NOT BE A STEP WHICH SATISFIES THE C SUFFICIENT DECREASE AND CURVATURE CONDITIONS. C TOLERANCES MAY BE TOO SMALL. cc cc write(44,*) 'step-cvs ', stp cc IF (INFO .EQ. -1) THEN C RETURN TO FETCH FUNCTION AND GRADIENT IFLAG=1 rewind 24 write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN, . STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD, . DGOUT,DG,DG1,d1, . MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES, . NEW,FINISH ,IM,NDES RETURN ENDIF c wq INFO5 IF(INFO.EQ.4 .and. ITER.lt.33) GO TO 91 IF(INFO.EQ.5 .and. ITER.lt.33) GO TO 91 IF(INFO.EQ.6 .and. ITER.lt.33) GO TO 91 IF(INFO.NE.1) GO TO 90 91 continue c wq C C TEST IF DESCENT DIRECTION IS OBTAINED FOR METHODS 2 AND 3 C --------------------------------------------------------- C GG= DDOT(N,G,1,G,1) GG0= DDOT(N,G,1,GOLD,1) BETAPR= (GG-GG0)/GNORM**2 IF (IREST.EQ.1.AND.NRST.GT.N) THEN NRST=0 NEW=.TRUE. GO TO 75 ENDIF C IF (METHOD.EQ.1) THEN GO TO 75 ELSE DG1=-GG + BETAPR*DGOUT IF (DG1.lt. 0.0d0 ) GO TO 75 IF (IPRINT(1).GE.0) write(6,*) 'no descent' IF (IPRINT(1).GE.0) write(44,*) 'ERROR, no descent' IDES= IDES + 1 IF(IDES.GT.5) GO TO 95 GO TO 72 ENDIF C C DETERMINE CORRECT BETA VALUE FOR METHOD CHOSEN C ---------------------------------------------- C C IM = NUMBER OF TIMES BETAPR WAS NEGATIVE FOR METHOD 2 OR C NUMBER OF TIMES BETAPR WAS 0 FOR METHOD 3 C C NDES = NUMBER OF LINE SEARCH ITERATIONS AFTER WOLFE CONDITIONS C WERE SATISFIED C 75 NFUN= NFUN + NFEV NDES= NDES + IDES BETAFR= GG/GNORM**2 IF (NRST.EQ.0) THEN BETA= ZERO ELSE IF (METHOD.EQ.1) BETA=BETAFR IF (METHOD.EQ.2) BETA=BETAPR IF ((METHOD.EQ.2.OR.METHOD.EQ.3).AND.BETAPR.LT.0) IM=IM+1 IF (METHOD.EQ.3) BETA=MAX(ZERO,BETAPR) ENDIF C C COMPUTE THE NEW DIRECTION C -------------------------- C DO 78 I=1,N c W(I) = G(I) 78 D(I) = -G(I) +BETA*D(I) c wq here, in projection, use is old N c oldN=N-1 call PROJECTION(oldN,D,ProMat) c c IF(DDOT(N,G,1,D,1) .gt. 0.0d0) then c d1=-1.0d0 c write(44,*) 'D - turned by ', d1 c else c d1=1.0d0 c endif c DO 778 I=1,N c778 D(I) = d1*D(I) c wq end c DG0= DGOLD*STP ccc xxx IF(DDOT(N,G,1,D,1) .gt. 0.0d0) then C C RETURN TO DRIVER FOR TERMINATION TEST C ------------------------------------- C cc wq ??? N-1 GNORM=DSQRT(DDOT(N,G,1,G,1)) cc c call PROJECTION(66,W,ProMat) c GNORM=DSQRT(DDOT(66,W,1,W,1)) cc IFLAG=2 rewind 24 write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN, . STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD, . DGOUT,DG,DG1,d1, . MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES, . NEW,FINISH ,IM,NDES RETURN 80 CONTINUE C C Call subroutine for printing output C IF(IPRINT(1).GE.0) CALL CGBD(IPRINT,ITER,NFUN, * GNORM,N,X,F,G,STP,FINISH,NDES,IM,BETAFR,BETAPR,BETA) IF (FINISH) THEN IFLAG = 0 rewind 24 write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN, . STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD, . DGOUT,DG,DG1,d1, . MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES, . NEW,FINISH ,IM,NDES RETURN END IF GOTO 8 C C ---------------------------------------- C END OF MAIN ITERATION LOOP. ERROR EXITS. C ---------------------------------------- C 90 IFLAG=-1 IF(LP.GT.0) WRITE(LP,100) INFO rewind 24 write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN, . STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD, . DGOUT,DG,DG1,d1, . MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES, . NEW,FINISH ,IM,NDES RETURN 95 IFLAG=-2 IF(LP.GT.0) WRITE(LP,135) I rewind 24 write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN, . STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD, . DGOUT,DG,DG1,d1, . MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES, . NEW,FINISH ,IM,NDES RETURN 96 IFLAG= -3 IF(LP.GT.0) WRITE(LP,140) rewind 24 write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN, . STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD, . DGOUT,DG,DG1,d1, . MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES, . NEW,FINISH ,IM,NDES RETURN C 100 FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED - ERROR RETURN', . ' OF LINE SEARCH: INFO= ',I2) cc . ,/,' POSSIBLE CAUSE: FUNCTION OR GRADIENT ARE INCORRECT') 135 FORMAT(/' IFLAG= -2',/' DESCENT WAS NOT OBTAINED') 140 FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N', . ' IS NOT POSITIVE)') 145 FORMAT(/' GTOL IS LESS THAN OR EQUAL TO 1.D-04', . / ' IT HAS BEEN RESET TO 1.D-02') END C C LAST LINE OF ROUTINE CGFAM C C************************************************************************** SUBROUTINE CGBD(IPRINT,ITER,NFUN, * GNORM,N,X,F,G,STP,FINISH,NDES,IM,BETAFR,BETAPR,BETA) C C --------------------------------------------------------------------- C THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND AMOUNT C OF OUTPUT ARE CONTROLLED BY IPRINT. C --------------------------------------------------------------------- C INTEGER N DOUBLE PRECISION X(N),G(N),F,GNORM,STP,BETAFR,BETAPR,BETA,QQ INTEGER IPRINT(2),ITER,NFUN,LP,MP,NDES,IM,Iqq,I LOGICAL FINISH COMMON /CGDD/MP,LP C ??: QQ=betafr QQ=betapr Iqq=NDES Iqq= IM + Iqq* QQ QQ = Iqq c wq : ?? employ never used variables IF (ITER.EQ.0)THEN PRINT* WRITE(MP,10) WRITE(MP,20) N WRITE(MP,30) F,GNORM IF (IPRINT(2).GE.1)THEN WRITE(MP,40) WRITE(MP,50) (X(I),I=1,N) WRITE(MP,60) WRITE(MP,50) (G(I),I=1,N) ENDIF WRITE(MP,10) WRITE(MP,70) ELSE IF ((IPRINT(1).EQ.0).AND.(ITER.NE.1.AND..NOT.FINISH))RETURN IF (IPRINT(1).NE.0)THEN IF(MOD(ITER-1,IPRINT(1)).EQ.0.OR.FINISH)THEN IF(IPRINT(2).GT.1.AND.ITER.GT.1) WRITE(MP,70) WRITE(MP,80)ITER,NFUN,F,GNORM,STP,BETA ELSE RETURN ENDIF ELSE IF( IPRINT(2).GT.1.AND.FINISH) WRITE(MP,70) WRITE(MP,80)ITER,NFUN,F,GNORM,STP,BETA ENDIF IF (IPRINT(2).EQ.2.OR.IPRINT(2).EQ.3)THEN WRITE(MP,40) WRITE(MP,50)(X(I),I=1,N) IF (IPRINT(2).EQ.3)THEN WRITE(MP,60) WRITE(MP,50)(G(I),I=1,N) ENDIF ENDIF IF (FINISH) WRITE(MP,100) ENDIF C 10 FORMAT('*************************************************') 20 FORMAT(' N=',I5, / , ' INITIAL VALUES:') 30 FORMAT(' F= ',1PD10.3,' GNORM= ',1PD10.3) 40 FORMAT(/,' VECTOR X= ') 50 FORMAT(6(2X,1PD10.3/)) 60 FORMAT(' GRADIENT VECTOR G= ') 70 FORMAT(/' I NFN',4X,'FUNC',7X,'GNORM',6X, * 'STEPLEN',4x,'BETA',/, * ' ----------------------------------------------------') 80 FORMAT(I4,1X,I3,2X,2(1PD10.3,2X),1PD8.1,2x,1PD8.1) 100 FORMAT(/' SUCCESSFUL CONVERGENCE (NO ERRORS).' * ,/,' IFLAG = 0') C RETURN END C LAST LINE OF CGBD C************************************************************************* C DATA C MP = UNIT NUMBER WHICH DETERMINES WHERE TO WRITE REGULAR OUTPUT C LP = UNIT NUMBER WHICH DETERMINES WHERE TO WRITE ERROR OUPUT BLOCK DATA CGCD COMMON /CGDD/MP,LP INTEGER LP,MP DATA MP,LP/44,44/ END C LAST LINE OF CGFAM c*********************************************************************** c BLAS parts: c*********************************************************************** c ddot c daxpy c dasum c dcopy c*********************************************************************** c matmult c*********************************************************************** double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end c*********************************************************************** subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end c*********************************************************************** double precision function dasum(n,dx,incx) c c takes the sum of the absolute values. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dtemp integer i,incx,m,mp1,n,nincx c dasum = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dtemp = dtemp + dabs(dx(i)) 10 continue dasum = dtemp return c c code for increment equal to 1 c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dabs(dx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) 50 continue 60 dasum = dtemp return end c*********************************************************************** subroutine dcopy(n,dx,incx,dy,incy) c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 dy(i) = dx(i) dy(i + 1) = dx(i + 1) dy(i + 2) = dx(i + 2) dy(i + 3) = dx(i + 3) dy(i + 4) = dx(i + 4) dy(i + 5) = dx(i + 5) dy(i + 6) = dx(i + 6) 50 continue return end c *********************************************************** subroutine matmult(m1,d1,m2,d2,mout,dim) c matrix multiplication: c m1[d1,dim]*m2[dim,d2]=mout[d1,d2] real*8 m1,m2,mout,swap integer d1,d2,dim,i,j,k dimension m1(d1,dim),m2(dim,d2),mout(d1,d2),swap(d1,d2) call mat_init(swap,d1,d2,0.0d0) do 10,i=1,d1 do 20,j=1,d2 do 30,k=1,dim swap(i,j)=swap(i,j)+m1(i,k)*m2(k,j) 30 continue 20 continue 10 continue call mat_copy(swap,mout,d1,d2) return end cc subroutine mat_init(mat,nrow,ncol,value) c fills the NROW * NCOL matrix MAT with VALUE real*8 mat,value integer nrow,ncol,i,j dimension mat(nrow,ncol) do 10,i=1,nrow do 20,j=1,ncol mat(i,j)=value 20 continue 10 continue return end cc subroutine mat_copy(matin,matout,nrow,ncol) c copy MATIN(NROW,NCOL) to MATOUT(NROW,NCOL) real*8 matin,matout integer nrow,ncol,i,j dimension matin(nrow,ncol),matout(nrow,ncol) do 10,i=1,nrow do 20,j=1,ncol matout(i,j)=matin(i,j) 20 continue 10 continue return end