Line data Source code
1 : !!****m* ABINIT/m_lbfgs
2 : !! NAME
3 : !! m_lbfgs
4 : !!
5 : !! FUNCTION
6 : !! This module provides several routines for the application of a
7 : !! Limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) minimization algorithm.
8 : !! The working routines were based on the original implementation of J. Nocera available on netlib.org
9 : !! They have been reshaped and translated into modern fortran here.
10 : !!
11 : !! COPYRIGHT
12 : !! Copyright (C) 2012-2026 ABINIT group (FB)
13 : !! This file is distributed under the terms of the
14 : !! GNU General Public License, see ~abinit/COPYING
15 : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : !!
17 : !! SOURCE
18 :
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
61 :
62 : !!***
63 :
64 : contains
65 :
66 : !----------------------------------------------------------------------
67 :
68 : !!****f* m_lbfgs/lbfgs_init
69 : !! NAME
70 : !! lbfgs_init
71 : !!
72 : !! FUNCTION
73 : !! Initialize the internal object lbfgs_internal for LBFGS minimization
74 : !!
75 : !! INPUTS
76 : !!
77 : !! OUTPUT
78 : !!
79 : !! SIDE EFFECTS
80 : !!
81 : !! SOURCE
82 :
83 1 : 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 1 : lbfgs_plan%lbfgs_status = 0
92 1 : lbfgs_plan%iter = 0
93 1 : lbfgs_plan%ndim = ndim
94 1 : lbfgs_plan%history_record = history_record
95 3 : ABI_MALLOC(lbfgs_plan%diag,(ndim))
96 :
97 1 : nwork = ndim * ( 2 * history_record + 1 ) + 2 * history_record
98 3 : ABI_MALLOC(lbfgs_plan%work,(nwork))
99 :
100 1 : lbfgs_plan%gtol = 0.9
101 1 : lbfgs_plan%line_stpmin = 1.0e-20
102 1 : lbfgs_plan%line_stpmax = 1.0e+20
103 1 : lbfgs_plan%line_stp = 1.0
104 :
105 28 : lbfgs_plan%diag(:) = diag_guess(:)
106 :
107 1 : end subroutine lbfgs_init
108 : !!***
109 :
110 :
111 : !----------------------------------------------------------------------
112 :
113 : !!****f* m_lbfgs/lbfgs_destroy
114 : !! NAME
115 : !! lbfgs_destroy
116 : !!
117 : !! FUNCTION
118 : !! Free the memory of the internal object lbfgs_internal for LBFGS minimization
119 : !!
120 : !! INPUTS
121 : !!
122 : !! OUTPUT
123 : !!
124 : !! SOURCE
125 :
126 1 : subroutine lbfgs_destroy()
127 :
128 1 : ABI_SFREE(lbfgs_plan%work)
129 1 : ABI_SFREE(lbfgs_plan%diag)
130 :
131 1 : end subroutine lbfgs_destroy
132 : !!***
133 :
134 : !----------------------------------------------------------------------
135 :
136 : !!****f* m_lbfgs/lbfgs_execute
137 : !! NAME
138 : !! lbfgs_execute
139 : !!
140 : !! FUNCTION
141 : !! Perform one-step of LBFGS minimization
142 : !! all the internal information are stored in the lbfgs_internal object
143 : !!
144 : !! INPUTS
145 : !! x: input and output position vector (atomic reduced coordinates + cell parameters)
146 : !! f: total energy
147 : !! gradf: gradient of the total energy (=negative forces)
148 : !!
149 : !! OUTPUT
150 : !!
151 : !! SOURCE
152 :
153 11 : 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 11 : lbfgs_plan%line_bracket, lbfgs_plan%line_stage1, lbfgs_plan%line_infoc)
168 :
169 :
170 : !lbfgs_execute = lbfgs_plan%lbfgs_status
171 11 : lbfgs_execute = lbfgs_plan%line_info
172 :
173 11 : end function lbfgs_execute
174 : !!***
175 :
176 :
177 : !----------------------------------------------------------------------
178 :
179 : !!****f* m_lbfgs/lbfgs
180 : !! NAME
181 : !! lbfgs
182 : !!
183 : !! FUNCTION
184 : !! Perform the LBFGS step
185 : !! Fortran90 rewritting of the original subroutine by J. Nocera
186 : !!
187 : !! INPUTS
188 : !!
189 : !! OUTPUT
190 : !!
191 : !! SIDE EFFECTS
192 : !!
193 : !! SOURCE
194 :
195 11 : 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 11 : FTOL = 1.0D-4
236 11 : MAXFEV = 20
237 :
238 11 : ISPT = N + 2 * M
239 11 : IYPT = ISPT + N * M
240 11 : POINT = MAX( 0 , MOD(ITER-1,M) )
241 11 : NPT = POINT * N
242 11 : ITER = ITER + 1
243 11 : 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 11 : if( IFLAG /= 1 ) then
250 28 : 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 10 : LINE_BRACKT,LINE_STAGE1,LINE_INFOC)
260 : !
261 : ! Compute the new step and gradient change
262 : !
263 10 : NPT = POINT * N
264 280 : W(ISPT+NPT+1:ISPT+NPT+N) = STP * W(ISPT+NPT+1:ISPT+NPT+N)
265 550 : W(IYPT+NPT+1:IYPT+NPT+N) = G(1:N) - W(1:N)
266 :
267 280 : YS = DOT_PRODUCT( W(IYPT+NPT+1:IYPT+NPT+N) , W(ISPT+NPT+1:ISPT+NPT+N) )
268 280 : YY = DOT_PRODUCT( W(IYPT+NPT+1:IYPT+NPT+N) , W(IYPT+NPT+1:IYPT+NPT+N) )
269 280 : 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 10 : POINT = MODULO(ITER - 1,M)
278 10 : CP = POINT
279 10 : if (POINT == 0) CP = M
280 :
281 10 : W(N+CP) = one / YS
282 280 : W(1:N) = -G(1:N)
283 :
284 : CP = POINT
285 50 : do I= 1,BOUND
286 40 : CP = CP - 1
287 40 : if (CP == -1) CP = M - 1
288 1120 : SQ = DOT_PRODUCT(W(ISPT+CP*N+1:ISPT+CP*N+N),W(1:N))
289 40 : INMC = N + M + CP + 1
290 40 : IYCN = IYPT + CP * N
291 40 : W(INMC)= W(N+CP+1) * SQ
292 1130 : W(1:N) = W(1:N) - W(INMC) * W(IYCN+1:IYCN+N)
293 : enddo
294 :
295 280 : W(1:N) = DIAG(1:N) * W(1:N)
296 :
297 50 : do I=1,BOUND
298 1120 : YR = DOT_PRODUCT(W(IYPT+CP*N+1:IYPT+CP*N+N),W(1:N))
299 40 : BETA = W(N+CP+1) * YR
300 40 : INMC = N + M + CP + 1
301 40 : BETA = W(INMC) - BETA
302 40 : ISCN = ISPT + CP * N
303 1120 : W(1:N) = W(1:N) + BETA * W(ISCN+1:ISCN+N)
304 40 : CP = CP + 1
305 50 : if (CP == M) CP = 0
306 : enddo
307 :
308 : !
309 : ! STORE THE NEW SEARCH DIRECTION
310 550 : 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 11 : NFEV = 0
319 11 : STP = one
320 308 : W(1:N) = G(1:N)
321 :
322 11 : 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 11 : LINE_BRACKT,LINE_STAGE1,LINE_INFOC)
330 :
331 11 : if (INFO == -1) then
332 11 : IFLAG = 1
333 11 : return
334 : else
335 0 : IFLAG = -1
336 0 : return
337 : endif
338 :
339 : end subroutine lbfgs
340 : !!***
341 :
342 : !----------------------------------------------------------------------
343 :
344 : !!****f* m_lbfgs/mcsrch
345 : !! NAME
346 : !! mcsrch
347 : !!
348 : !! FUNCTION
349 : !! Perform the line minimization step
350 : !! Fortran90 rewritting of the original subroutine by J. Nocera
351 : !!
352 : !! INPUTS
353 : !!
354 : !! OUTPUT
355 : !!
356 : !! SIDE EFFECTS
357 : !!
358 : !! SOURCE
359 :
360 21 : 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 21 : DGTEST = FTOL * DGINIT
392 21 : WIDTH = STPMAX - STPMIN
393 21 : WIDTH1 = WIDTH / P5
394 :
395 : ! Is it a first entry (info == 0)
396 : ! or a second entry (info == -1)?
397 21 : if( INFO == -1 ) then
398 :
399 : ! Reset INFO
400 10 : INFO = 0
401 :
402 10 : NFEV = NFEV + 1
403 280 : DG = SUM( G(:) * S(:) )
404 10 : FTEST1 = FINIT + STP * DGTEST
405 : !
406 : ! TEST FOR CONVERGENCE.
407 : !
408 : if ((BRACKT .AND. (STP <= STMIN .OR. STP >= STMAX)) &
409 10 : .OR. INFOC == 0) INFO = 6
410 : if (STP == STPMAX .AND. &
411 10 : F <= FTEST1 .AND. DG <= DGTEST) INFO = 5
412 10 : if (STP == STPMIN .AND. &
413 0 : (F > FTEST1 .OR. DG >= DGTEST)) INFO = 4
414 10 : if (NFEV >= MAXFEV) INFO = 3
415 10 : if (BRACKT .AND. STMAX-STMIN <= XTOL*STMAX) INFO = 2
416 10 : if (F <= FTEST1 .AND. ABS(DG) <= GTOL*(-DGINIT)) INFO = 1
417 : !
418 : ! CHECK FOR TERMINATION.
419 : !
420 10 : 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 1 : if (STAGE1 .AND. F <= FTEST1 .AND. &
426 0 : 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 1 : if (STAGE1 .AND. F <= FX .AND. F > FTEST1) then
435 : !
436 : ! DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
437 : !
438 0 : FM = F - STP * DGTEST
439 0 : FXM = FX - STX * DGTEST
440 0 : FYM = FY - STY * DGTEST
441 0 : DGM = DG - DGTEST
442 0 : DGXM = DGX - DGTEST
443 0 : DGYM = DGY - DGTEST
444 : !
445 : ! CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
446 : ! AND TO COMPUTE THE NEW STEP.
447 : !
448 0 : 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 0 : FX = FXM + STX * DGTEST
453 0 : FY = FYM + STY * DGTEST
454 0 : DGX = DGXM + DGTEST
455 0 : DGY = DGYM + DGTEST
456 : else
457 : !
458 : ! CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
459 : ! AND TO COMPUTE THE NEW STEP.
460 : !
461 1 : 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 1 : if (BRACKT) then
468 1 : if (ABS(STY-STX) >= P66 * WIDTH1) STP = STX + P5 * (STY - STX)
469 12 : WIDTH1 = WIDTH
470 12 : WIDTH = ABS(STY-STX)
471 : endif
472 :
473 : else
474 :
475 11 : 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 11 : .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 308 : DGINIT = DOT_PRODUCT( G , S )
487 :
488 11 : if (DGINIT > zero) then
489 : return
490 : endif
491 : !
492 : ! INITIALIZE LOCAL VARIABLES.
493 : !
494 :
495 11 : BRACKT = .FALSE.
496 11 : STAGE1 = .TRUE.
497 11 : NFEV = 0
498 11 : FINIT = F
499 11 : DGTEST = FTOL * DGINIT
500 308 : 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 11 : STX = zero
513 11 : FX = FINIT
514 11 : DGX = DGINIT
515 11 : STY = zero
516 11 : FY = FINIT
517 11 : DGY = DGINIT
518 : endif
519 :
520 : !
521 : !SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
522 : !TO THE PRESENT INTERVAL OF UNCERTAINTY.
523 : !
524 12 : if (BRACKT) then
525 1 : STMIN = MIN(STX,STY)
526 1 : STMAX = MAX(STX,STY)
527 : else
528 11 : STMIN = STX
529 11 : STMAX = STP + XTRAPF*(STP - STX)
530 : endif
531 : !
532 : !FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
533 : !
534 12 : STP = MAX(STPMIN,STP)
535 12 : 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 12 : .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 336 : X(:) = WA(:) + STP * S(:)
550 :
551 12 : INFO = -1
552 :
553 : end subroutine mcsrch
554 : !!***
555 :
556 : !----------------------------------------------------------------------
557 :
558 : !!****f* m_lbfgs/mcstep
559 : !! NAME
560 : !! mcstep
561 : !!
562 : !! FUNCTION
563 : !! Perform the step choice in line minimization
564 : !! Fortran90 rewritting of the original subroutine by J. Nocera
565 : !!
566 : !! INPUTS
567 : !!
568 : !! OUTPUT
569 : !!
570 : !! SIDE EFFECTS
571 : !!
572 : !! SOURCE
573 :
574 1 : 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 1 : 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 1 : DX*(STP-STX) >= 0.0 .OR. STPMAX < STPMIN) RETURN
596 : !
597 : ! Determine if the derivatives have opposite sign
598 : !
599 1 : 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 1 : IF (FP > FX) THEN
607 1 : INFO = 1
608 1 : BOUND = .TRUE.
609 1 : THETA = 3*(FX - FP)/(STP - STX) + DX + DG
610 1 : S = MAX(ABS(THETA),ABS(DX),ABS(DG))
611 1 : GAM = S * SQRT( (THETA/S)**2 - (DX/S)*(DG/S) )
612 1 : IF (STP < STX) GAM = -GAM
613 1 : P = (GAM - DX) + THETA
614 1 : Q = ((GAM - DX) + GAM) + DG
615 1 : R = P / Q
616 1 : STPC = STX + R*(STP - STX)
617 1 : STPQ = STX + ( ( DX / ( ( FX - FP ) / ( STP - STX ) + DX ) ) / 2 ) * ( STP - STX )
618 1 : IF (ABS(STPC-STX) < ABS(STPQ-STX)) THEN
619 : STPF = STPC
620 : ELSE
621 0 : STPF = STPC + (STPQ - STPC) / 2
622 : END IF
623 1 : 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 0 : ELSE IF (SGND < 0.0) THEN
631 0 : INFO = 2
632 0 : BOUND = .FALSE.
633 0 : THETA = 3*(FX - FP)/(STP - STX) + DX + DG
634 0 : S = MAX(ABS(THETA),ABS(DX),ABS(DG))
635 0 : GAM = S * SQRT( (THETA/S)**2 - (DX/S)*(DG/S) )
636 0 : IF (STP > STX) GAM = -GAM
637 0 : P = (GAM - DG) + THETA
638 0 : Q = ((GAM - DG) + GAM) + DX
639 0 : R = P/Q
640 0 : STPC = STP + R*(STX - STP)
641 0 : STPQ = STP + (DG/(DG-DX))*(STX - STP)
642 0 : IF (ABS(STPC-STP) > ABS(STPQ-STP)) THEN
643 : STPF = STPC
644 : ELSE
645 0 : STPF = STPQ
646 : END IF
647 0 : 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 0 : ELSE IF (ABS(DG) < ABS(DX)) THEN
659 0 : INFO = 3
660 0 : BOUND = .TRUE.
661 0 : THETA = 3*(FX - FP)/(STP - STX) + DX + DG
662 0 : 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 0 : GAM = S * SQRT( MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DG/S)) )
668 0 : IF (STP > STX) GAM = -GAM
669 0 : P = (GAM - DG) + THETA
670 0 : Q = (GAM + (DX - DG)) + GAM
671 0 : R = P/Q
672 0 : IF (R < 0.0 .AND. GAM .NE. 0.0) THEN
673 0 : STPC = STP + R*(STX - STP)
674 0 : ELSE IF (STP > STX) THEN
675 : STPC = STPMAX
676 : ELSE
677 0 : STPC = STPMIN
678 : END IF
679 0 : STPQ = STP + (DG/(DG-DX))*(STX - STP)
680 0 : IF (BRACKT) THEN
681 0 : IF (ABS(STP-STPC) < ABS(STP-STPQ)) THEN
682 : STPF = STPC
683 : ELSE
684 0 : STPF = STPQ
685 : END IF
686 : ELSE
687 0 : IF (ABS(STP-STPC) > ABS(STP-STPQ)) THEN
688 : STPF = STPC
689 : ELSE
690 0 : 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 0 : INFO = 4
701 0 : BOUND = .FALSE.
702 0 : IF (BRACKT) THEN
703 0 : THETA = 3*(FP - FY)/(STY - STP) + DY + DG
704 0 : S = MAX(ABS(THETA),ABS(DY),ABS(DG))
705 0 : GAM = S * SQRT( (THETA/S)**2 - (DY/S)*(DG/S) )
706 0 : IF (STP > STY) GAM = -GAM
707 0 : P = (GAM - DG) + THETA
708 0 : Q = ((GAM - DG) + GAM) + DY
709 0 : R = P/Q
710 0 : STPC = STP + R*(STY - STP)
711 0 : STPF = STPC
712 0 : ELSE IF (STP > STX) THEN
713 : STPF = STPMAX
714 : ELSE
715 0 : 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 1 : IF (FP > FX) THEN
724 1 : STY = STP
725 1 : FY = FP
726 1 : DY = DG
727 : ELSE
728 0 : IF (SGND < 0.0) THEN
729 0 : STY = STX
730 0 : FY = FX
731 0 : DY = DX
732 : END IF
733 0 : STX = STP
734 0 : FX = FP
735 0 : DX = DG
736 : END IF
737 :
738 : !
739 : ! Compute the new step and safeguard it.
740 : !
741 1 : STPF = MIN(STPMAX,STPF)
742 1 : STPF = MAX(STPMIN,STPF)
743 1 : STP = STPF
744 1 : IF (BRACKT .AND. BOUND) THEN
745 1 : IF (STY > STX) THEN
746 1 : STP = MIN( STX + 0.66 * (STY-STX) , STP)
747 : ELSE
748 0 : STP = MAX( STX + 0.66 * (STY-STX) , STP)
749 : END IF
750 : END IF
751 :
752 :
753 : end subroutine mcstep
754 : !!***
755 :
756 0 : end module m_lbfgs
757 : !!***
|
