home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / EXTRMFCT.ZIP / FMFP.FOR < prev   
Text File  |  1985-11-29  |  10KB  |  316 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE FMFP
  5. C
  6. C        PURPOSE
  7. C           TO FIND A LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES
  8. C           BY THE METHOD OF FLETCHER AND POWELL
  9. C
  10. C        USAGE
  11. C           CALL FMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           FUNCT  - USER-WRITTEN SUBROUTINE CONCERNING THE FUNCTION TO
  15. C                    BE MINIMIZED. IT MUST BE OF THE FORM
  16. C                    SUBROUTINE FUNCT(N,ARG,VAL,GRAD)
  17. C                    AND MUST SERVE THE FOLLOWING PURPOSE
  18. C                    FOR EACH N-DIMENSIONAL ARGUMENT VECTOR  ARG,
  19. C                    FUNCTION VALUE AND GRADIENT VECTOR MUST BE COMPUTED
  20. C                    AND, ON RETURN, STORED IN VAL AND GRAD RESPECTIVELY
  21. C           N      - NUMBER OF VARIABLES
  22. C           X      - VECTOR OF DIMENSION N CONTAINING THE INITIAL
  23. C                    ARGUMENT WHERE THE ITERATION STARTS. ON RETURN,
  24. C                    X HOLDS THE ARGUMENT CORRESPONDING TO THE
  25. C                    COMPUTED MINIMUM FUNCTION VALUE
  26. C           F      - SINGLE VARIABLE CONTAINING THE MINIMUM FUNCTION
  27. C                    VALUE ON RETURN, I.E. F=F(X).
  28. C           G      - VECTOR OF DIMENSION N CONTAINING THE GRADIENT
  29. C                    VECTOR CORRESPONDING TO THE MINIMUM ON RETURN,
  30. C                    I.E. G=G(X).
  31. C           EST    - IS AN ESTIMATE OF THE MINIMUM FUNCTION VALUE.
  32. C           EPS    - TESTVALUE REPRESENTING THE EXPECTED ABSOLUTE ERROR.
  33. C                    A REASONABLE CHOICE IS 10**(-6), I.E.
  34. C                    SOMEWHAT GREATER THAN 10**(-D), WHERE D IS THE
  35. C                    NUMBER OF SIGNIFICANT DIGITS IN FLOATING POINT
  36. C                    REPRESENTATION.
  37. C           LIMIT  - MAXIMUM NUMBER OF ITERATIONS.
  38. C           IER    - ERROR PARAMETER
  39. C                    IER = 0 MEANS CONVERGENCE WAS OBTAINED
  40. C                    IER = 1 MEANS NO CONVERGENCE IN LIMIT ITERATIONS
  41. C                    IER =-1 MEANS ERRORS IN GRADIENT CALCULATION
  42. C                    IER = 2 MEANS LINEAR SEARCH TECHNIQUE INDICATES
  43. C                    IT IS LIKELY THAT THERE EXISTS NO MINIMUM.
  44. C           H      - WORKING STORAGE OF DIMENSION N*(N+7)/2.
  45. C
  46. C        REMARKS
  47. C            I) THE SUBROUTINE NAME REPLACING THE DUMMY ARGUMENT  FUNCT
  48. C               MUST BE DECLARED AS EXTERNAL IN THE CALLING PROGRAM.
  49. C           II) IER IS SET TO 2 IF , STEPPING IN ONE OF THE COMPUTED
  50. C               DIRECTIONS, THE FUNCTION WILL NEVER INCREASE WITHIN
  51. C               A TOLERABLE RANGE OF ARGUMENT.
  52. C               IER = 2 MAY OCCUR ALSO IF THE INTERVAL WHERE F
  53. C               INCREASES IS SMALL AND THE INITIAL ARGUMENT WAS
  54. C               RELATIVELY FAR AWAY FROM THE MINIMUM SUCH THAT THE
  55. C               MINIMUM WAS OVERLEAPED. THIS IS DUE TO THE SEARCH
  56. C               TECHNIQUE WHICH DOUBLES THE STEPSIZE UNTIL A POINT
  57. C               IS FOUND WHERE THE FUNCTION INCREASES.
  58. C
  59. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  60. C           FUNCT
  61. C
  62. C        METHOD
  63. C           THE METHOD IS DESCRIBED IN THE FOLLOWING ARTICLE
  64. C           R. FLETCHER AND M.J.D. POWELL, A RAPID DESCENT METHOD FOR
  65. C           MINIMIZATION,
  66. C           COMPUTER JOURNAL VOL.6, ISS. 2, 1963, PP.163-168.
  67. C
  68. C     ..................................................................
  69. C
  70.       SUBROUTINE FMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)
  71. C
  72. C        DIMENSIONED DUMMY VARIABLES
  73.       DIMENSION H(1),X(1),G(1)
  74. C
  75. C        COMPUTE FUNCTION VALUE AND GRADIENT VECTOR FOR INITIAL ARGUMENT
  76.       CALL FUNCT(N,X,F,G)
  77. C
  78. C        RESET ITERATION COUNTER AND GENERATE IDENTITY MATRIX
  79.       IER=0
  80.       KOUNT=0
  81.       N2=N+N
  82.       N3=N2+N
  83.       N31=N3+1
  84.     1 K=N31
  85.       DO 4 J=1,N
  86.       H(K)=1.
  87.       NJ=N-J
  88.       IF(NJ)5,5,2
  89.     2 DO 3 L=1,NJ
  90.       KL=K+L
  91.     3 H(KL)=0.
  92.     4 K=KL+1
  93. C
  94. C        START ITERATION LOOP
  95.     5 KOUNT=KOUNT +1
  96. C
  97. C        SAVE FUNCTION VALUE, ARGUMENT VECTOR AND GRADIENT VECTOR
  98.       OLDF=F
  99.       DO 9 J=1,N
  100.       K=N+J
  101.       H(K)=G(J)
  102.       K=K+N
  103.       H(K)=X(J)
  104. C
  105. C        DETERMINE DIRECTION VECTOR H
  106.       K=J+N3
  107.       T=0.
  108.       DO 8 L=1,N
  109.       T=T-G(L)*H(K)
  110.       IF(L-J)6,7,7
  111.     6 K=K+N-L
  112.       GO TO 8
  113.     7 K=K+1
  114.     8 CONTINUE
  115.     9 H(J)=T
  116. C
  117. C        CHECK WHETHER FUNCTION WILL DECREASE STEPPING ALONG H.
  118.       DY=0.
  119.       HNRM=0.
  120.       GNRM=0.
  121. C
  122. C        CALCULATE DIRECTIONAL DERIVATIVE AND TESTVALUES FOR DIRECTION
  123. C        VECTOR H AND GRADIENT VECTOR G.
  124.       DO 10 J=1,N
  125.       HNRM=HNRM+ABS(H(J))
  126.       GNRM=GNRM+ABS(G(J))
  127.    10 DY=DY+H(J)*G(J)
  128. C
  129. C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTIONAL
  130. C        DERIVATIVE APPEARS TO BE POSITIVE OR ZERO.
  131.       IF(DY)11,51,51
  132. C
  133. C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTION
  134. C        VECTOR H IS SMALL COMPARED TO GRADIENT VECTOR G.
  135.    11 IF(HNRM/GNRM-EPS)51,51,12
  136. C
  137. C        SEARCH MINIMUM ALONG DIRECTION H
  138. C
  139. C        SEARCH ALONG H FOR POSITIVE DIRECTIONAL DERIVATIVE
  140.    12 FY=F
  141.       ALFA=2.*(EST-F)/DY
  142.       AMBDA=1.
  143. C
  144. C        USE ESTIMATE FOR STEPSIZE ONLY IF IT IS POSITIVE AND LESS THAN
  145. C        1. OTHERWISE TAKE 1. AS STEPSIZE
  146.       IF(ALFA)15,15,13
  147.    13 IF(ALFA-AMBDA)14,15,15
  148.    14 AMBDA=ALFA
  149.    15 ALFA=0.
  150. C
  151. C        SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT
  152.    16 FX=FY
  153.       DX=DY
  154. C
  155. C        STEP ARGUMENT ALONG H
  156.       DO 17 I=1,N
  157.    17 X(I)=X(I)+AMBDA*H(I)
  158. C
  159. C        COMPUTE FUNCTION VALUE AND GRADIENT FOR NEW ARGUMENT
  160.       CALL FUNCT(N,X,F,G)
  161.       FY=F
  162. C
  163. C        COMPUTE DIRECTIONAL DERIVATIVE DY FOR NEW ARGUMENT.  TERMINATE
  164. C        SEARCH, IF DY IS POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND
  165.       DY=0.
  166.       DO 18 I=1,N
  167.    18 DY=DY+G(I)*H(I)
  168.       IF(DY)19,36,22
  169. C
  170. C        TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT
  171. C        A MINIMUM HAS BEEN PASSED
  172.    19 IF(FY-FX)20,22,22
  173. C
  174. C        REPEAT SEARCH AND DOUBLE STEPSIZE FOR FURTHER SEARCHES
  175.    20 AMBDA=AMBDA+ALFA
  176.       ALFA=AMBDA
  177. C        END OF SEARCH LOOP
  178. C
  179. C        TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE
  180.       IF(HNRM*AMBDA-1.E10)16,16,21
  181. C
  182. C        LINEAR SEARCH TECHNIQUE INDICATES THAT NO MINIMUM EXISTS
  183.    21 IER=2
  184.       RETURN
  185. C
  186. C        INTERPOLATE CUBICALLY IN THE INTERVAL DEFINED BY THE SEARCH
  187. C        ABOVE AND COMPUTE THE ARGUMENT X FOR WHICH THE INTERPOLATION
  188. C        POLYNOMIAL IS MINIMIZED
  189.    22 T=0.
  190.    23 IF(AMBDA)24,36,24
  191.    24 Z=3.*(FX-FY)/AMBDA+DX+DY
  192.       ALFA=AMAX1(ABS(Z),ABS(DX),ABS(DY))
  193.       DALFA=Z/ALFA
  194.       DALFA=DALFA*DALFA-DX/ALFA*DY/ALFA
  195.       IF(DALFA)51,25,25
  196.    25 W=ALFA*SQRT(DALFA)
  197.       ALFA=DY-DX+W+W
  198.       IF(ALFA) 250,251,250
  199.   250 ALFA=(DY-Z+W)/ALFA
  200.       GO TO 252
  201.   251 ALFA=(Z+DY-W)/(Z+DX+Z+DY)
  202.   252 ALFA=ALFA*AMBDA
  203.       DO 26 I=1,N
  204.    26 X(I)=X(I)+(T-ALFA)*H(I)
  205. C
  206. C        TERMINATE, IF THE VALUE OF THE ACTUAL FUNCTION AT X IS LESS
  207. C        THAN THE FUNCTION VALUES AT THE INTERVAL ENDS. OTHERWISE REDUCE
  208. C        THE INTERVAL BY CHOOSING ONE END-POINT EQUAL TO X AND REPEAT
  209. C        THE INTERPOLATION.  WHICH END-POINT IS CHOOSEN DEPENDS ON THE
  210. C        VALUE OF THE FUNCTION AND ITS GRADIENT AT X
  211. C
  212.       CALL FUNCT(N,X,F,G)
  213.       IF(F-FX)27,27,28
  214.    27 IF(F-FY)36,36,28
  215.    28 DALFA=0.
  216.       DO 29 I=1,N
  217.    29 DALFA=DALFA+G(I)*H(I)
  218.       IF(DALFA)30,33,33
  219.    30 IF(F-FX)32,31,33
  220.    31 IF(DX-DALFA)32,36,32
  221.    32 FX=F
  222.       DX=DALFA
  223.       T=ALFA
  224.       AMBDA=ALFA
  225.       GO TO 23
  226.    33 IF(FY-F)35,34,35
  227.    34 IF(DY-DALFA)35,36,35
  228.    35 FY=F
  229.       DY=DALFA
  230.       AMBDA=AMBDA-ALFA
  231.       GO TO 22
  232. C
  233. C        TERMINATE, IF FUNCTION HAS NOT DECREASED DURING LAST ITERATION
  234.    36 IF(OLDF-F+EPS)51,38,38
  235. C
  236. C        COMPUTE DIFFERENCE VECTORS OF ARGUMENT AND GRADIENT FROM
  237. C        TWO CONSECUTIVE ITERATIONS
  238.    38 DO 37 J=1,N
  239.       K=N+J
  240.       H(K)=G(J)-H(K)
  241.       K=N+K
  242.    37 H(K)=X(J)-H(K)
  243. C
  244. C        TEST LENGTH OF ARGUMENT DIFFERENCE VECTOR AND DIRECTION VECTOR
  245. C        IF AT LEAST N ITERATIONS HAVE BEEN EXECUTED. TERMINATE, IF
  246. C        BOTH ARE LESS THAN  EPS
  247.       IER=0
  248.       IF(KOUNT-N)42,39,39
  249.    39 T=0.
  250.       Z=0.
  251.       DO 40 J=1,N
  252.       K=N+J
  253.       W=H(K)
  254.       K=K+N
  255.       T=T+ABS(H(K))
  256.    40 Z=Z+W*H(K)
  257.       IF(HNRM-EPS)41,41,42
  258.    41 IF(T-EPS)56,56,42
  259. C
  260. C        TERMINATE, IF NUMBER OF ITERATIONS WOULD EXCEED  LIMIT
  261.    42 IF(KOUNT-LIMIT)43,50,50
  262. C
  263. C        PREPARE UPDATING OF MATRIX H
  264.    43 ALFA=0.
  265.       DO 47 J=1,N
  266.       K=J+N3
  267.       W=0.
  268.       DO 46 L=1,N
  269.       KL=N+L
  270.       W=W+H(KL)*H(K)
  271.       IF(L-J)44,45,45
  272.    44 K=K+N-L
  273.       GO TO 46
  274.    45 K=K+1
  275.    46 CONTINUE
  276.       K=N+J
  277.       ALFA=ALFA+W*H(K)
  278.    47 H(J)=W
  279. C
  280. C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF RESULTS
  281. C        ARE NOT SATISFACTORY
  282.       IF(Z*ALFA)48,1,48
  283. C
  284. C        UPDATE MATRIX H
  285.    48 K=N31
  286.       DO 49 L=1,N
  287.       KL=N2+L
  288.       DO 49 J=L,N
  289.       NJ=N2+J
  290.       H(K)=H(K)+H(KL)*H(NJ)/Z-H(L)*H(J)/ALFA
  291.    49 K=K+1
  292.       GO TO 5
  293. C        END OF ITERATION LOOP
  294. C
  295. C        NO CONVERGENCE AFTER  LIMIT  ITERATIONS
  296.    50 IER=1
  297.       RETURN
  298. C
  299. C        RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS
  300.    51 DO 52 J=1,N
  301.       K=N2+J
  302.    52 X(J)=H(K)
  303.       CALL FUNCT(N,X,F,G)
  304. C
  305. C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DERIVATIVE
  306. C        FAILS TO BE SUFFICIENTLY SMALL
  307.       IF(GNRM-EPS)55,55,53
  308. C
  309. C        TEST FOR REPEATED FAILURE OF ITERATION
  310.    53 IF(IER)56,54,54
  311.    54 IER=-1
  312.       GOTO 1
  313.    55 IER=0
  314.    56 RETURN
  315.       END
  316.