Branch data 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 : : !!***
|
