TABLE OF CONTENTS


ABINIT/m_lbfgs [ Modules ]

[ Top ] [ Modules ]

NAME

  m_lbfgs

FUNCTION

  This module provides several routines for the application of a
  Limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) minimization algorithm.
  The working routines were based on the original implementation of J. Nocera available on netlib.org
  They have been reshaped and translated into modern fortran here.

COPYRIGHT

 Copyright (C) 2012-2018 ABINIT group (FB)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 #include "abi_common.h"
29 
30 module m_lbfgs
31 
32  use defs_basis
33  use m_abicore
34 
35 type,public :: lbfgs_internal
36  integer              :: lbfgs_status
37  integer              :: ndim
38  integer              :: history_record
39  integer              :: iter
40  real(dp)             :: gtol
41  real(dp),allocatable :: diag(:)
42  real(dp),allocatable :: work(:)
43  real(dp)             :: line_stp
44  real(dp)             :: line_stpmin
45  real(dp)             :: line_stpmax
46  integer              :: line_info
47  integer              :: line_infoc
48  integer              :: line_nfev
49  real(dp)             :: line_dginit
50  real(dp)             :: line_finit
51  real(dp)             :: line_stx
52  real(dp)             :: line_fx
53  real(dp)             :: line_dgx
54  real(dp)             :: line_sty
55  real(dp)             :: line_fy
56  real(dp)             :: line_dgy
57  real(dp)             :: line_stmin
58  real(dp)             :: line_stmax
59  logical              :: line_bracket
60  logical              :: line_stage1
61 end type lbfgs_internal
62 
63 type(lbfgs_internal),save,public :: lbfgs_plan

m_lbfgs/lbfgs [ Functions ]

[ Top ] [ m_lbfgs ] [ Functions ]

NAME

 lbfgs

FUNCTION

   Perform the LBFGS step
   Fortran90 rewritting of the original subroutine by J. Nocera

INPUTS

OUTPUT

SIDE EFFECTS

PARENTS

      m_lbfgs

CHILDREN

SOURCE

247 subroutine lbfgs(N,M,X,F,G,DIAG,W,IFLAG,      &
248                  GTOL,STPMIN,STPMAX,STP,ITER, &
249                  INFO, NFEV,                  &
250                  LINE_DGINIT,LINE_FINIT,      &
251                  LINE_STX,LINE_FX,LINE_DGX,   &
252                  LINE_STY,LINE_FY,LINE_DGY,   &
253                  LINE_STMIN,LINE_STMAX,       &
254                  LINE_BRACKT,LINE_STAGE1,LINE_INFOC)
255 
256 
257 !This section has been created automatically by the script Abilint (TD).
258 !Do not modify the following lines by hand.
259 #undef ABI_FUNC
260 #define ABI_FUNC 'lbfgs'
261 !End of the abilint section
262 
263  implicit none
264 
265 !Arguments ------------------------------------
266 !scalars
267  integer,intent(inout) :: LINE_INFOC
268  integer,intent(inout)  :: ITER,IFLAG,INFO,NFEV
269  integer,intent(in)     :: N,M
270  real(dp),intent(inout) :: GTOL
271  real(dp),intent(in)    :: STPMIN,STPMAX
272  real(dp),intent(inout) :: STP
273  real(dp),intent(in)    :: F
274  real(dp),intent(inout) :: LINE_DGINIT,LINE_FINIT
275  real(dp),intent(inout) :: LINE_STX,LINE_FX,LINE_DGX
276  real(dp),intent(inout) :: LINE_STY,LINE_FY,LINE_DGY
277  real(dp),intent(inout) :: LINE_STMIN,LINE_STMAX
278  logical,intent(inout)  :: LINE_BRACKT,LINE_STAGE1
279 !arrays
280  real(dp),intent(inout) :: X(N),DIAG(N),W(N*(2*M+1)+2*M)
281  real(dp),intent(in)    :: G(N)
282 
283 !Local variables-------------------------------
284 !scalars
285  real(dp) :: FTOL,YS,YY,SQ,YR,BETA
286  integer :: POINT,ISPT,IYPT,MAXFEV, &
287          BOUND,NPT,CP,I,INMC,IYCN,ISCN
288 !***************************************************************************
289 
290 
291 !
292 ! Initialize
293 !-----------
294 
295 ! Parameters for line search routine
296  FTOL = 1.0D-4
297  MAXFEV = 20
298 
299  ISPT = N + 2 * M
300  IYPT = ISPT + N * M
301  POINT = MAX( 0 , MOD(ITER-1,M) )
302  NPT = POINT * N
303  ITER  = ITER + 1
304  BOUND = MIN( ITER-1 , M)
305 
306 
307  !
308  ! Entering the subroutine with a new position and gradient
309  ! or entering for the first time ever
310  if( IFLAG /= 1 ) then
311    W(ISPT+1:ISPT+N) = -G(1:N) * DIAG(1:N)
312 
313  else
314 
315    call MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,MAXFEV,INFO,NFEV, &
316                DIAG,GTOL,STPMIN,STPMAX,LINE_DGINIT,LINE_FINIT, &
317                LINE_STX,LINE_FX,LINE_DGX, &
318                LINE_STY,LINE_FY,LINE_DGY, &
319                LINE_STMIN,LINE_STMAX, &
320                LINE_BRACKT,LINE_STAGE1,LINE_INFOC)
321    !
322    ! Compute the new step and gradient change
323    !
324    NPT = POINT * N
325    W(ISPT+NPT+1:ISPT+NPT+N) = STP * W(ISPT+NPT+1:ISPT+NPT+N)
326    W(IYPT+NPT+1:IYPT+NPT+N) = G(1:N) - W(1:N)
327 
328    YS = DOT_PRODUCT( W(IYPT+NPT+1:IYPT+NPT+N) , W(ISPT+NPT+1:ISPT+NPT+N) )
329    YY = DOT_PRODUCT( W(IYPT+NPT+1:IYPT+NPT+N) , W(IYPT+NPT+1:IYPT+NPT+N) )
330    DIAG(1:N)= YS / YY
331 
332 !
333 !  COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980,
334 !  "Updating quasi-Newton matrices with limited storage",
335 !  Mathematics of Computation, Vol.24, No.151, pp. 773-782.
336 !  ---------------------------------------------------------
337 !
338    POINT = MODULO(ITER - 1,M)
339    CP = POINT
340    if (POINT == 0) CP = M
341 
342    W(N+CP) = one / YS
343    W(1:N)  = -G(1:N)
344 
345    CP = POINT
346    do I= 1,BOUND
347      CP = CP - 1
348      if (CP ==  -1) CP = M - 1
349      SQ = DOT_PRODUCT(W(ISPT+CP*N+1:ISPT+CP*N+N),W(1:N))
350      INMC = N + M + CP + 1
351      IYCN = IYPT + CP * N
352      W(INMC)= W(N+CP+1) * SQ
353      W(1:N) = W(1:N) - W(INMC) * W(IYCN+1:IYCN+N)
354    enddo
355 
356    W(1:N) = DIAG(1:N) * W(1:N)
357 
358    do I=1,BOUND
359      YR = DOT_PRODUCT(W(IYPT+CP*N+1:IYPT+CP*N+N),W(1:N))
360      BETA = W(N+CP+1) * YR
361      INMC = N + M + CP + 1
362      BETA = W(INMC) - BETA
363      ISCN = ISPT + CP * N
364      W(1:N) = W(1:N) + BETA * W(ISCN+1:ISCN+N)
365      CP = CP + 1
366      if (CP == M) CP = 0
367    enddo
368 
369 !
370 !  STORE THE NEW SEARCH DIRECTION
371    W(ISPT+POINT*N+1:ISPT+POINT*N+N) = W(1:N)
372 
373  endif
374 
375 !
376 ! Obtain the one-dimensional minimizer of the function
377 ! by using the line search routine mcsrch
378 !----------------------------------------------------
379  NFEV = 0
380  STP = one
381  W(1:N) = G(1:N)
382 
383  INFO  = 0
384 
385  call MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,MAXFEV,INFO,NFEV, &
386              DIAG,GTOL,STPMIN,STPMAX,LINE_DGINIT,LINE_FINIT, &
387              LINE_STX,LINE_FX,LINE_DGX, &
388              LINE_STY,LINE_FY,LINE_DGY, &
389              LINE_STMIN,LINE_STMAX, &
390              LINE_BRACKT,LINE_STAGE1,LINE_INFOC)
391 
392  if (INFO  ==  -1) then
393    IFLAG = 1
394    return
395  else
396    IFLAG = -1
397    return
398  endif
399 
400 end subroutine lbfgs

m_lbfgs/lbfgs_destroy [ Functions ]

[ Top ] [ m_lbfgs ] [ Functions ]

NAME

 lbfgs_destroy

FUNCTION

   Free the memory of the internal object lbfgs_internal for LBFGS minimization

INPUTS

OUTPUT

PARENTS

      pred_lbfgs

CHILDREN

SOURCE

148 subroutine lbfgs_destroy()
149 
150 
151 !This section has been created automatically by the script Abilint (TD).
152 !Do not modify the following lines by hand.
153 #undef ABI_FUNC
154 #define ABI_FUNC 'lbfgs_destroy'
155 !End of the abilint section
156 
157 implicit none
158 
159  if(allocated (lbfgs_plan%work)) then
160    ABI_DEALLOCATE(lbfgs_plan%work)
161  end if
162  if(allocated (lbfgs_plan%diag)) then
163    ABI_DEALLOCATE(lbfgs_plan%diag)
164  end if
165 
166 end subroutine lbfgs_destroy

m_lbfgs/lbfgs_execute [ Functions ]

[ Top ] [ m_lbfgs ] [ Functions ]

NAME

 lbfgs_execute

FUNCTION

   Perform one-step of LBFGS minimization
   all the internal information are stored in the lbfgs_internal object

INPUTS

   x: input and output position vector (atomic reduced coordinates + cell parameters)
   f: total energy
   gradf: gradient of the total energy (=negative forces)

OUTPUT

PARENTS

      pred_lbfgs

CHILDREN

SOURCE

193 function lbfgs_execute(x,f,gradf)
194 
195 
196 !This section has been created automatically by the script Abilint (TD).
197 !Do not modify the following lines by hand.
198 #undef ABI_FUNC
199 #define ABI_FUNC 'lbfgs_execute'
200 !End of the abilint section
201 
202 implicit none
203 real(dp),intent(inout) :: x(lbfgs_plan%ndim)
204 real(dp),intent(in)    :: f
205 real(dp),intent(in)    :: gradf(lbfgs_plan%ndim)
206 integer                :: lbfgs_execute
207 
208  call lbfgs(lbfgs_plan%ndim, lbfgs_plan%history_record, x, f, gradf, lbfgs_plan%diag, lbfgs_plan%work, lbfgs_plan%lbfgs_status, &
209        lbfgs_plan%gtol, lbfgs_plan%line_stpmin, lbfgs_plan%line_stpmax, lbfgs_plan%line_stp, lbfgs_plan%iter,                   &
210        lbfgs_plan%line_info, lbfgs_plan%line_nfev,                                                   &
211        lbfgs_plan%line_dginit, lbfgs_plan%line_finit,                                                &
212        lbfgs_plan%line_stx,  lbfgs_plan%line_fx,  lbfgs_plan%line_dgx,                               &
213        lbfgs_plan%line_sty,  lbfgs_plan%line_fy,  lbfgs_plan%line_dgy,                               &
214        lbfgs_plan%line_stmin,  lbfgs_plan%line_stmax,                                                &
215        lbfgs_plan%line_bracket, lbfgs_plan%line_stage1, lbfgs_plan%line_infoc)
216 
217 
218  lbfgs_execute = lbfgs_plan%lbfgs_status
219 
220 end function lbfgs_execute

m_lbfgs/lbfgs_init [ Functions ]

[ Top ] [ m_lbfgs ] [ Functions ]

NAME

 lbfgs_init

FUNCTION

   Initialize the internal object lbfgs_internal for LBFGS minimization

INPUTS

OUTPUT

SIDE EFFECTS

PARENTS

      pred_lbfgs

CHILDREN

SOURCE

 91 subroutine lbfgs_init(ndim,history_record,diag_guess)
 92 
 93 
 94 !This section has been created automatically by the script Abilint (TD).
 95 !Do not modify the following lines by hand.
 96 #undef ABI_FUNC
 97 #define ABI_FUNC 'lbfgs_init'
 98 !End of the abilint section
 99 
100 implicit none
101 
102 integer,intent(in)   :: ndim
103 integer,intent(in)   :: history_record
104 real(dp),intent(in)  :: diag_guess(ndim)
105 
106 integer :: nwork
107 
108  lbfgs_plan%lbfgs_status = 0
109  lbfgs_plan%iter   = 0
110  lbfgs_plan%ndim = ndim
111  lbfgs_plan%history_record = history_record
112  ABI_ALLOCATE(lbfgs_plan%diag,(ndim))
113 
114  nwork = ndim * ( 2 * history_record + 1 ) + 2 * history_record
115  ABI_ALLOCATE(lbfgs_plan%work,(nwork))
116 
117  lbfgs_plan%gtol = 0.9
118  lbfgs_plan%line_stpmin = 1.0e-20
119  lbfgs_plan%line_stpmax = 1.0e+20
120  lbfgs_plan%line_stp    = 1.0
121 
122  lbfgs_plan%diag(:) = diag_guess(:)
123 
124 end subroutine lbfgs_init

m_lbfgs/mcsrch [ Functions ]

[ Top ] [ m_lbfgs ] [ Functions ]

NAME

 mcsrch

FUNCTION

   Perform the line minimization step
   Fortran90 rewritting of the original subroutine by J. Nocera

INPUTS

OUTPUT

SIDE EFFECTS

PARENTS

      m_lbfgs

CHILDREN

SOURCE

426 subroutine mcsrch(N,X,F,G,S,STP,FTOL,MAXFEV,INFO,NFEV,WA, &
427                   GTOL,STPMIN,STPMAX,DGINIT,FINIT, &
428                   STX,FX,DGX,STY,FY,DGY,STMIN,STMAX, &
429                   BRACKT,STAGE1,INFOC)
430 
431 
432 !This section has been created automatically by the script Abilint (TD).
433 !Do not modify the following lines by hand.
434 #undef ABI_FUNC
435 #define ABI_FUNC 'mcsrch'
436 !End of the abilint section
437 
438  implicit none
439 
440 !Arguments ------------------------------------
441 !scalars
442  integer,intent(in)     :: N,MAXFEV
443  integer,intent(inout)  :: INFO,NFEV
444  integer,intent(inout)  :: INFOC
445  real(dp),intent(in)    :: GTOL,STPMIN,STPMAX
446  real(dp),intent(in)    :: F,FTOL
447  real(dp),intent(inout) :: STP,DGINIT,FINIT
448  real(dp),intent(inout) :: STX,FX,DGX
449  real(dp),intent(inout) :: STY,FY,DGY
450  real(dp),intent(inout) :: STMIN,STMAX
451  logical,intent(inout) :: BRACKT,STAGE1
452 !arrays
453  real(dp),intent(in)     :: G(N)
454  real(dp),intent(inout)  :: X(N),S(N),WA(N)
455 
456 !Local variables-------------------------------
457 !scalars
458  real(dp),parameter :: XTOL=1.0e-17_dp
459  real(dp),parameter :: P5     = 0.50_dp
460  real(dp),parameter :: P66    = 0.66_dp
461  real(dp),parameter :: XTRAPF = 4.00_dp
462  real(dp) :: DG,DGM,DGTEST,DGXM,DGYM, &
463         FTEST1,FM,FXM,FYM,WIDTH,WIDTH1
464 !***************************************************************************
465 
466  DGTEST = FTOL * DGINIT
467  WIDTH = STPMAX - STPMIN
468  WIDTH1 = WIDTH / P5
469 
470  ! Is it a first entry (info == 0)
471  ! or a second entry (info == -1)?
472  if( INFO == -1 ) then
473 
474    ! Reset INFO
475    INFO = 0
476 
477    NFEV = NFEV + 1
478    DG = SUM( G(:) * S(:) )
479    FTEST1 = FINIT + STP * DGTEST
480 !
481 !  TEST FOR CONVERGENCE.
482 !
483    if ((BRACKT .AND. (STP <= STMIN .OR. STP >= STMAX)) &
484       .OR. INFOC  ==  0) INFO = 6
485    if (STP  ==  STPMAX .AND. &
486        F <= FTEST1 .AND. DG <= DGTEST) INFO = 5
487    if (STP  ==  STPMIN .AND.  &
488        (F > FTEST1 .OR. DG >= DGTEST)) INFO = 4
489    if (NFEV >= MAXFEV) INFO = 3
490    if (BRACKT .AND. STMAX-STMIN <= XTOL*STMAX) INFO = 2
491    if (F <= FTEST1 .AND. ABS(DG) <= GTOL*(-DGINIT)) INFO = 1
492 !
493 !  CHECK FOR TERMINATION.
494 !
495    if (INFO /= 0) return
496 !
497 !  IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
498 !  FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
499 !
500    if (STAGE1 .AND. F <= FTEST1 .AND. &
501        DG >= MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE.
502 !
503 !  A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
504 !  WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
505 !  FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
506 !  DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
507 !  OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
508 !
509    if (STAGE1 .AND. F <= FX .AND. F > FTEST1) then
510 !
511 !     DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
512 !
513       FM = F - STP * DGTEST
514       FXM = FX - STX * DGTEST
515       FYM = FY - STY * DGTEST
516       DGM = DG - DGTEST
517       DGXM = DGX - DGTEST
518       DGYM = DGY - DGTEST
519 !
520 !     CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
521 !     AND TO COMPUTE THE NEW STEP.
522 !
523       call mcstep(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,BRACKT,STMIN,STMAX,INFOC)
524 !
525 !     RESET THE FUNCTION AND GRADIENT VALUES FOR F.
526 !
527       FX = FXM + STX * DGTEST
528       FY = FYM + STY * DGTEST
529       DGX = DGXM + DGTEST
530       DGY = DGYM + DGTEST
531    else
532 !
533 !     CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
534 !     AND TO COMPUTE THE NEW STEP.
535 !
536       call mcstep(STX,FX,DGX,STY,FY,DGY,STP,F,DG,BRACKT,STMIN,STMAX,INFOC)
537    endif
538 !
539 !  FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
540 !  INTERVAL OF UNCERTAINTY.
541 !
542    if (BRACKT) then
543       if (ABS(STY-STX) >= P66 * WIDTH1) STP = STX + P5 * (STY - STX)
544       WIDTH1 = WIDTH
545       WIDTH = ABS(STY-STX)
546    endif
547 
548  else
549 
550    INFOC = 1
551 !
552 !  CHECK THE INPUT PARAMETERS FOR ERRORS.
553 !
554    if ( STP <= zero .OR. FTOL < zero .OR.  &
555        GTOL < zero .OR. XTOL < zero .OR. STPMIN < zero &
556        .OR. STPMAX < STPMIN ) return
557 !
558 !  COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
559 !  AND CHECK THAT S IS A DESCENT DIRECTION.
560 !
561    DGINIT = DOT_PRODUCT( G , S )
562 
563    if (DGINIT > zero) then
564      return
565    endif
566 !
567 !  INITIALIZE LOCAL VARIABLES.
568 !
569 
570    BRACKT = .FALSE.
571    STAGE1 = .TRUE.
572    NFEV = 0
573    FINIT = F
574    DGTEST = FTOL * DGINIT
575    WA(:) = X(:)
576 
577 
578 !
579 !  THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
580 !  FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
581 !  THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
582 !  FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
583 !  THE INTERVAL OF UNCERTAINTY.
584 !  THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
585 !  FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
586 !
587    STX = zero
588    FX = FINIT
589    DGX = DGINIT
590    STY = zero
591    FY = FINIT
592    DGY = DGINIT
593  endif
594 
595 !
596 !SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
597 !TO THE PRESENT INTERVAL OF UNCERTAINTY.
598 !
599  if (BRACKT) then
600     STMIN = MIN(STX,STY)
601     STMAX = MAX(STX,STY)
602  else
603     STMIN = STX
604     STMAX = STP + XTRAPF*(STP - STX)
605  endif
606 !
607 !FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
608 !
609  STP = MAX(STPMIN,STP)
610  STP = MIN(STP,STPMAX)
611 !
612 !IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
613 !STP BE THE LOWEST POINT OBTAINED SO FAR.
614 !
615  if ((BRACKT .AND. (STP <= STMIN .OR. STP >= STMAX)) &
616     .OR. NFEV >= MAXFEV-1 .OR. INFOC  ==  0 &
617     .OR. (BRACKT .AND. STMAX-STMIN <= XTOL*STMAX)) STP = STX
618 
619 !
620 !Evaluate the function and gradient at STP
621 !and compute the directional derivative.
622 !We return to main program to obtain F and G.
623 !
624  X(:) = WA(:) + STP * S(:)
625 
626  INFO = -1
627 
628 end subroutine mcsrch

m_lbfgs/mcstep [ Functions ]

[ Top ] [ m_lbfgs ] [ Functions ]

NAME

 mcstep

FUNCTION

   Perform the step choice in line minimization
   Fortran90 rewritting of the original subroutine by J. Nocera

INPUTS

OUTPUT

SIDE EFFECTS

PARENTS

      m_lbfgs

CHILDREN

SOURCE

654 subroutine mcstep(STX,FX,DX,STY,FY,DY,STP,FP,DG,BRACKT,STPMIN,STPMAX,INFO)
655 
656 
657 !This section has been created automatically by the script Abilint (TD).
658 !Do not modify the following lines by hand.
659 #undef ABI_FUNC
660 #define ABI_FUNC 'mcstep'
661 !End of the abilint section
662 
663  implicit none
664 
665 !Arguments ------------------------------------
666 !scalars
667  integer,intent(inout)  :: INFO
668  real(dp),intent(in)     :: FP
669  real(dp),intent(inout)  :: STX,FX,DX,STY,FY,DY,STP,DG,STPMIN,STPMAX
670  logical,intent(inout) :: BRACKT
671 
672 !Local variables-------------------------------
673 !scalars
674  logical BOUND
675  real(dp) GAM,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA
676 !***************************************************************************
677 
678  INFO = 0
679 !
680 ! CHECK THE INPUT PARAMETERS FOR ERRORS.
681 !
682  IF ((BRACKT .AND. (STP <= MIN(STX,STY) .OR. &
683      STP >= MAX(STX,STY))) .OR.  &
684      DX*(STP-STX) >= 0.0 .OR. STPMAX < STPMIN) RETURN
685 !
686 ! Determine if the derivatives have opposite sign
687 !
688  SGND = DG * ( DX / ABS(DX) )
689 
690 ! FIRST CASE. A HIGHER FUNCTION VALUE.
691 ! THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
692 ! TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
693 ! ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
694 !
695  IF (FP > FX) THEN
696     INFO = 1
697     BOUND = .TRUE.
698     THETA = 3*(FX - FP)/(STP - STX) + DX + DG
699     S = MAX(ABS(THETA),ABS(DX),ABS(DG))
700     GAM = S * SQRT( (THETA/S)**2 - (DX/S)*(DG/S) )
701     IF (STP < STX) GAM = -GAM
702     P = (GAM - DX) + THETA
703     Q = ((GAM - DX) + GAM) + DG
704     R = P / Q
705     STPC = STX + R*(STP - STX)
706     STPQ = STX + ( ( DX / ( ( FX - FP ) / ( STP - STX ) + DX ) ) / 2 ) * ( STP - STX )
707     IF (ABS(STPC-STX) < ABS(STPQ-STX)) THEN
708        STPF = STPC
709     ELSE
710       STPF = STPC + (STPQ - STPC) / 2
711     END IF
712     BRACKT = .TRUE.
713 !
714 ! SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
715 ! OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
716 ! STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
717 ! THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
718 !
719  ELSE IF (SGND < 0.0) THEN
720     INFO = 2
721     BOUND = .FALSE.
722     THETA = 3*(FX - FP)/(STP - STX) + DX + DG
723     S = MAX(ABS(THETA),ABS(DX),ABS(DG))
724     GAM = S * SQRT( (THETA/S)**2 - (DX/S)*(DG/S) )
725     IF (STP > STX) GAM = -GAM
726     P = (GAM - DG) + THETA
727     Q = ((GAM - DG) + GAM) + DX
728     R = P/Q
729     STPC = STP + R*(STX - STP)
730     STPQ = STP + (DG/(DG-DX))*(STX - STP)
731     IF (ABS(STPC-STP) > ABS(STPQ-STP)) THEN
732        STPF = STPC
733     ELSE
734        STPF = STPQ
735     END IF
736     BRACKT = .TRUE.
737 !
738 ! THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
739 ! SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
740 ! THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
741 ! IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
742 ! IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
743 ! EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
744 ! COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
745 ! CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
746 !
747  ELSE IF (ABS(DG) < ABS(DX)) THEN
748     INFO = 3
749     BOUND = .TRUE.
750     THETA = 3*(FX - FP)/(STP - STX) + DX + DG
751     S = MAX(ABS(THETA),ABS(DX),ABS(DG))
752 !
753 !   THE CASE GAM = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
754 !   TO INFINITY IN THE DIRECTION OF THE STEP.
755 !
756     GAM = S * SQRT( MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DG/S)) )
757     IF (STP > STX) GAM = -GAM
758     P = (GAM - DG) + THETA
759     Q = (GAM + (DX - DG)) + GAM
760     R = P/Q
761     IF (R < 0.0 .AND. GAM .NE. 0.0) THEN
762        STPC = STP + R*(STX - STP)
763     ELSE IF (STP > STX) THEN
764        STPC = STPMAX
765     ELSE
766        STPC = STPMIN
767     END IF
768     STPQ = STP + (DG/(DG-DX))*(STX - STP)
769     IF (BRACKT) THEN
770        IF (ABS(STP-STPC) < ABS(STP-STPQ)) THEN
771           STPF = STPC
772        ELSE
773           STPF = STPQ
774        END IF
775     ELSE
776        IF (ABS(STP-STPC) > ABS(STP-STPQ)) THEN
777           STPF = STPC
778        ELSE
779           STPF = STPQ
780        END IF
781     END IF
782 !
783 ! FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
784 ! SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
785 ! NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
786 ! IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
787 !
788  ELSE
789     INFO = 4
790     BOUND = .FALSE.
791     IF (BRACKT) THEN
792        THETA = 3*(FP - FY)/(STY - STP) + DY + DG
793        S = MAX(ABS(THETA),ABS(DY),ABS(DG))
794        GAM = S * SQRT( (THETA/S)**2 - (DY/S)*(DG/S) )
795        IF (STP > STY) GAM = -GAM
796        P = (GAM - DG) + THETA
797        Q = ((GAM - DG) + GAM) + DY
798        R = P/Q
799        STPC = STP + R*(STY - STP)
800        STPF = STPC
801     ELSE IF (STP > STX) THEN
802        STPF = STPMAX
803     ELSE
804        STPF = STPMIN
805     END IF
806  END IF
807 
808 !
809 ! Update the interval of uncertainty. this update does not
810 ! depend on the new step or the case analysis above.
811 !
812  IF (FP > FX) THEN
813     STY = STP
814     FY = FP
815     DY = DG
816  ELSE
817     IF (SGND < 0.0) THEN
818        STY = STX
819        FY = FX
820        DY = DX
821     END IF
822     STX = STP
823     FX = FP
824     DX = DG
825  END IF
826 
827 !
828 ! Compute the new step and safeguard it.
829 !
830  STPF = MIN(STPMAX,STPF)
831  STPF = MAX(STPMIN,STPF)
832  STP = STPF
833  IF (BRACKT .AND. BOUND) THEN
834     IF (STY > STX) THEN
835        STP = MIN( STX + 0.66 * (STY-STX) , STP)
836     ELSE
837        STP = MAX( STX + 0.66 * (STY-STX) , STP)
838     END IF
839  END IF
840 
841 
842 end subroutine mcstep