C
C     ..................................................................
C
C        SUBROUTINE DFMCG
C
C        PURPOSE
C           TO FIND A LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES
C           BY THE METHOD OF CONJUGATE GRADIENTS
C
C        USAGE
C           CALL DFMCG(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 2*N.
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 C.M.REEVES, FUNCTION MINIMIZATION BY
C           CONJUGATE GRADIENTS,
C           COMPUTER JOURNAL VOL.7, ISS.2, 1964, PP.149-154.
C
C     ..................................................................
C
      SUBROUTINE DFMCG(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)
C
C        DIMENSIONED DUMMY VARIABLES
      DIMENSION X(1),G(1),H(1)
      DOUBLE PRECISION X,G,GNRM,H,HNRM,F,FX,FY,OLDF,OLDG,SNRM,AMBDA,
     1ALFA,DALFA,T,Z,W,DX,DY
C
C        COMPUTE FUNCTION VALUE AND GRADIENT VECTOR FOR INITIAL ARGUMENT
      CALL FUNCT(N,X,F,G)
C
C        RESET ITERATION COUNTER
      KOUNT=0
      IER=0
      N1=N+1
C
C        START ITERATION CYCLE FOR EVERY N+1 ITERATIONS
    1 DO 43 II=1,N1
C
C        STEP ITERATION COUNTER AND SAVE FUNCTION VALUE
      KOUNT=KOUNT+1
      OLDF=F
C
C        COMPUTE SQUARE OF GRADIENT AND TERMINATE IF ZERO
      GNRM=0.D0
      DO 2 J=1,N
    2 GNRM=GNRM+G(J)*G(J)
      IF(GNRM)46,46,3
C
C        EACH TIME THE ITERATION LOOP IS EXECUTED , THE FIRST STEP WILL
C        BE IN DIRECTION OF STEEPEST DESCENT
    3 IF(II-1)4,4,6
    4 DO 5 J=1,N
    5 H(J)=-G(J)
      GO TO 8
C
C        FURTHER DIRECTION VECTORS H WILL BE CHOOSEN CORRESPONDING
C        TO THE CONJUGATE GRADIENT METHOD
    6 AMBDA=GNRM/OLDG
      DO 7 J=1,N
    7 H(J)=AMBDA*H(J)-G(J)
C
C        COMPUTE TESTVALUE FOR DIRECTIONAL VECTOR AND DIRECTIONAL
C        DERIVATIVE
    8 DY=0.D0
      HNRM=0.D0
      DO 9 J=1,N
      K=J+N
C
C        SAVE ARGUMENT VECTOR
      H(K)=X(J)
      HNRM=HNRM+DABS(H(J))
    9 DY=DY+H(J)*G(J)
C
C        CHECK WHETHER FUNCTION WILL DECREASE STEPPING ALONG H AND
C        SKIP LINEAR SEARCH ROUTINE IF NOT
      IF(DY)10,42,42
C
C        COMPUTE SCALE FACTOR USED IN LINEAR SEARCH SUBROUTINE
   10 SNRM=1.D0/HNRM
C
C        SEARCH MINIMUM ALONG DIRECTION H
C
C        SEARCH ALONG H FOR POSITIVE DIRECTIONAL DERIVATIVE
      FY=F
      ALFA=2.D0*(EST-F)/DY
      AMBDA=SNRM
C
C        USE ESTIMATE FOR STEPSIZE ONLY IF IT IS POSITIVE AND LESS THAN
C        SNRM. OTHERWISE TAKE SNRM AS STEPSIZE.
      IF(ALFA)13,13,11
   11 IF(ALFA-AMBDA)12,13,13
   12 AMBDA=ALFA
   13 ALFA=0.D0
C
C        SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT
   14 FX=FY
      DX=DY
C
C        STEP ARGUMENT ALONG H
      DO 15 I=1,N
   15 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 POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND
      DY=0.D0
      DO 16 I=1,N
   16 DY=DY+G(I)*H(I)
      IF(DY)17,38,20
C
C        TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT
C        A MINIMUM HAS BEEN PASSED
   17 IF(FY-FX)18,20,20
C
C        REPEAT SEARCH AND DOUBLE STEPSIZE FOR FURTHER SEARCHES
   18 AMBDA=AMBDA+ALFA
      ALFA=AMBDA
C
C        TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE
      IF(HNRM*AMBDA-1.D10)14,14,19
C
C        LINEAR SEARCH TECHNIQUE INDICATES THAT NO MINIMUM EXISTS
   19 IER=2
C
C        RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS
      F=OLDF
      DO 100 J=1,N
      G(J)=H(J)
      K=N+J
  100 X(J)=H(K)
      RETURN
C        END OF SEARCH LOOP
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
C
   20 T=0.
   21 IF(AMBDA)22,38,22
   22 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)23,27,27
C
C        RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS
   23 DO 24 J=1,N
      K=N+J
   24 X(J)=H(K)
      CALL FUNCT(N,X,F,G)
C
C        TEST FOR REPEATED FAILURE OF ITERATION
   25 IF(IER)47,26,47
   26 IER=-1
      GOTO 1
   27 W=ALFA*DSQRT(DALFA)
      ALFA=DY-DX+W+W
      IF(ALFA)270,271,270
  270 ALFA=(DY-Z+W)/ALFA
      GO TO 272
  271 ALFA=(Z+DY-W)/(Z+DX+Z+DY)
  272 ALFA=ALFA*AMBDA
      DO 28 I=1,N
   28 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)29,29,30
   29 IF(F-FY)38,38,30
C
C        COMPUTE DIRECTIONAL DERIVATIVE
   30 DALFA=0.D0
      DO 31 I=1,N
   31 DALFA=DALFA+G(I)*H(I)
      IF(DALFA)32,35,35
   32 IF(F-FX)34,33,35
   33 IF(DX-DALFA)34,38,34
   34 FX=F
      DX=DALFA
      T=ALFA
      AMBDA=ALFA
      GO TO 21
   35 IF(FY-F)37,36,37
   36 IF(DY-DALFA)37,38,37
   37 FY=F
      DY=DALFA
      AMBDA=AMBDA-ALFA
      GO TO 20
C
C        TERMINATE, IF FUNCTION HAS NOT DECREASED DURING LAST ITERATION
C        OTHERWISE SAVE GRADIENT NORM
   38 IF(OLDF-F+EPS)19,25,39
   39 OLDG=GNRM
C
C        COMPUTE DIFFERENCE OF NEW AND OLD ARGUMENT VECTOR
      T=0.D0
      DO 40 J=1,N
      K=J+N
      H(K)=X(J)-H(K)
   40 T=T+DABS(H(K))
C
C        TEST LENGTH OF DIFFERENCE VECTOR IF AT LEAST N+1 ITERATIONS
C        HAVE BEEN EXECUTED. TERMINATE, IF LENGTH IS LESS THAN EPS
      IF(KOUNT-N1)42,41,41
   41 IF(T-EPS)45,45,42
C
C        TERMINATE, IF NUMBER OF ITERATIONS WOULD EXCEED  LIMIT
   42 IF(KOUNT-LIMIT)43,44,44
   43 IER=0
C        END OF ITERATION CYCLE
C
C        START NEXT ITERATION CYCLE
      GO TO 1
C
C        NO CONVERGENCE AFTER  LIMIT  ITERATIONS
   44 IER=1
      IF(GNRM-EPS)46,46,47
C
C        TEST FOR SUFFICIENTLY SMALL GRADIENT
   45 IF(GNRM-EPS)46,46,25
   46 IER=0
   47 RETURN
      END
