home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / EXTRMFCT.ZIP / FMCG.FOR < prev    next >
Text File  |  1985-11-29  |  9KB  |  280 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE FMCG
  5. C
  6. C        PURPOSE
  7. C           TO FIND A LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES
  8. C           BY THE METHOD OF CONJUGATE GRADIENTS
  9. C
  10. C        USAGE
  11. C           CALL FMCG(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 2*N.
  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 C.M.REEVES, FUNCTION MINIMIZATION BY
  65. C           CONJUGATE GRADIENTS,
  66. C           COMPUTER JOURNAL VOL.7, ISS.2, 1964, PP.149-154.
  67. C
  68. C     ..................................................................
  69. C
  70.       SUBROUTINE FMCG(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)
  71. C
  72. C        DIMENSIONED DUMMY VARIABLES
  73.       DIMENSION X(1),G(1),H(1)
  74. C
  75. C
  76. C        COMPUTE FUNCTION VALUE AND GRADIENT VECTOR FOR INITIAL ARGUMENT
  77.       CALL FUNCT(N,X,F,G)
  78. C
  79. C        RESET ITERATION COUNTER
  80.       KOUNT=0
  81.       IER=0
  82.       N1=N+1
  83. C
  84. C        START ITERATION CYCLE FOR EVERY N+1 ITERATIONS
  85.     1 DO 43 II=1,N1
  86. C
  87. C        STEP ITERATION COUNTER AND SAVE FUNCTION VALUE
  88.       KOUNT=KOUNT+1
  89.       OLDF=F
  90. C
  91. C        COMPUTE SQUARE OF GRADIENT AND TERMINATE IF ZERO
  92.       GNRM=0.
  93.       DO 2 J=1,N
  94.     2 GNRM=GNRM+G(J)*G(J)
  95.       IF(GNRM)46,46,3
  96. C
  97. C        EACH TIME THE ITERATION LOOP IS EXECUTED , THE FIRST STEP WILL
  98. C        BE IN DIRECTION OF STEEPEST DESCENT
  99.     3 IF(II-1)4,4,6
  100.     4 DO 5 J=1,N
  101.     5 H(J)=-G(J)
  102.       GO TO 8
  103. C
  104. C        FURTHER DIRECTION VECTORS H WILL BE CHOOSEN CORRESPONDING
  105. C        TO THE CONJUGATE GRADIENT METHOD
  106.     6 AMBDA=GNRM/OLDG
  107.       DO 7 J=1,N
  108.     7 H(J)=AMBDA*H(J)-G(J)
  109. C
  110. C        COMPUTE TESTVALUE FOR DIRECTIONAL VECTOR AND DIRECTIONAL
  111. C        DERIVATIVE
  112.     8 DY=0.
  113.       HNRM=0.
  114.       DO 9 J=1,N
  115.       K=J+N
  116. C
  117. C        SAVE ARGUMENT VECTOR
  118.       H(K)=X(J)
  119.       HNRM=HNRM+ABS(H(J))
  120.     9 DY=DY+H(J)*G(J)
  121. C
  122. C        CHECK WHETHER FUNCTION WILL DECREASE STEPPING ALONG H AND
  123. C        SKIP LINEAR SEARCH ROUTINE IF NOT
  124.       IF(DY)10,42,42
  125. C
  126. C        COMPUTE SCALE FACTOR USED IN LINEAR SEARCH SUBROUTINE
  127.    10 SNRM=1./HNRM
  128. C
  129. C        SEARCH MINIMUM ALONG DIRECTION H
  130. C
  131. C        SEARCH ALONG H FOR POSITIVE DIRECTIONAL DERIVATIVE
  132.       FY=F
  133.       ALFA=2.*(EST-F)/DY
  134.       AMBDA=SNRM
  135. C
  136. C        USE ESTIMATE FOR STEPSIZE ONLY IF IT IS POSITIVE AND LESS THAN
  137. C        SNRM. OTHERWISE TAKE SNRM AS STEPSIZE.
  138.       IF(ALFA)13,13,11
  139.    11 IF(ALFA-AMBDA)12,13,13
  140.    12 AMBDA=ALFA
  141.    13 ALFA=0.
  142. C
  143. C        SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT
  144.    14 FX=FY
  145.       DX=DY
  146. C
  147. C        STEP ARGUMENT ALONG H
  148.       DO 15 I=1,N
  149.    15 X(I)=X(I)+AMBDA*H(I)
  150. C
  151. C        COMPUTE FUNCTION VALUE AND GRADIENT FOR NEW ARGUMENT
  152.       CALL FUNCT(N,X,F,G)
  153.       FY=F
  154. C
  155. C        COMPUTE DIRECTIONAL DERIVATIVE DY FOR NEW ARGUMENT.  TERMINATE
  156. C        SEARCH, IF DY POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND
  157.       DY=0.
  158.       DO 16 I=1,N
  159.    16 DY=DY+G(I)*H(I)
  160.       IF(DY)17,38,20
  161. C
  162. C        TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT
  163. C        A MINIMUM HAS BEEN PASSED
  164.    17 IF(FY-FX)18,20,20
  165. C
  166. C        REPEAT SEARCH AND DOUBLE STEPSIZE FOR FURTHER SEARCHES
  167.    18 AMBDA=AMBDA+ALFA
  168.       ALFA=AMBDA
  169. C
  170. C        TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE
  171.       IF(HNRM*AMBDA-1.E10)14,14,19
  172. C
  173. C        LINEAR SEARCH TECHNIQUE INDICATES THAT NO MINIMUM EXISTS
  174.    19 IER=2
  175. C
  176. C        RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS
  177.       F=OLDF
  178.       DO 100 J=1,N
  179.       G(J)=H(J)
  180.       K=N+J
  181.   100 X(J)=H(K)
  182.       RETURN
  183. C        END OF SEARCH LOOP
  184. C
  185. C        INTERPOLATE CUBICALLY IN THE INTERVAL DEFINED BY THE SEARCH
  186. C        ABOVE AND COMPUTE THE ARGUMENT X FOR WHICH THE INTERPOLATION
  187. C        POLYNOMIAL IS MINIMIZED
  188. C
  189.    20 T=0.
  190.    21 IF(AMBDA)22,38,22
  191.    22 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)23,27,27
  196. C
  197. C        RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS
  198.    23 DO 24 J=1,N
  199.       K=N+J
  200.    24 X(J)=H(K)
  201.       CALL FUNCT(N,X,F,G)
  202. C
  203. C        TEST FOR REPEATED FAILURE OF ITERATION
  204.    25 IF(IER)47,26,47
  205.    26 IER=-1
  206.       GOTO 1
  207.    27 W=ALFA*SQRT(DALFA)
  208.       ALFA=DY-DX+W+W
  209.       IF(ALFA)270,271,270
  210.   270 ALFA=(DY-Z+W)/ALFA
  211.       GO TO 272
  212.   271 ALFA=(Z+DY-W)/(Z+DX+Z+DY)
  213.   272 ALFA=ALFA*AMBDA
  214.       DO 28 I=1,N
  215.    28 X(I)=X(I)+(T-ALFA)*H(I)
  216. C
  217. C        TERMINATE, IF THE VALUE OF THE ACTUAL FUNCTION AT X IS LESS
  218. C        THAN THE FUNCTION VALUES AT THE INTERVAL ENDS. OTHERWISE REDUCE
  219. C        THE INTERVAL BY CHOOSING ONE END-POINT EQUAL TO X AND REPEAT
  220. C        THE INTERPOLATION.  WHICH END-POINT IS CHOOSEN DEPENDS ON THE
  221. C        VALUE OF THE FUNCTION AND ITS GRADIENT AT X
  222. C
  223.       CALL FUNCT(N,X,F,G)
  224.       IF(F-FX)29,29,30
  225.    29 IF(F-FY)38,38,30
  226. C
  227. C        COMPUTE DIRECTIONAL DERIVATIVE
  228.    30 DALFA=0.
  229.       DO 31 I=1,N
  230.    31 DALFA=DALFA+G(I)*H(I)
  231.       IF(DALFA)32,35,35
  232.    32 IF(F-FX)34,33,35
  233.    33 IF(DX-DALFA)34,38,34
  234.    34 FX=F
  235.       DX=DALFA
  236.       T=ALFA
  237.       AMBDA=ALFA
  238.       GO TO 21
  239.    35 IF(FY-F)37,36,37
  240.    36 IF(DY-DALFA)37,38,37
  241.    37 FY=F
  242.       DY=DALFA
  243.       AMBDA=AMBDA-ALFA
  244.       GO TO 20
  245. C
  246. C        TERMINATE, IF FUNCTION HAS NOT DECREASED DURING LAST ITERATION
  247. C        OTHERWISE SAVE GRADIENT NORM
  248.    38 IF(OLDF-F+EPS)19,25,39
  249.    39 OLDG=GNRM
  250. C
  251. C        COMPUTE DIFFERENCE OF NEW AND OLD ARGUMENT VECTOR
  252.       T=0.
  253.       DO 40 J=1,N
  254.       K=J+N
  255.       H(K)=X(J)-H(K)
  256.    40 T=T+ABS(H(K))
  257. C
  258. C        TEST LENGTH OF DIFFERENCE VECTOR IF AT LEAST N+1 ITERATIONS
  259. C        HAVE BEEN EXECUTED. TERMINATE, IF LENGTH IS LESS THAN EPS
  260.       IF(KOUNT-N1)42,41,41
  261.    41 IF(T-EPS)45,45,42
  262. C
  263. C        TERMINATE, IF NUMBER OF ITERATIONS WOULD EXCEED  LIMIT
  264.    42 IF(KOUNT-LIMIT)43,44,44
  265.    43 IER=0
  266. C        END OF ITERATION CYCLE
  267. C
  268. C        START NEXT ITERATION CYCLE
  269.       GO TO 1
  270. C
  271. C        NO CONVERGENCE AFTER  LIMIT  ITERATIONS
  272.    44 IER=1
  273.       IF(GNRM-EPS)46,46,47
  274. C
  275. C        TEST FOR SUFFICIENTLY SMALL GRADIENT
  276.    45 IF(GNRM-EPS)46,46,25
  277.    46 IER=0
  278.    47 RETURN
  279.       END
  280.