mcsrch.f90

! This file is part of libmin
#define XTOL 1.0e-17
#define DEBUG .FALSE.
#define UNITDEBUG 101

!
!
!                 SUBROUTINE MCSRCH
!            
! A slight modification of the subroutine CSRCH of More' and Thuente.
! The changes are to allow reverse communication, and do not affect
! the performance of the routine. 
!
! THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES
! A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION.
!
! AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF
! UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF
! UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A
! MINIMIZER OF THE MODIFIED FUNCTION
!
!      F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
!
! IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION
! HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE,
! THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT
! CONTAINS A MINIMIZER OF F(X+STP*S).
!
! THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES
! THE SUFFICIENT DECREASE CONDITION
!
!       F(X+STP*S) <= F(X) + FTOL*STP*(GRADF(X)'S),
!
! AND THE CURVATURE CONDITION
!
!       ABS(GRADF(X+STP*S)'S)) <= GTOL*ABS(GRADF(X)'S).
!
! IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION
! IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES
! BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH
! CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING
! ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY
! SATISFIES THE SUFFICIENT DECREASE CONDITION.
!
! THE SUBROUTINE STATEMENT IS
!
!    SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA)
! WHERE
!
!   N IS A POSITIVE integer(C_INT) INPUT VARIABLE SET TO THE NUMBER
!     OF VARIABLES.
!
!   X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
!     BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS
!     X + STP*S.
!
!   F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F
!     AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S.
!
!   G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
!     GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT
!     OF F AT X + STP*S.
!
!   S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE
!     SEARCH DIRECTION.
!
!   STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN
!     INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT
!     STP CONTAINS THE FINAL ESTIMATE.
!
!   FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. (In this reverse
!     communication implementation GTOL is defined in a COMMON
!     statement.) TERMINATION OCCURS WHEN THE SUFFICIENT DECREASE
!     CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
!     SATISFIED.
!
!   XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
!     WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
!     IS AT MOST XTOL.
!
!   STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH
!     SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. (In this reverse
!     communication implementatin they are defined in a COMMON
!     statement).
!
!   MAXFEV IS A POSITIVE integer(C_INT) INPUT VARIABLE. TERMINATION
!     OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST
!     MAXFEV BY THE END OF AN ITERATION.
!
!   INFO IS AN integer(C_INT) OUTPUT VARIABLE SET AS FOLLOWS:
!
!     INFO = 0  IMPROPER INPUT PARAMETERS.
!
!     INFO =-1  A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT.
!
!     INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE
!               DIRECTIONAL DERIVATIVE CONDITION HOLD.
!
!     INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
!               IS AT MOST XTOL.
!
!     INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
!
!     INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN.
!
!     INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX.
!
!     INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.
!               THERE MAY NOT BE A STEP WHICH SATISFIES THE
!               SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
!               TOLERANCES MAY BE TOO SMALL.
!
!   NFEV IS AN integer(C_INT) OUTPUT VARIABLE SET TO THE NUMBER OF
!     CALLS TO FCN.
!
!   WA IS A WORK ARRAY OF LENGTH N.
!
! SUBPROGRAMS CALLED
!
!   MCSTEP
!
!   FORTRAN-SUPPLIED...ABS,MAX,MIN
!
! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
! JORGE J. MORE', DAVID J. THUENTE
!
! **********
subroutine mcsrch(N,X,F,G,S,STP,FTOL,MAXFEV,INFO,NFEV,WA, &
                  GTOL,STPMIN,STPMAX,DGINIT,FINIT, &
                  STX,FX,DGX,STY,FY,DGY,STMIN,STMAX, &
                  BRACKT,STAGE1,INFOC)
 use, intrinsic :: iso_c_binding
 use, intrinsic :: iso_fortran_env, only : ERROR_UNIT
 implicit none

 real(C_DOUBLE),intent(in)     :: GTOL,STPMIN,STPMAX
 integer(C_INT),intent(in)     :: N,MAXFEV
 integer,intent(inout)         :: INFO,NFEV
 real(C_DOUBLE),intent(in)     :: F,FTOL
 real(C_DOUBLE),intent(inout)  :: STP,DGINIT,FINIT
 real(C_DOUBLE),intent(in)     :: G(N)
 real(C_DOUBLE),intent(inout)  :: X(N),S(N),WA(N)
 logical(C_BOOL),intent(inout) :: BRACKT,STAGE1
 integer(C_INT),intent(inout)  :: INFOC
 real(C_DOUBLE),intent(inout)  :: STX,FX,DGX
 real(C_DOUBLE),intent(inout)  :: STY,FY,DGY
 real(C_DOUBLE),intent(inout)  :: STMIN,STMAX
!==================================================
 real(C_DOUBLE),parameter :: P5     = 0.50_C_DOUBLE
 real(C_DOUBLE),parameter :: P66    = 0.66_C_DOUBLE
 real(C_DOUBLE),parameter :: XTRAPF = 4.00_C_DOUBLE
 real(C_DOUBLE),parameter :: ZERO   = 0.00_C_DOUBLE
 integer(C_INT) :: J
 real(C_DOUBLE) :: DG,DGM,DGTEST,DGXM,DGYM, &
        FTEST1,FM,FXM,FYM,WIDTH,WIDTH1
!==================================================

 if( DEBUG ) then
   write(UNITDEBUG,*) 'INFO INFOC',INFO,INFOC
   write(UNITDEBUG,*) 'STP       ',STP
   write(UNITDEBUG,*) 'F,FTOL:',F,FTOL
   write(UNITDEBUG,*) 'X:',X(:)
   write(UNITDEBUG,*) 'G:',G(:)
   write(UNITDEBUG,*) 'S:',S(:)
   write(UNITDEBUG,*) 'WA:',WA(:)
 endif

 DGTEST = FTOL * DGINIT
 WIDTH = STPMAX - STPMIN
 WIDTH1 = WIDTH / P5

 ! Is it a first entry (info == 0)
 ! or a second entry (info == -1)?
 if( INFO == -1 ) then

   ! Reset INFO
   INFO = 0

   NFEV = NFEV + 1
   DG = SUM( G(:) * S(:) )
   FTEST1 = FINIT + STP * DGTEST
!
!  TEST FOR CONVERGENCE.
!
   if ((BRACKT .AND. (STP <= STMIN .OR. STP >= STMAX)) &
      .OR. INFOC  ==  0) INFO = 6
   if (STP  ==  STPMAX .AND. &
       F <= FTEST1 .AND. DG <= DGTEST) INFO = 5
   if (STP  ==  STPMIN .AND.  &
       (F > FTEST1 .OR. DG >= DGTEST)) INFO = 4
   if (NFEV >= MAXFEV) INFO = 3
   if (BRACKT .AND. STMAX-STMIN <= XTOL*STMAX) INFO = 2
   if (F <= FTEST1 .AND. ABS(DG) <= GTOL*(-DGINIT)) INFO = 1
!
!  CHECK FOR TERMINATION.
!
   if (INFO /= 0) return
!
!  IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
!  FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
!
   if (STAGE1 .AND. F <= FTEST1 .AND. &
       DG >= MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE.
!
!  A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
!  WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
!  FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
!  DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
!  OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
!
   if (STAGE1 .AND. F <= FX .AND. F > FTEST1) then
!
!     DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
!
      FM = F - STP * DGTEST
      FXM = FX - STX * DGTEST
      FYM = FY - STY * DGTEST
      DGM = DG - DGTEST
      DGXM = DGX - DGTEST
      DGYM = DGY - DGTEST
!
!     CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
!     AND TO COMPUTE THE NEW STEP.
!
      call mcstep(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,BRACKT,STMIN,STMAX,INFOC)
!
!     RESET THE FUNCTION AND GRADIENT VALUES FOR F.
!
      FX = FXM + STX * DGTEST
      FY = FYM + STY * DGTEST
      DGX = DGXM + DGTEST
      DGY = DGYM + DGTEST
   else
!
!     CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
!     AND TO COMPUTE THE NEW STEP.
!
      call mcstep(STX,FX,DGX,STY,FY,DGY,STP,F,DG,BRACKT,STMIN,STMAX,INFOC)
   endif
!
!  FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
!  INTERVAL OF UNCERTAINTY.
!
   if (BRACKT) then
      if (ABS(STY-STX) >= P66 * WIDTH1) STP = STX + P5 * (STY - STX)
      WIDTH1 = WIDTH
      WIDTH = ABS(STY-STX)
   endif

 else

   INFOC = 1
!
!  CHECK THE INPUT PARAMETERS FOR ERRORS.
!
   if ( STP <= ZERO .OR. FTOL < ZERO .OR.  &
       GTOL < ZERO .OR. XTOL < ZERO .OR. STPMIN < ZERO &
       .OR. STPMAX < STPMIN ) return
!
!  COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
!  AND CHECK THAT S IS A DESCENT DIRECTION.
!
   DGINIT = DOT_PRODUCT( G , S )

   if (DGINIT > ZERO) then
     write(ERROR_UNIT,'(x,a)') 'SEARCH DIRECTION IS NOT A DESCENT'
     return
   endif
!
!  INITIALIZE LOCAL VARIABLES.
!

   BRACKT = .FALSE.
   STAGE1 = .TRUE.
   NFEV = 0
   FINIT = F
   DGTEST = FTOL * DGINIT
   WA(:) = X(:)


!
!  THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
!  FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
!  THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
!  FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
!  THE INTERVAL OF UNCERTAINTY.
!  THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
!  FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
!
   STX = ZERO
   FX = FINIT
   DGX = DGINIT
   STY = ZERO
   FY = FINIT
   DGY = DGINIT
 endif

!
!SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
!TO THE PRESENT INTERVAL OF UNCERTAINTY.
!
 if (BRACKT) then
    STMIN = MIN(STX,STY)
    STMAX = MAX(STX,STY)
 else
    STMIN = STX
    STMAX = STP + XTRAPF*(STP - STX)
 endif
!
!FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
!
 STP = MAX(STPMIN,STP)
 STP = MIN(STP,STPMAX)
!
!IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
!STP BE THE LOWEST POINT OBTAINED SO FAR.
!
 if ((BRACKT .AND. (STP <= STMIN .OR. STP >= STMAX)) &
    .OR. NFEV >= MAXFEV-1 .OR. INFOC  ==  0 &
    .OR. (BRACKT .AND. STMAX-STMIN <= XTOL*STMAX)) STP = STX

!
!Evaluate the function and gradient at STP
!and compute the directional derivative.
!We return to main program to obtain F and G.
!
 X(:) = WA(:) + STP * S(:)

 INFO = -1


end subroutine mcsrch



! SUBROUTINE MCSTEP
!
! THE PURPOSE OF MCSTEP IS TO COMPUTE A SAFEGUARDED STEP FOR
! A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR
! A MINIMIZER OF THE FUNCTION.
!
! THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION
! VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS
! ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE
! DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A
! MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY
! WITH ENDPOINTS STX AND STY.
!
! THE SUBROUTINE STATEMENT IS
!
!   SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
!                    STPMIN,STPMAX,INFO)
!
! WHERE
!
!   STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP,
!     THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED
!     SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION
!     OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE
!     SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY.
!
!   STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP,
!     THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF
!     THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE
!     UPDATED APPROPRIATELY.
!
!   STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP,
!     THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP.
!     IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE
!     BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP.
!
!   BRACKT IS A logical(C_BOOL) VARIABLE WHICH SPECIFIES IF A MINIMIZER
!     HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED
!     THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER
!     IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE.
!
!   STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER
!     AND UPPER BOUNDS FOR THE STEP.
!
!   INFO IS AN integer(C_INT) OUTPUT VARIABLE SET AS FOLLOWS:
!     IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED
!     ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE
!     INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS.
!
!
! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
! JORGE J. MORE', DAVID J. THUENTE
!
subroutine mcstep(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,STPMIN,STPMAX,INFO)
 use, intrinsic :: iso_c_binding
 implicit none

 integer(C_INT),intent(inout)  :: INFO
 real(C_DOUBLE),intent(in)     :: FP
 real(C_DOUBLE),intent(inout)  :: STX,FX,DX,STY,FY,DY,STP,DP,STPMIN,STPMAX
 logical(C_BOOL),intent(inout) :: BRACKT
!==================================
 logical(C_BOOL) BOUND
 real(C_DOUBLE) GAM,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA
!==================================

 INFO = 0
!
! CHECK THE INPUT PARAMETERS FOR ERRORS.
!
 IF ((BRACKT .AND. (STP <= MIN(STX,STY) .OR. &
     STP >= MAX(STX,STY))) .OR.  &
     DX*(STP-STX) >= 0.0 .OR. STPMAX < STPMIN) RETURN
!
! Determine if the derivatives have opposite sign
!
 SGND = DP * ( DX / ABS(DX) )

! FIRST CASE. A HIGHER FUNCTION VALUE.
! THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
! TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
! ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
!
 IF (FP > FX) THEN
    INFO = 1
    BOUND = .TRUE.
    THETA = 3*(FX - FP)/(STP - STX) + DX + DP
    S = MAX(ABS(THETA),ABS(DX),ABS(DP))
    GAM = S * SQRT( (THETA/S)**2 - (DX/S)*(DP/S) )
    IF (STP < STX) GAM = -GAM
    P = (GAM - DX) + THETA
    Q = ((GAM - DX) + GAM) + DP
    R = P / Q
    STPC = STX + R*(STP - STX)
    STPQ = STX + ( ( DX / ( ( FX - FP ) / ( STP - STX ) + DX ) ) / 2 ) * ( STP - STX )
    IF (ABS(STPC-STX) < ABS(STPQ-STX)) THEN
       STPF = STPC
    ELSE
      STPF = STPC + (STPQ - STPC) / 2
    END IF
    BRACKT = .TRUE.
!
! SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
! OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
! STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
! THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
!
 ELSE IF (SGND < 0.0) THEN
    INFO = 2
    BOUND = .FALSE.
    THETA = 3*(FX - FP)/(STP - STX) + DX + DP
    S = MAX(ABS(THETA),ABS(DX),ABS(DP))
    GAM = S * SQRT( (THETA/S)**2 - (DX/S)*(DP/S) )
    IF (STP > STX) GAM = -GAM
    P = (GAM - DP) + THETA
    Q = ((GAM - DP) + GAM) + DX
    R = P/Q
    STPC = STP + R*(STX - STP)
    STPQ = STP + (DP/(DP-DX))*(STX - STP)
    IF (ABS(STPC-STP) > ABS(STPQ-STP)) THEN
       STPF = STPC
    ELSE
       STPF = STPQ
    END IF
    BRACKT = .TRUE.
!
! THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
! SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
! THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
! IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
! IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
! EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
! COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
! CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
!
 ELSE IF (ABS(DP) < ABS(DX)) THEN
    INFO = 3
    BOUND = .TRUE.
    THETA = 3*(FX - FP)/(STP - STX) + DX + DP
    S = MAX(ABS(THETA),ABS(DX),ABS(DP))
!
!   THE CASE GAM = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
!   TO INFINITY IN THE DIRECTION OF THE STEP.
!
    GAM = S * SQRT( MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S)) )
    IF (STP > STX) GAM = -GAM
    P = (GAM - DP) + THETA
    Q = (GAM + (DX - DP)) + GAM
    R = P/Q
    IF (R < 0.0 .AND. GAM .NE. 0.0) THEN
       STPC = STP + R*(STX - STP)
    ELSE IF (STP > STX) THEN
       STPC = STPMAX
    ELSE
       STPC = STPMIN
    END IF
    STPQ = STP + (DP/(DP-DX))*(STX - STP)
    IF (BRACKT) THEN
       IF (ABS(STP-STPC) < ABS(STP-STPQ)) THEN
          STPF = STPC
       ELSE
          STPF = STPQ
       END IF
    ELSE
       IF (ABS(STP-STPC) > ABS(STP-STPQ)) THEN
          STPF = STPC
       ELSE
          STPF = STPQ
       END IF
    END IF
!
! FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
! SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
! NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
! IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
!
 ELSE
    INFO = 4
    BOUND = .FALSE.
    IF (BRACKT) THEN
       THETA = 3*(FP - FY)/(STY - STP) + DY + DP
       S = MAX(ABS(THETA),ABS(DY),ABS(DP))
       GAM = S * SQRT( (THETA/S)**2 - (DY/S)*(DP/S) )
       IF (STP > STY) GAM = -GAM
       P = (GAM - DP) + THETA
       Q = ((GAM - DP) + GAM) + DY
       R = P/Q
       STPC = STP + R*(STY - STP)
       STPF = STPC
    ELSE IF (STP > STX) THEN
       STPF = STPMAX
    ELSE
       STPF = STPMIN
    END IF
 END IF

!
! Update the interval of uncertainty. this update does not
! depend on the new step or the case analysis above.
!
 IF (FP > FX) THEN
    STY = STP
    FY = FP
    DY = DP
 ELSE
    IF (SGND < 0.0) THEN
       STY = STX
       FY = FX
       DY = DX
    END IF
    STX = STP
    FX = FP
    DX = DP
 END IF

!
! Compute the new step and safeguard it.
!
 STPF = MIN(STPMAX,STPF)
 STPF = MAX(STPMIN,STPF)
 STP = STPF
 IF (BRACKT .AND. BOUND) THEN
    IF (STY > STX) THEN
       STP = MIN( STX + 0.66 * (STY-STX) , STP)
    ELSE
       STP = MAX( STX + 0.66 * (STY-STX) , STP)
    END IF
 END IF


end subroutine mcstep