C C .................................................................. C C SUBROUTINE FMFP C C PURPOSE C TO FIND A LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES C BY THE METHOD OF FLETCHER AND POWELL C C USAGE C CALL FMFP(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 N*(N+7)/2. 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 M.J.D. POWELL, A RAPID DESCENT METHOD FOR C MINIMIZATION, C COMPUTER JOURNAL VOL.6, ISS. 2, 1963, PP.163-168. C C .................................................................. C SUBROUTINE FMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H) C C DIMENSIONED DUMMY VARIABLES DIMENSION H(1),X(1),G(1) C C COMPUTE FUNCTION VALUE AND GRADIENT VECTOR FOR INITIAL ARGUMENT CALL FUNCT(N,X,F,G) C C RESET ITERATION COUNTER AND GENERATE IDENTITY MATRIX IER=0 KOUNT=0 N2=N+N N3=N2+N N31=N3+1 1 K=N31 DO 4 J=1,N H(K)=1. NJ=N-J IF(NJ)5,5,2 2 DO 3 L=1,NJ KL=K+L 3 H(KL)=0. 4 K=KL+1 C C START ITERATION LOOP 5 KOUNT=KOUNT +1 C C SAVE FUNCTION VALUE, ARGUMENT VECTOR AND GRADIENT VECTOR OLDF=F DO 9 J=1,N K=N+J H(K)=G(J) K=K+N H(K)=X(J) C C DETERMINE DIRECTION VECTOR H K=J+N3 T=0. DO 8 L=1,N T=T-G(L)*H(K) IF(L-J)6,7,7 6 K=K+N-L GO TO 8 7 K=K+1 8 CONTINUE 9 H(J)=T C C CHECK WHETHER FUNCTION WILL DECREASE STEPPING ALONG H. DY=0. HNRM=0. GNRM=0. C C CALCULATE DIRECTIONAL DERIVATIVE AND TESTVALUES FOR DIRECTION C VECTOR H AND GRADIENT VECTOR G. DO 10 J=1,N HNRM=HNRM+ABS(H(J)) GNRM=GNRM+ABS(G(J)) 10 DY=DY+H(J)*G(J) C C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTIONAL C DERIVATIVE APPEARS TO BE POSITIVE OR ZERO. IF(DY)11,51,51 C C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTION C VECTOR H IS SMALL COMPARED TO GRADIENT VECTOR G. 11 IF(HNRM/GNRM-EPS)51,51,12 C C SEARCH MINIMUM ALONG DIRECTION H C C SEARCH ALONG H FOR POSITIVE DIRECTIONAL DERIVATIVE 12 FY=F ALFA=2.*(EST-F)/DY AMBDA=1. C C USE ESTIMATE FOR STEPSIZE ONLY IF IT IS POSITIVE AND LESS THAN C 1. OTHERWISE TAKE 1. AS STEPSIZE IF(ALFA)15,15,13 13 IF(ALFA-AMBDA)14,15,15 14 AMBDA=ALFA 15 ALFA=0. C C SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT 16 FX=FY DX=DY C C STEP ARGUMENT ALONG H DO 17 I=1,N 17 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 IS POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND DY=0. DO 18 I=1,N 18 DY=DY+G(I)*H(I) IF(DY)19,36,22 C C TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT C A MINIMUM HAS BEEN PASSED 19 IF(FY-FX)20,22,22 C C REPEAT SEARCH AND DOUBLE STEPSIZE FOR FURTHER SEARCHES 20 AMBDA=AMBDA+ALFA ALFA=AMBDA C END OF SEARCH LOOP C C TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE IF(HNRM*AMBDA-1.E10)16,16,21 C C LINEAR SEARCH TECHNIQUE INDICATES THAT NO MINIMUM EXISTS 21 IER=2 RETURN 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 22 T=0. 23 IF(AMBDA)24,36,24 24 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)51,25,25 25 W=ALFA*SQRT(DALFA) ALFA=DY-DX+W+W IF(ALFA) 250,251,250 250 ALFA=(DY-Z+W)/ALFA GO TO 252 251 ALFA=(Z+DY-W)/(Z+DX+Z+DY) 252 ALFA=ALFA*AMBDA DO 26 I=1,N 26 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)27,27,28 27 IF(F-FY)36,36,28 28 DALFA=0. DO 29 I=1,N 29 DALFA=DALFA+G(I)*H(I) IF(DALFA)30,33,33 30 IF(F-FX)32,31,33 31 IF(DX-DALFA)32,36,32 32 FX=F DX=DALFA T=ALFA AMBDA=ALFA GO TO 23 33 IF(FY-F)35,34,35 34 IF(DY-DALFA)35,36,35 35 FY=F DY=DALFA AMBDA=AMBDA-ALFA GO TO 22 C C TERMINATE, IF FUNCTION HAS NOT DECREASED DURING LAST ITERATION 36 IF(OLDF-F+EPS)51,38,38 C C COMPUTE DIFFERENCE VECTORS OF ARGUMENT AND GRADIENT FROM C TWO CONSECUTIVE ITERATIONS 38 DO 37 J=1,N K=N+J H(K)=G(J)-H(K) K=N+K 37 H(K)=X(J)-H(K) C C TEST LENGTH OF ARGUMENT DIFFERENCE VECTOR AND DIRECTION VECTOR C IF AT LEAST N ITERATIONS HAVE BEEN EXECUTED. TERMINATE, IF C BOTH ARE LESS THAN EPS IER=0 IF(KOUNT-N)42,39,39 39 T=0. Z=0. DO 40 J=1,N K=N+J W=H(K) K=K+N T=T+ABS(H(K)) 40 Z=Z+W*H(K) IF(HNRM-EPS)41,41,42 41 IF(T-EPS)56,56,42 C C TERMINATE, IF NUMBER OF ITERATIONS WOULD EXCEED LIMIT 42 IF(KOUNT-LIMIT)43,50,50 C C PREPARE UPDATING OF MATRIX H 43 ALFA=0. DO 47 J=1,N K=J+N3 W=0. DO 46 L=1,N KL=N+L W=W+H(KL)*H(K) IF(L-J)44,45,45 44 K=K+N-L GO TO 46 45 K=K+1 46 CONTINUE K=N+J ALFA=ALFA+W*H(K) 47 H(J)=W C C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF RESULTS C ARE NOT SATISFACTORY IF(Z*ALFA)48,1,48 C C UPDATE MATRIX H 48 K=N31 DO 49 L=1,N KL=N2+L DO 49 J=L,N NJ=N2+J H(K)=H(K)+H(KL)*H(NJ)/Z-H(L)*H(J)/ALFA 49 K=K+1 GO TO 5 C END OF ITERATION LOOP C C NO CONVERGENCE AFTER LIMIT ITERATIONS 50 IER=1 RETURN C C RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS 51 DO 52 J=1,N K=N2+J 52 X(J)=H(K) CALL FUNCT(N,X,F,G) C C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DERIVATIVE C FAILS TO BE SUFFICIENTLY SMALL IF(GNRM-EPS)55,55,53 C C TEST FOR REPEATED FAILURE OF ITERATION 53 IF(IER)56,54,54 54 IER=-1 GOTO 1 55 IER=0 56 RETURN END