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