home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Frostbyte's 1980s DOS Shareware Collection
/
floppyshareware.zip
/
floppyshareware
/
DOOG
/
PCSSP1.ZIP
/
EXTRMFCT.ZIP
/
DFMFP.FOR
< prev
next >
Wrap
Text File
|
1985-11-29
|
10KB
|
325 lines
C
C ..................................................................
C
C SUBROUTINE DFMFP
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 DFMFP(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 ARG,VAL AND GRAD MUST BE OF DOUBLE PRECISION.
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 DOUBLE PRECISION VECTOR.
C F - SINGLE VARIABLE CONTAINING THE MINIMUM FUNCTION
C VALUE ON RETURN, I.E. F=F(X).
C DOUBLE PRECISION VARIABLE.
C G - VECTOR OF DIMENSION N CONTAINING THE GRADIENT
C VECTOR CORRESPONDING TO THE MINIMUM ON RETURN,
C I.E. G=G(X).
C DOUBLE PRECISION VECTOR.
C EST - IS AN ESTIMATE OF THE MINIMUM FUNCTION VALUE.
C SINGLE PRECISION VARIABLE.
C EPS - TESTVALUE REPRESENTING THE EXPECTED ABSOLUTE ERROR.
C A REASONABLE CHOICE IS 10**(-16), I.E.
C SOMEWHAT GREATER THAN 10**(-D), WHERE D IS THE
C NUMBER OF SIGNIFICANT DIGITS IN FLOATING POINT
C REPRESENTATION.
C SINGLE PRECISION VARIABLE.
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 DOUBLE PRECISION ARRAY.
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 DFMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)
C
C DIMENSIONED DUMMY VARIABLES
DIMENSION H(1),X(1),G(1)
DOUBLE PRECISION X,F,FX,FY,OLDF,HNRM,GNRM,H,G,DX,DY,ALFA,DALFA,
1AMBDA,T,Z,W
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.D0
NJ=N-J
IF(NJ)5,5,2
2 DO 3 L=1,NJ
KL=K+L
3 H(KL)=0.D0
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.D0
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.D0
HNRM=0.D0
GNRM=0.D0
C
C CALCULATE DIRECTIONAL DERIVATIVE AND TESTVALUES FOR DIRECTION
C VECTOR H AND GRADIENT VECTOR G.
DO 10 J=1,N
HNRM=HNRM+DABS(H(J))
GNRM=GNRM+DABS(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.D0*(EST-F)/DY
AMBDA=1.D0
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.D0
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.D0
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.D10)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.D0
23 IF(AMBDA)24,36,24
24 Z=3.D0*(FX-FY)/AMBDA+DX+DY
ALFA=DMAX1(DABS(Z),DABS(DX),DABS(DY))
DALFA=Z/ALFA
DALFA=DALFA*DALFA-DX/ALFA*DY/ALFA
IF(DALFA)51,25,25
25 W=ALFA*DSQRT(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.D0
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.D0
Z=0.D0
DO 40 J=1,N
K=N+J
W=H(K)
K=K+N
T=T+DABS(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.D0
DO 47 J=1,N
K=J+N3
W=0.D0
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