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