TABLE OF CONTENTS
- ABINIT/m_lbfgs
- m_lbfgs/lbfgs
- m_lbfgs/lbfgs_destroy
- m_lbfgs/lbfgs_execute
- m_lbfgs/lbfgs_init
- m_lbfgs/mcsrch
- m_lbfgs/mcstep
ABINIT/m_lbfgs [ 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-2024 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 .
SOURCE
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 module m_lbfgs 26 27 use defs_basis 28 use m_abicore 29 30 implicit none 31 32 type,public :: lbfgs_internal 33 integer :: lbfgs_status 34 integer :: ndim 35 integer :: history_record 36 integer :: iter 37 real(dp) :: gtol 38 real(dp),allocatable :: diag(:) 39 real(dp),allocatable :: work(:) 40 real(dp) :: line_stp 41 real(dp) :: line_stpmin 42 real(dp) :: line_stpmax 43 integer :: line_info 44 integer :: line_infoc 45 integer :: line_nfev 46 real(dp) :: line_dginit 47 real(dp) :: line_finit 48 real(dp) :: line_stx 49 real(dp) :: line_fx 50 real(dp) :: line_dgx 51 real(dp) :: line_sty 52 real(dp) :: line_fy 53 real(dp) :: line_dgy 54 real(dp) :: line_stmin 55 real(dp) :: line_stmax 56 logical :: line_bracket 57 logical :: line_stage1 58 end type lbfgs_internal 59 60 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
SOURCE
195 subroutine lbfgs(N,M,X,F,G,DIAG,W,IFLAG, & 196 GTOL,STPMIN,STPMAX,STP,ITER, & 197 INFO, NFEV, & 198 LINE_DGINIT,LINE_FINIT, & 199 LINE_STX,LINE_FX,LINE_DGX, & 200 LINE_STY,LINE_FY,LINE_DGY, & 201 LINE_STMIN,LINE_STMAX, & 202 LINE_BRACKT,LINE_STAGE1,LINE_INFOC) 203 204 !Arguments ------------------------------------ 205 !scalars 206 integer,intent(inout) :: LINE_INFOC 207 integer,intent(inout) :: ITER,IFLAG,INFO,NFEV 208 integer,intent(in) :: N,M 209 real(dp),intent(inout) :: GTOL 210 real(dp),intent(in) :: STPMIN,STPMAX 211 real(dp),intent(inout) :: STP 212 real(dp),intent(in) :: F 213 real(dp),intent(inout) :: LINE_DGINIT,LINE_FINIT 214 real(dp),intent(inout) :: LINE_STX,LINE_FX,LINE_DGX 215 real(dp),intent(inout) :: LINE_STY,LINE_FY,LINE_DGY 216 real(dp),intent(inout) :: LINE_STMIN,LINE_STMAX 217 logical,intent(inout) :: LINE_BRACKT,LINE_STAGE1 218 !arrays 219 real(dp),intent(inout) :: X(N),DIAG(N),W(N*(2*M+1)+2*M) 220 real(dp),intent(in) :: G(N) 221 222 !Local variables------------------------------- 223 !scalars 224 real(dp) :: FTOL,YS,YY,SQ,YR,BETA 225 integer :: POINT,ISPT,IYPT,MAXFEV, & 226 BOUND,NPT,CP,I,INMC,IYCN,ISCN 227 !*************************************************************************** 228 229 230 ! 231 ! Initialize 232 !----------- 233 234 ! Parameters for line search routine 235 FTOL = 1.0D-4 236 MAXFEV = 20 237 238 ISPT = N + 2 * M 239 IYPT = ISPT + N * M 240 POINT = MAX( 0 , MOD(ITER-1,M) ) 241 NPT = POINT * N 242 ITER = ITER + 1 243 BOUND = MIN( ITER-1 , M) 244 245 246 ! 247 ! Entering the subroutine with a new position and gradient 248 ! or entering for the first time ever 249 if( IFLAG /= 1 ) then 250 W(ISPT+1:ISPT+N) = -G(1:N) * DIAG(1:N) 251 252 else 253 254 call MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,MAXFEV,INFO,NFEV, & 255 DIAG,GTOL,STPMIN,STPMAX,LINE_DGINIT,LINE_FINIT, & 256 LINE_STX,LINE_FX,LINE_DGX, & 257 LINE_STY,LINE_FY,LINE_DGY, & 258 LINE_STMIN,LINE_STMAX, & 259 LINE_BRACKT,LINE_STAGE1,LINE_INFOC) 260 ! 261 ! Compute the new step and gradient change 262 ! 263 NPT = POINT * N 264 W(ISPT+NPT+1:ISPT+NPT+N) = STP * W(ISPT+NPT+1:ISPT+NPT+N) 265 W(IYPT+NPT+1:IYPT+NPT+N) = G(1:N) - W(1:N) 266 267 YS = DOT_PRODUCT( W(IYPT+NPT+1:IYPT+NPT+N) , W(ISPT+NPT+1:ISPT+NPT+N) ) 268 YY = DOT_PRODUCT( W(IYPT+NPT+1:IYPT+NPT+N) , W(IYPT+NPT+1:IYPT+NPT+N) ) 269 DIAG(1:N)= YS / YY 270 271 ! 272 ! COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, 273 ! "Updating quasi-Newton matrices with limited storage", 274 ! Mathematics of Computation, Vol.24, No.151, pp. 773-782. 275 ! --------------------------------------------------------- 276 ! 277 POINT = MODULO(ITER - 1,M) 278 CP = POINT 279 if (POINT == 0) CP = M 280 281 W(N+CP) = one / YS 282 W(1:N) = -G(1:N) 283 284 CP = POINT 285 do I= 1,BOUND 286 CP = CP - 1 287 if (CP == -1) CP = M - 1 288 SQ = DOT_PRODUCT(W(ISPT+CP*N+1:ISPT+CP*N+N),W(1:N)) 289 INMC = N + M + CP + 1 290 IYCN = IYPT + CP * N 291 W(INMC)= W(N+CP+1) * SQ 292 W(1:N) = W(1:N) - W(INMC) * W(IYCN+1:IYCN+N) 293 enddo 294 295 W(1:N) = DIAG(1:N) * W(1:N) 296 297 do I=1,BOUND 298 YR = DOT_PRODUCT(W(IYPT+CP*N+1:IYPT+CP*N+N),W(1:N)) 299 BETA = W(N+CP+1) * YR 300 INMC = N + M + CP + 1 301 BETA = W(INMC) - BETA 302 ISCN = ISPT + CP * N 303 W(1:N) = W(1:N) + BETA * W(ISCN+1:ISCN+N) 304 CP = CP + 1 305 if (CP == M) CP = 0 306 enddo 307 308 ! 309 ! STORE THE NEW SEARCH DIRECTION 310 W(ISPT+POINT*N+1:ISPT+POINT*N+N) = W(1:N) 311 312 endif 313 314 ! 315 ! Obtain the one-dimensional minimizer of the function 316 ! by using the line search routine mcsrch 317 !---------------------------------------------------- 318 NFEV = 0 319 STP = one 320 W(1:N) = G(1:N) 321 322 INFO = 0 323 324 call MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,MAXFEV,INFO,NFEV, & 325 DIAG,GTOL,STPMIN,STPMAX,LINE_DGINIT,LINE_FINIT, & 326 LINE_STX,LINE_FX,LINE_DGX, & 327 LINE_STY,LINE_FY,LINE_DGY, & 328 LINE_STMIN,LINE_STMAX, & 329 LINE_BRACKT,LINE_STAGE1,LINE_INFOC) 330 331 if (INFO == -1) then 332 IFLAG = 1 333 return 334 else 335 IFLAG = -1 336 return 337 endif 338 339 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
SOURCE
126 subroutine lbfgs_destroy() 127 128 ABI_SFREE(lbfgs_plan%work) 129 ABI_SFREE(lbfgs_plan%diag) 130 131 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
SOURCE
153 function lbfgs_execute(x,f,gradf) 154 155 real(dp),intent(inout) :: x(lbfgs_plan%ndim) 156 real(dp),intent(in) :: f 157 real(dp),intent(in) :: gradf(lbfgs_plan%ndim) 158 integer :: lbfgs_execute 159 160 call lbfgs(lbfgs_plan%ndim, lbfgs_plan%history_record, x, f, gradf, lbfgs_plan%diag, lbfgs_plan%work, lbfgs_plan%lbfgs_status, & 161 lbfgs_plan%gtol, lbfgs_plan%line_stpmin, lbfgs_plan%line_stpmax, lbfgs_plan%line_stp, lbfgs_plan%iter, & 162 lbfgs_plan%line_info, lbfgs_plan%line_nfev, & 163 lbfgs_plan%line_dginit, lbfgs_plan%line_finit, & 164 lbfgs_plan%line_stx, lbfgs_plan%line_fx, lbfgs_plan%line_dgx, & 165 lbfgs_plan%line_sty, lbfgs_plan%line_fy, lbfgs_plan%line_dgy, & 166 lbfgs_plan%line_stmin, lbfgs_plan%line_stmax, & 167 lbfgs_plan%line_bracket, lbfgs_plan%line_stage1, lbfgs_plan%line_infoc) 168 169 170 !lbfgs_execute = lbfgs_plan%lbfgs_status 171 lbfgs_execute = lbfgs_plan%line_info 172 173 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
SOURCE
83 subroutine lbfgs_init(ndim,history_record,diag_guess) 84 85 integer,intent(in) :: ndim 86 integer,intent(in) :: history_record 87 real(dp),intent(in) :: diag_guess(ndim) 88 89 integer :: nwork 90 91 lbfgs_plan%lbfgs_status = 0 92 lbfgs_plan%iter = 0 93 lbfgs_plan%ndim = ndim 94 lbfgs_plan%history_record = history_record 95 ABI_MALLOC(lbfgs_plan%diag,(ndim)) 96 97 nwork = ndim * ( 2 * history_record + 1 ) + 2 * history_record 98 ABI_MALLOC(lbfgs_plan%work,(nwork)) 99 100 lbfgs_plan%gtol = 0.9 101 lbfgs_plan%line_stpmin = 1.0e-20 102 lbfgs_plan%line_stpmax = 1.0e+20 103 lbfgs_plan%line_stp = 1.0 104 105 lbfgs_plan%diag(:) = diag_guess(:) 106 107 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
SOURCE
360 subroutine mcsrch(N,X,F,G,S,STP,FTOL,MAXFEV,INFO,NFEV,WA, & 361 GTOL,STPMIN,STPMAX,DGINIT,FINIT, & 362 STX,FX,DGX,STY,FY,DGY,STMIN,STMAX, & 363 BRACKT,STAGE1,INFOC) 364 365 !Arguments ------------------------------------ 366 !scalars 367 integer,intent(in) :: N,MAXFEV 368 integer,intent(inout) :: INFO,NFEV 369 integer,intent(inout) :: INFOC 370 real(dp),intent(in) :: GTOL,STPMIN,STPMAX 371 real(dp),intent(in) :: F,FTOL 372 real(dp),intent(inout) :: STP,DGINIT,FINIT 373 real(dp),intent(inout) :: STX,FX,DGX 374 real(dp),intent(inout) :: STY,FY,DGY 375 real(dp),intent(inout) :: STMIN,STMAX 376 logical,intent(inout) :: BRACKT,STAGE1 377 !arrays 378 real(dp),intent(in) :: G(N) 379 real(dp),intent(inout) :: X(N),S(N),WA(N) 380 381 !Local variables------------------------------- 382 !scalars 383 real(dp),parameter :: XTOL=1.0e-17_dp 384 real(dp),parameter :: P5 = 0.50_dp 385 real(dp),parameter :: P66 = 0.66_dp 386 real(dp),parameter :: XTRAPF = 4.00_dp 387 real(dp) :: DG,DGM,DGTEST,DGXM,DGYM, & 388 FTEST1,FM,FXM,FYM,WIDTH,WIDTH1 389 !*************************************************************************** 390 391 DGTEST = FTOL * DGINIT 392 WIDTH = STPMAX - STPMIN 393 WIDTH1 = WIDTH / P5 394 395 ! Is it a first entry (info == 0) 396 ! or a second entry (info == -1)? 397 if( INFO == -1 ) then 398 399 ! Reset INFO 400 INFO = 0 401 402 NFEV = NFEV + 1 403 DG = SUM( G(:) * S(:) ) 404 FTEST1 = FINIT + STP * DGTEST 405 ! 406 ! TEST FOR CONVERGENCE. 407 ! 408 if ((BRACKT .AND. (STP <= STMIN .OR. STP >= STMAX)) & 409 .OR. INFOC == 0) INFO = 6 410 if (STP == STPMAX .AND. & 411 F <= FTEST1 .AND. DG <= DGTEST) INFO = 5 412 if (STP == STPMIN .AND. & 413 (F > FTEST1 .OR. DG >= DGTEST)) INFO = 4 414 if (NFEV >= MAXFEV) INFO = 3 415 if (BRACKT .AND. STMAX-STMIN <= XTOL*STMAX) INFO = 2 416 if (F <= FTEST1 .AND. ABS(DG) <= GTOL*(-DGINIT)) INFO = 1 417 ! 418 ! CHECK FOR TERMINATION. 419 ! 420 if (INFO /= 0) return 421 ! 422 ! IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED 423 ! FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE. 424 ! 425 if (STAGE1 .AND. F <= FTEST1 .AND. & 426 DG >= MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE. 427 ! 428 ! A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF 429 ! WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED 430 ! FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE 431 ! DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN 432 ! OBTAINED BUT THE DECREASE IS NOT SUFFICIENT. 433 ! 434 if (STAGE1 .AND. F <= FX .AND. F > FTEST1) then 435 ! 436 ! DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES. 437 ! 438 FM = F - STP * DGTEST 439 FXM = FX - STX * DGTEST 440 FYM = FY - STY * DGTEST 441 DGM = DG - DGTEST 442 DGXM = DGX - DGTEST 443 DGYM = DGY - DGTEST 444 ! 445 ! CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY 446 ! AND TO COMPUTE THE NEW STEP. 447 ! 448 call mcstep(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,BRACKT,STMIN,STMAX,INFOC) 449 ! 450 ! RESET THE FUNCTION AND GRADIENT VALUES FOR F. 451 ! 452 FX = FXM + STX * DGTEST 453 FY = FYM + STY * DGTEST 454 DGX = DGXM + DGTEST 455 DGY = DGYM + DGTEST 456 else 457 ! 458 ! CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY 459 ! AND TO COMPUTE THE NEW STEP. 460 ! 461 call mcstep(STX,FX,DGX,STY,FY,DGY,STP,F,DG,BRACKT,STMIN,STMAX,INFOC) 462 endif 463 ! 464 ! FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE 465 ! INTERVAL OF UNCERTAINTY. 466 ! 467 if (BRACKT) then 468 if (ABS(STY-STX) >= P66 * WIDTH1) STP = STX + P5 * (STY - STX) 469 WIDTH1 = WIDTH 470 WIDTH = ABS(STY-STX) 471 endif 472 473 else 474 475 INFOC = 1 476 ! 477 ! CHECK THE INPUT PARAMETERS FOR ERRORS. 478 ! 479 if ( STP <= zero .OR. FTOL < zero .OR. & 480 GTOL < zero .OR. XTOL < zero .OR. STPMIN < zero & 481 .OR. STPMAX < STPMIN ) return 482 ! 483 ! COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION 484 ! AND CHECK THAT S IS A DESCENT DIRECTION. 485 ! 486 DGINIT = DOT_PRODUCT( G , S ) 487 488 if (DGINIT > zero) then 489 return 490 endif 491 ! 492 ! INITIALIZE LOCAL VARIABLES. 493 ! 494 495 BRACKT = .FALSE. 496 STAGE1 = .TRUE. 497 NFEV = 0 498 FINIT = F 499 DGTEST = FTOL * DGINIT 500 WA(:) = X(:) 501 502 503 ! 504 ! THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP, 505 ! FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP. 506 ! THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP, 507 ! FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF 508 ! THE INTERVAL OF UNCERTAINTY. 509 ! THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP, 510 ! FUNCTION, AND DERIVATIVE AT THE CURRENT STEP. 511 ! 512 STX = zero 513 FX = FINIT 514 DGX = DGINIT 515 STY = zero 516 FY = FINIT 517 DGY = DGINIT 518 endif 519 520 ! 521 !SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND 522 !TO THE PRESENT INTERVAL OF UNCERTAINTY. 523 ! 524 if (BRACKT) then 525 STMIN = MIN(STX,STY) 526 STMAX = MAX(STX,STY) 527 else 528 STMIN = STX 529 STMAX = STP + XTRAPF*(STP - STX) 530 endif 531 ! 532 !FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN. 533 ! 534 STP = MAX(STPMIN,STP) 535 STP = MIN(STP,STPMAX) 536 ! 537 !IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET 538 !STP BE THE LOWEST POINT OBTAINED SO FAR. 539 ! 540 if ((BRACKT .AND. (STP <= STMIN .OR. STP >= STMAX)) & 541 .OR. NFEV >= MAXFEV-1 .OR. INFOC == 0 & 542 .OR. (BRACKT .AND. STMAX-STMIN <= XTOL*STMAX)) STP = STX 543 544 ! 545 !Evaluate the function and gradient at STP 546 !and compute the directional derivative. 547 !We return to main program to obtain F and G. 548 ! 549 X(:) = WA(:) + STP * S(:) 550 551 INFO = -1 552 553 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
SOURCE
574 subroutine mcstep(STX,FX,DX,STY,FY,DY,STP,FP,DG,BRACKT,STPMIN,STPMAX,INFO) 575 576 !Arguments ------------------------------------ 577 !scalars 578 integer,intent(inout) :: INFO 579 real(dp),intent(in) :: FP 580 real(dp),intent(inout) :: STX,FX,DX,STY,FY,DY,STP,DG,STPMIN,STPMAX 581 logical,intent(inout) :: BRACKT 582 583 !Local variables------------------------------- 584 !scalars 585 logical BOUND 586 real(dp) GAM,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA 587 !*************************************************************************** 588 589 INFO = 0 590 ! 591 ! CHECK THE INPUT PARAMETERS FOR ERRORS. 592 ! 593 IF ((BRACKT .AND. (STP <= MIN(STX,STY) .OR. & 594 STP >= MAX(STX,STY))) .OR. & 595 DX*(STP-STX) >= 0.0 .OR. STPMAX < STPMIN) RETURN 596 ! 597 ! Determine if the derivatives have opposite sign 598 ! 599 SGND = DG * ( DX / ABS(DX) ) 600 601 ! FIRST CASE. A HIGHER FUNCTION VALUE. 602 ! THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER 603 ! TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN, 604 ! ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN. 605 ! 606 IF (FP > FX) THEN 607 INFO = 1 608 BOUND = .TRUE. 609 THETA = 3*(FX - FP)/(STP - STX) + DX + DG 610 S = MAX(ABS(THETA),ABS(DX),ABS(DG)) 611 GAM = S * SQRT( (THETA/S)**2 - (DX/S)*(DG/S) ) 612 IF (STP < STX) GAM = -GAM 613 P = (GAM - DX) + THETA 614 Q = ((GAM - DX) + GAM) + DG 615 R = P / Q 616 STPC = STX + R*(STP - STX) 617 STPQ = STX + ( ( DX / ( ( FX - FP ) / ( STP - STX ) + DX ) ) / 2 ) * ( STP - STX ) 618 IF (ABS(STPC-STX) < ABS(STPQ-STX)) THEN 619 STPF = STPC 620 ELSE 621 STPF = STPC + (STPQ - STPC) / 2 622 END IF 623 BRACKT = .TRUE. 624 ! 625 ! SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF 626 ! OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC 627 ! STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP, 628 ! THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN. 629 ! 630 ELSE IF (SGND < 0.0) THEN 631 INFO = 2 632 BOUND = .FALSE. 633 THETA = 3*(FX - FP)/(STP - STX) + DX + DG 634 S = MAX(ABS(THETA),ABS(DX),ABS(DG)) 635 GAM = S * SQRT( (THETA/S)**2 - (DX/S)*(DG/S) ) 636 IF (STP > STX) GAM = -GAM 637 P = (GAM - DG) + THETA 638 Q = ((GAM - DG) + GAM) + DX 639 R = P/Q 640 STPC = STP + R*(STX - STP) 641 STPQ = STP + (DG/(DG-DX))*(STX - STP) 642 IF (ABS(STPC-STP) > ABS(STPQ-STP)) THEN 643 STPF = STPC 644 ELSE 645 STPF = STPQ 646 END IF 647 BRACKT = .TRUE. 648 ! 649 ! THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE 650 ! SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES. 651 ! THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY 652 ! IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC 653 ! IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE 654 ! EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO 655 ! COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP 656 ! CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN. 657 ! 658 ELSE IF (ABS(DG) < ABS(DX)) THEN 659 INFO = 3 660 BOUND = .TRUE. 661 THETA = 3*(FX - FP)/(STP - STX) + DX + DG 662 S = MAX(ABS(THETA),ABS(DX),ABS(DG)) 663 ! 664 ! THE CASE GAM = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND 665 ! TO INFINITY IN THE DIRECTION OF THE STEP. 666 ! 667 GAM = S * SQRT( MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DG/S)) ) 668 IF (STP > STX) GAM = -GAM 669 P = (GAM - DG) + THETA 670 Q = (GAM + (DX - DG)) + GAM 671 R = P/Q 672 IF (R < 0.0 .AND. GAM .NE. 0.0) THEN 673 STPC = STP + R*(STX - STP) 674 ELSE IF (STP > STX) THEN 675 STPC = STPMAX 676 ELSE 677 STPC = STPMIN 678 END IF 679 STPQ = STP + (DG/(DG-DX))*(STX - STP) 680 IF (BRACKT) THEN 681 IF (ABS(STP-STPC) < ABS(STP-STPQ)) THEN 682 STPF = STPC 683 ELSE 684 STPF = STPQ 685 END IF 686 ELSE 687 IF (ABS(STP-STPC) > ABS(STP-STPQ)) THEN 688 STPF = STPC 689 ELSE 690 STPF = STPQ 691 END IF 692 END IF 693 ! 694 ! FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE 695 ! SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES 696 ! NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP 697 ! IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN. 698 ! 699 ELSE 700 INFO = 4 701 BOUND = .FALSE. 702 IF (BRACKT) THEN 703 THETA = 3*(FP - FY)/(STY - STP) + DY + DG 704 S = MAX(ABS(THETA),ABS(DY),ABS(DG)) 705 GAM = S * SQRT( (THETA/S)**2 - (DY/S)*(DG/S) ) 706 IF (STP > STY) GAM = -GAM 707 P = (GAM - DG) + THETA 708 Q = ((GAM - DG) + GAM) + DY 709 R = P/Q 710 STPC = STP + R*(STY - STP) 711 STPF = STPC 712 ELSE IF (STP > STX) THEN 713 STPF = STPMAX 714 ELSE 715 STPF = STPMIN 716 END IF 717 END IF 718 719 ! 720 ! Update the interval of uncertainty. this update does not 721 ! depend on the new step or the case analysis above. 722 ! 723 IF (FP > FX) THEN 724 STY = STP 725 FY = FP 726 DY = DG 727 ELSE 728 IF (SGND < 0.0) THEN 729 STY = STX 730 FY = FX 731 DY = DX 732 END IF 733 STX = STP 734 FX = FP 735 DX = DG 736 END IF 737 738 ! 739 ! Compute the new step and safeguard it. 740 ! 741 STPF = MIN(STPMAX,STPF) 742 STPF = MAX(STPMIN,STPF) 743 STP = STPF 744 IF (BRACKT .AND. BOUND) THEN 745 IF (STY > STX) THEN 746 STP = MIN( STX + 0.66 * (STY-STX) , STP) 747 ELSE 748 STP = MAX( STX + 0.66 * (STY-STX) , STP) 749 END IF 750 END IF 751 752 753 end subroutine mcstep