PROGRAM A22cgGS3 C CCC for ALA potential N=66 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C W.Quapp 30.11.2004/11.08.2006 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 (modified by wq) CG+ of Nocedal et al. C the actual point, to start the correction, is Y(LA,..) c c version for use of z_mat coordinates additionally ! Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Integer istatus,ITALL,N3 Integer LA,NDIM,NDIM1,L PARAMETER(NDIM=66,L=32,NDIM1=67) c !! NDIM auch weiter unten einstellen!! c steplength-suche mit '2308' !! 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) REAL*8 Y(L,NDIM),VV(NDIM),D1,Yold(NDIM) REAL*8 Rs(NDIM),energy,Cx,gnorm,Pena,xyzmass(NDIM/3) REAL*8 ProMat(NDIM,NDIM),DiStart,DiCurr logical finish, DiLog integer iprint(2),iflag,icall,n,method,iter,nfun,IREST Integer JJ,I,J,Np,Np1,IterMax,J1,Natoms(NDIM/3) Character*2 CharAtoms(NDIM/3) 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'/ Data Natoms/1,6,1,6,1,8,7,6,1,6,1,6,8,7,6,1,1,1,1,1,1,1/ c OPEN(9,FILE='ala22chain.dat') OPEN(10,FILE='ala22Rs.dat') OPEN(11,FILE='ala22param.dat') OPEN(12,FILE='ala22proj.dat') OPEN(13,FILE='ala22point.dat') OPEN(14,FILE='ala22ener.dat') OPEN(15,FILE='ala22grad.dat') OPEN(16,FILE='ala22engs.weg') OPEN(17,FILE='ala22lamb.dat') OPEN(19,FILE='ala22grno.dat') OPEN(21,FILE='ala22iflag.dat') OPEN(22,FILE='ala22cgplFi.dat') OPEN(23,FILE='ala22cggsW.dat') OPEN(25,FILE='ala22dddd.dat') OPEN(26,FILE='ala22gold.dat') OPEN(28,FILE='ala22yold.dat') OPEN(33,FILE='ala22poig.dat') OPEN(35,FILE='ala22Distart.dat') OPEN(38,FILE='ala22chain.mol',status='unknown',access='append') OPEN(44,FILE='ala22proto.txt',status='unknown',access='append') c Read problem input information Pena= 1.0000001d+0 N = NDIM N1 = N+1 N3 = N/3 C N1 for Lagrange formulation of one additional constraint C N3: number of atoms, Pena: now useless constant IPRINT(1)= 0 IPRINT(2)= 0 XTOL= 11.11111D-11 ITERMAX= 2* N c Note: EPS1= 1.0D-5 of CG+ is very to small here !!! C CG+ method = 3 IREST = 1 rewind 21 read(21,*) IFLAG c write( 6,*) ' iflag is ', IFLAG c write(44,*) ' 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,N1) rewind 26 read(26,50) (Gold(J),J=1,N1) rewind 28 read(28,50) (Yold(J),J=1,N) rewind 17 read(17,50) XX(N1) endif istatus=0 rewind 11 read(11,*) JJ, LA, eps ,ITALL, NFUN ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc itGS=itGS+1 -->> now other count: by NFUN of CG+ cc JJ=L !! rewind 11 write(11,*) L, LA, eps ,ITALL, NFUN write(44,*)' L=',L,' LA= ',LA,' ITall = ',ITALL,'+',NFUN c Print parameters if (iprint(1).ge.0 .AND. LA.lt.3 .AND. IFLAG.eq.0)then write(*,820) write(*,840) N1, method, irest write(44,820) write(44,840) N1, method, irest endif c c 20 CONTINUE c NOTE: CG+ cycle replaced by script organisation c Call the function and gradient values here c 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,*) (X(I),I=1,N) rewind 14 READ(14,*) energy WRITE(44,*)'energy=',energy rewind 9 Do 4, J=1,L Read(9,*) J1 Read(9,*) Do 4, JJ=1,N3 4 Read(9,555) CharAtoms(JJ),(Y(J,3*(JJ-1)+I),I=1,3) 555 Format(1X,A3, 3F16.8) If(IFLAG.EQ.0) then Do 444, I=1,N 444 Y(LA,I)=X(I) endif c rewind 15 Do 1, J=1,N3 1 READ(15,*,END=111) J1, JJ, (VV( 3*(J-1)+I ),I=1,3) cc g3 call PERROR("I/O processing") cc g3 write(6,*) ' for gradient ' IF(J1 .ne. N3 ) WRITE(44,*) ' ERROR in gradient reading ! ' Do 5 I=1,N 5 VV(I)=-VV(I) c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c for using the special file with Chartesian Grad output c from Gaussian03 !! Note: - sign from gaussian gradient for VV !! c and: pull back the mass weighting of Gaussian03 gradient c we use all (coordinates, PES, Gradient) without mass weighting ! c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c do 31,I=1,N3 JJ=3*(I-1) xyzmass(I)=Dsqrt( Natoms(I)*2.0d0 ) IF(Natoms(I).ne.1) then VV(JJ+1) = VV(JJ+1)/ xyzmass(I) VV(JJ+2) = VV(JJ+2)/ xyzmass(I) VV(JJ+3) = VV(JJ+3)/ xyzmass(I) ENDIF 31 continue 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 project out search direction If(ICALL.eq.1) then Do 133 J=1,N 133 Proj(J)=RS(J) call PROJECTION(N,Proj,ProMat) D1= DSQRT(DDOT(N,Proj,1,Proj,1)) WRITE(6,*) ' projected RS ', D1 WRITE(44,*)' projected RS ', D1 Endif Do 13 J=1,N Proj(J)=VV(J) 13 XX(J) =X(J) cc enforce before every CG+ step the "true" lambda XX(N1)= gnorm / Pena c estimation for lambda of Lagrange ansatz (with + sign): 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 cccccccccccccccccccccccccccccccccccccccc C Newton trajectory with RGF Projector c and Projector of rotations: C Projektion of Gradient cccccccccccccccccccccccccccccccccccccccc call PROJECTION(N,Proj,ProMat) D1= DSQRT(DDOT(N,Proj,1,Proj,1)) WRITE(6,*) ' projected grad norm ', D1 If(D1 .lt. eps) THEN WRITE(44,*)'at IteGS=',NFUN,' ok: projNorm',D1,' less then', 1 ' tolerance ', eps WRITE(44,*)' SUCCESSFUL CONVERGENCE (no ERROR) ' WRITE(6, *)' IteGS=',NFUN,' == ok - proj norm ',D1 ITALL = ITALL+NFUN NFUN=0 IFLAG=0 istatus=1 c GOTO 111 ELSE WRITE(44,*)' projected grad norm ', D1 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 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, put to 1 usually cc Cx=0.0d0 Do 23 J=1,N 23 Cx= Cx + (X(J)-Yold(J))*Rs(J) * Pena write(44,*) ' C(x) = ' , Cx c Energy and Gradient for constrained problem: Fun= energy - XX(N1)*Cx Do 18 J=1,N CCCC ultimate TEST: use proj gradient !!! C18 G(J)= VV(J) - XX(N1)*Rs(J)*Pena 18 G(J)= Proj(J) CCCC wq 18.07.2006 G(N1)= -Cx c TESTs c IF(ICALL.eq.3) c 1 write(44,*) 'G(J) = VV(J) - XX(N1)*Rs(J)',(G(J),J=1,N1) c D1=0.0d0 c Do 1388 J=1,N 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,N 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(N1,XX,Fun,G,D,GOLD,IPRINT,W, * IFLAG,IREST,METHOD,FINISH,ProMat) 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 c write(44,*) ' NFUN= ', NFUN,' at IFLAG',IFLAG, c 1 ': lambda given by CG+ ',XX(N1) rewind 11 write(11,*) L, LA, eps ,ITALL, NFUN ccccccccccccccccccccccccccccccccccccccccccccc C step normalization if(IFLAG.eq.2) then Read(35,*) Distart,DiLog DiCurr=0.0d0 Do 717, I=1,N 717 DiCurr=DiCurr+(Y(LA-1,I)-XX(I))**2 DiCurr=Dsqrt(DiCurr) if(DiLog) then Do 718, I=1,N 718 XX(I)=Y(LA-1,I)*(1.0d0-Distart/DiCurr)+Distart/DiCurr*XX(I) Write(44,*)'DistanceCorrection: ', Distart/DiCurr else Write(44,*)'NormedDistanceLastNode: ',Distart/DiCurr Write(44,*)'AbsDistanceLastNode : ',DiCurr endif endif cccccccccccccccccccccccccccccccccccccccccccccc Do 201,J=1,N 201 X(J) = XX(J) if(IFLAG.ne.2) write(44,*) c actual point approximation: rewind 13 WRITE(13,50) (XX(I),I=1,N) rewind 17 WRITE(17,50) XX(N1) c cc for gaussian input file c rewind 32 c READ(32,*) (Natoms(I),I=1,N3) rewind 33 Do 7, J=1,N3 7 WRITE(33,*) Natoms(J), (X(3*(J-1)+I ),I=1,3) WRITE(33,*) ' ' c store CG+ data for script organisation: rewind 23 write(23,50) (W(J),J=1,N1) rewind 25 write(25,50) (D(J),J=1,N1) rewind 26 write(26,50) (Gold(J),J=1,N1) rewind 21 write(21,*) IFLAG write( 6,*) ' iflag nach CG+ ', IFLAG c write(44,*) ' iflag nach CG+ ', IFLAG rewind 22 write(22,*) FINISH,ICALL,NFUN ccc 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 WRITE(44,*)' Line search LOOP, ICALL= ',ICALL WRITE(44,*)' lambda =',XX(N1) rewind 22 write(22,*) FINISH,ICALL,NFUN rewind 11 write(11,*) L, LA, eps ,ITALL, NFUN istatus=0 call exit(istatus) STOP cc GO TO 20 --> replaced by script organisation ENDIF IF(IFLAG.EQ.2) THEN c Termination Test. You 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 FINISH = .TRUE. write(44,*) ' finish = true ' ITALL=ITALL+NFUN NFUN=0 GO TO 30 ENDIF cc wq IF(DABS(Proj(I)).GT.TLEV) THEN write(44,*)'Line search finished' write(44,*)'Iter because ProjGradient(',I,') = ', Proj(I) GO TO 30 ELSE GO TO 41 ENDIF ENDIF IF(IFLAG.eq.0) then write(44,*)'ICALL= ',ICALL,' under succsessful IFLAG=',IFLAG istatus=1 ITALL = ITALL+NFUN NFUN=0 rewind 21 write(21,*) iflag ENDIF IF(IFLAG.lt.0) then write(6,*)' Cancel Loop: not convergent in CG+ ' write(44,*)' ERROR : Corr not convergent in CG+ ' write(44,*)' IFLAG: ', IFLAG ITALL = ITALL+NFUN NFUN=0 rewind 21 write(21,*) iflag istatus=0 endif IF(ICALL.ge.ITERMAX) then write(6,*)' ERROR : ITERMAX-Abbruch in CG+ ' write(44,*)' ERROR : ITERMAX-Abbruch in CG+ ' istatus=1 ITALL=ITALL+NFUN NFUN=0 write(44,*)'Stop Corr: ICALL=ITERMAX, IFLAG=',IFLAG endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 111 continue ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (iprint(1).ge.0.and.iflag.ge.0) then write (*,890) Fun write (*,*) ' at LA ', LA write (44,890) Fun end if write(44,890) energy write(44,50) (X(I),I=1,N) write(44,*) ' C(x)= ', Cx/Pena write(44,*) write(38,*)' ',N3 write(38,*) Do 16, J=1,N3 16 WRITE(38,555) CharAtoms(J),(X(3*(J-1)+I ),I=1,3) if(LA.lt.L) write(44,*) ' new step: ' If(LA.gt.L) then write(6,*)' ERROR : LA-Abbruch in CG+ ' write(44,*)' ERROR : LA-Abbruch in CG+ ' ITall = ITall + NFUN NFUN=0 istatus=1 endif rewind 21 write(21,*) iflag c Istatus controls the script loops, new chain point becomes: write(6,*) ' finish is ', finish write(44,*) ' finish is ', finish write(44,*) ' ' write(44,*) ' ' Do 117 I=1,N 117 Y(LA,I)= XX(I) rewind 9 Do 118, J=1,L write(9,*) ' ',N3 write(9,*) ' ' Do 118, JJ=1,N3 118 write(9,555) CharAtoms(JJ),(Y(J,3*(JJ-1)+I),I=1,3) WRITE(6,*)' energy= ',energy WRITE(44,*)' energy= ',energy WRITE(16,*) '{ ' IF(LA.lt.L) WRITE(16,*) ' ', energy,' , ' IF(LA.eq.L) WRITE(16,*) ' ', energy,' } ' istatus=1 write(44,*)' istatus = ', istatus c rewind 11 write(11,*) L, LA, eps ,ITALL, NFUN write(44,*)' at the end of frame of CGFAM ' write(44,*)' L=', L,' LA= ', LA ,' ITall = ',ITALL, ' + ',NFUN write(44,*)' istatus = ', istatus call exit(istatus) STOP 50 FORMAT(1X,3F20.10) 820 format ( ' Conjugate Gradient Minimization Routine') 840 format ( ' n =', i6, * ' CG+ 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) 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 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 Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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='ala22proto.txt',status='unknown',access='append') OPEN(31,FILE='ala22cvsRest.dat') C wq for script organisation 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 .LT. 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 c wq IF (NFEV .GE. MAXFEV) INFO = 3 IF (NFEV .GE. MAXFEV) then INFO = 1 write(44,*) ' ERROR by MAXFEV ' ENDIF 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. END 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) 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 ,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/ c wq store for script organisation: OPEN(24,FILE='ala22cgRest.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 = N-1 cc Projector use to enforce the C(x)=0 condition, automatically c call PROJECTION(oldN,D,ProMat) c GNORM= DSQRT(DDOT(N,D,1,D,1)) cc write(44,*) 'N1-norm of D ', GNORM cc here is used the shortened length of projected D !! 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-20 c STPMIN= 1.0D-20 c STPMAX= 1.0D+20 c bisher = 20 cc wq new: MAXFEV= 8 stpmin= 0.23081947 stpmax= 1.750000000023081947d+0 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 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 70 GOLD(I)= G(I) DG= DDOT(N,D,1,G,1) 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 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. c 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 shut down some test which are to strong: INFO5 IF(INFO.EQ.4 .or. INFO.EQ.5 .or. INFO.EQ.6) goto 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 = N-1 call PROJECTION(oldN,D,ProMat) c DG0= DGOLD*STP C C RETURN TO DRIVER FOR TERMINATION TEST C ------------------------------------- C GNORM=DSQRT(DDOT(N,G,1,G,1)) c 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