C C .................................................................. C C SUBROUTINE FMCG C C PURPOSE C TO FIND A LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES C BY THE METHOD OF CONJUGATE GRADIENTS C C USAGE C CALL FMCG(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H) C C DESCRIPTION OF PARAMETERS C FUNCT - USER-WRITTEN SUBROUTINE CONCERNING THE FUNCTION TO C BE MINIMIZED. IT MUST BE OF THE FORM C SUBROUTINE FUNCT(N,ARG,VAL,GRAD) C AND MUST SERVE THE FOLLOWING PURPOSE C FOR EACH N-DIMENSIONAL ARGUMENT VECTOR ARG, C FUNCTION VALUE AND GRADIENT VECTOR MUST BE COMPUTED C AND, ON RETURN, STORED IN VAL AND GRAD RESPECTIVELY C N - NUMBER OF VARIABLES C X - VECTOR OF DIMENSION N CONTAINING THE INITIAL C ARGUMENT WHERE THE ITERATION STARTS. ON RETURN, C X HOLDS THE ARGUMENT CORRESPONDING TO THE C COMPUTED MINIMUM FUNCTION VALUE C F - SINGLE VARIABLE CONTAINING THE MINIMUM FUNCTION C VALUE ON RETURN, I.E. F=F(X). C G - VECTOR OF DIMENSION N CONTAINING THE GRADIENT C VECTOR CORRESPONDING TO THE MINIMUM ON RETURN, C I.E. G=G(X). C EST - IS AN ESTIMATE OF THE MINIMUM FUNCTION VALUE. C EPS - TESTVALUE REPRESENTING THE EXPECTED ABSOLUTE ERROR. C A REASONABLE CHOICE IS 10**(-6), I.E. C SOMEWHAT GREATER THAN 10**(-D), WHERE D IS THE C NUMBER OF SIGNIFICANT DIGITS IN FLOATING POINT C REPRESENTATION. C LIMIT - MAXIMUM NUMBER OF ITERATIONS. C IER - ERROR PARAMETER C IER = 0 MEANS CONVERGENCE WAS OBTAINED C IER = 1 MEANS NO CONVERGENCE IN LIMIT ITERATIONS C IER =-1 MEANS ERRORS IN GRADIENT CALCULATION C IER = 2 MEANS LINEAR SEARCH TECHNIQUE INDICATES C IT IS LIKELY THAT THERE EXISTS NO MINIMUM. C H - WORKING STORAGE OF DIMENSION 2*N. C C REMARKS C I) THE SUBROUTINE NAME REPLACING THE DUMMY ARGUMENT FUNCT C MUST BE DECLARED AS EXTERNAL IN THE CALLING PROGRAM. C II) IER IS SET TO 2 IF , STEPPING IN ONE OF THE COMPUTED C DIRECTIONS, THE FUNCTION WILL NEVER INCREASE WITHIN C A TOLERABLE RANGE OF ARGUMENT. C IER = 2 MAY OCCUR ALSO IF THE INTERVAL WHERE F C INCREASES IS SMALL AND THE INITIAL ARGUMENT WAS C RELATIVELY FAR AWAY FROM THE MINIMUM SUCH THAT THE C MINIMUM WAS OVERLEAPED. THIS IS DUE TO THE SEARCH C TECHNIQUE WHICH DOUBLES THE STEPSIZE UNTIL A POINT C IS FOUND WHERE THE FUNCTION INCREASES. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C FUNCT C C METHOD C THE METHOD IS DESCRIBED IN THE FOLLOWING ARTICLE C R.FLETCHER AND C.M.REEVES, FUNCTION MINIMIZATION BY C CONJUGATE GRADIENTS, C COMPUTER JOURNAL VOL.7, ISS.2, 1964, PP.149-154. C C .................................................................. C SUBROUTINE FMCG(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H) C C DIMENSIONED DUMMY VARIABLES DIMENSION X(1),G(1),H(1) C C C COMPUTE FUNCTION VALUE AND GRADIENT VECTOR FOR INITIAL ARGUMENT CALL FUNCT(N,X,F,G) C C RESET ITERATION COUNTER KOUNT=0 IER=0 N1=N+1 C C START ITERATION CYCLE FOR EVERY N+1 ITERATIONS 1 DO 43 II=1,N1 C C STEP ITERATION COUNTER AND SAVE FUNCTION VALUE KOUNT=KOUNT+1 OLDF=F C C COMPUTE SQUARE OF GRADIENT AND TERMINATE IF ZERO GNRM=0. DO 2 J=1,N 2 GNRM=GNRM+G(J)*G(J) IF(GNRM)46,46,3 C C EACH TIME THE ITERATION LOOP IS EXECUTED , THE FIRST STEP WILL C BE IN DIRECTION OF STEEPEST DESCENT 3 IF(II-1)4,4,6 4 DO 5 J=1,N 5 H(J)=-G(J) GO TO 8 C C FURTHER DIRECTION VECTORS H WILL BE CHOOSEN CORRESPONDING C TO THE CONJUGATE GRADIENT METHOD 6 AMBDA=GNRM/OLDG DO 7 J=1,N 7 H(J)=AMBDA*H(J)-G(J) C C COMPUTE TESTVALUE FOR DIRECTIONAL VECTOR AND DIRECTIONAL C DERIVATIVE 8 DY=0. HNRM=0. DO 9 J=1,N K=J+N C C SAVE ARGUMENT VECTOR H(K)=X(J) HNRM=HNRM+ABS(H(J)) 9 DY=DY+H(J)*G(J) C C CHECK WHETHER FUNCTION WILL DECREASE STEPPING ALONG H AND C SKIP LINEAR SEARCH ROUTINE IF NOT IF(DY)10,42,42 C C COMPUTE SCALE FACTOR USED IN LINEAR SEARCH SUBROUTINE 10 SNRM=1./HNRM C C SEARCH MINIMUM ALONG DIRECTION H C C SEARCH ALONG H FOR POSITIVE DIRECTIONAL DERIVATIVE FY=F ALFA=2.*(EST-F)/DY AMBDA=SNRM C C USE ESTIMATE FOR STEPSIZE ONLY IF IT IS POSITIVE AND LESS THAN C SNRM. OTHERWISE TAKE SNRM AS STEPSIZE. IF(ALFA)13,13,11 11 IF(ALFA-AMBDA)12,13,13 12 AMBDA=ALFA 13 ALFA=0. C C SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT 14 FX=FY DX=DY C C STEP ARGUMENT ALONG H DO 15 I=1,N 15 X(I)=X(I)+AMBDA*H(I) C C COMPUTE FUNCTION VALUE AND GRADIENT FOR NEW ARGUMENT CALL FUNCT(N,X,F,G) FY=F C C COMPUTE DIRECTIONAL DERIVATIVE DY FOR NEW ARGUMENT. TERMINATE C SEARCH, IF DY POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND DY=0. DO 16 I=1,N 16 DY=DY+G(I)*H(I) IF(DY)17,38,20 C C TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT C A MINIMUM HAS BEEN PASSED 17 IF(FY-FX)18,20,20 C C REPEAT SEARCH AND DOUBLE STEPSIZE FOR FURTHER SEARCHES 18 AMBDA=AMBDA+ALFA ALFA=AMBDA C C TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE IF(HNRM*AMBDA-1.E10)14,14,19 C C LINEAR SEARCH TECHNIQUE INDICATES THAT NO MINIMUM EXISTS 19 IER=2 C C RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS F=OLDF DO 100 J=1,N G(J)=H(J) K=N+J 100 X(J)=H(K) RETURN C END OF SEARCH LOOP C C INTERPOLATE CUBICALLY IN THE INTERVAL DEFINED BY THE SEARCH C ABOVE AND COMPUTE THE ARGUMENT X FOR WHICH THE INTERPOLATION C POLYNOMIAL IS MINIMIZED C 20 T=0. 21 IF(AMBDA)22,38,22 22 Z=3.*(FX-FY)/AMBDA+DX+DY ALFA=AMAX1(ABS(Z),ABS(DX),ABS(DY)) DALFA=Z/ALFA DALFA=DALFA*DALFA-DX/ALFA*DY/ALFA IF(DALFA)23,27,27 C C RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS 23 DO 24 J=1,N K=N+J 24 X(J)=H(K) CALL FUNCT(N,X,F,G) C C TEST FOR REPEATED FAILURE OF ITERATION 25 IF(IER)47,26,47 26 IER=-1 GOTO 1 27 W=ALFA*SQRT(DALFA) ALFA=DY-DX+W+W IF(ALFA)270,271,270 270 ALFA=(DY-Z+W)/ALFA GO TO 272 271 ALFA=(Z+DY-W)/(Z+DX+Z+DY) 272 ALFA=ALFA*AMBDA DO 28 I=1,N 28 X(I)=X(I)+(T-ALFA)*H(I) C C TERMINATE, IF THE VALUE OF THE ACTUAL FUNCTION AT X IS LESS C THAN THE FUNCTION VALUES AT THE INTERVAL ENDS. OTHERWISE REDUCE C THE INTERVAL BY CHOOSING ONE END-POINT EQUAL TO X AND REPEAT C THE INTERPOLATION. WHICH END-POINT IS CHOOSEN DEPENDS ON THE C VALUE OF THE FUNCTION AND ITS GRADIENT AT X C CALL FUNCT(N,X,F,G) IF(F-FX)29,29,30 29 IF(F-FY)38,38,30 C C COMPUTE DIRECTIONAL DERIVATIVE 30 DALFA=0. DO 31 I=1,N 31 DALFA=DALFA+G(I)*H(I) IF(DALFA)32,35,35 32 IF(F-FX)34,33,35 33 IF(DX-DALFA)34,38,34 34 FX=F DX=DALFA T=ALFA AMBDA=ALFA GO TO 21 35 IF(FY-F)37,36,37 36 IF(DY-DALFA)37,38,37 37 FY=F DY=DALFA AMBDA=AMBDA-ALFA GO TO 20 C C TERMINATE, IF FUNCTION HAS NOT DECREASED DURING LAST ITERATION C OTHERWISE SAVE GRADIENT NORM 38 IF(OLDF-F+EPS)19,25,39 39 OLDG=GNRM C C COMPUTE DIFFERENCE OF NEW AND OLD ARGUMENT VECTOR T=0. DO 40 J=1,N K=J+N H(K)=X(J)-H(K) 40 T=T+ABS(H(K)) C C TEST LENGTH OF DIFFERENCE VECTOR IF AT LEAST N+1 ITERATIONS C HAVE BEEN EXECUTED. TERMINATE, IF LENGTH IS LESS THAN EPS IF(KOUNT-N1)42,41,41 41 IF(T-EPS)45,45,42 C C TERMINATE, IF NUMBER OF ITERATIONS WOULD EXCEED LIMIT 42 IF(KOUNT-LIMIT)43,44,44 43 IER=0 C END OF ITERATION CYCLE C C START NEXT ITERATION CYCLE GO TO 1 C C NO CONVERGENCE AFTER LIMIT ITERATIONS 44 IER=1 IF(GNRM-EPS)46,46,47 C C TEST FOR SUFFICIENTLY SMALL GRADIENT 45 IF(GNRM-EPS)46,46,25 46 IER=0 47 RETURN END