Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief LBFGS-B routine (version 3.0, April 25, 2011)
10 : !> \note
11 : !> L-BFGS-B (version 3.0, April 25, 2011) converted to Fortran 90 module
12 : !> \par History
13 : !> 02.2005 Update to the new version 2.4 and deleting the blas part of
14 : !> the code (Teodoro Laino)
15 : !> 11.2012 New version 3.0 converted to Fortran 90 (Matthias Krack)
16 : !> 12.2020 Implementation of Space Group Symmetry (Pierre-André Cazade)
17 : !> \author Fawzi Mohamed (first version)
18 : ! **************************************************************************************************
19 : MODULE cp_lbfgs
20 : USE bibliography, ONLY: Byrd1995,&
21 : cite_reference
22 : USE cp_files, ONLY: open_file
23 : USE kinds, ONLY: dp
24 : USE machine, ONLY: default_output_unit,&
25 : m_walltime
26 : USE space_groups, ONLY: spgr_apply_rotations_coord,&
27 : spgr_apply_rotations_force
28 : USE space_groups_types, ONLY: spgr_type
29 : #include "../base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 :
33 : PRIVATE
34 :
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_lbfgs'
36 :
37 : PUBLIC :: setulb
38 :
39 : CONTAINS
40 :
41 : !=========== L-BFGS-B (version 3.0, April 25, 2011) ================
42 : !
43 : ! This is a modified version of L-BFGS-B.
44 : !
45 : ! Major changes are described in the accompanying paper:
46 : !
47 : ! Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778:
48 : ! L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constraine
49 : ! Optimization" (2011). To appear in ACM Transactions on
50 : ! Mathematical Software,
51 : !
52 : ! The paper describes an improvement and a correction to Algorithm 7
53 : ! It is shown that the performance of the algorithm can be improved
54 : ! significantly by making a relatively simple modication to the subs
55 : ! minimization phase. The correction concerns an error caused by the
56 : ! of routine dpmeps to estimate machine precision.
57 : !
58 : ! The total work space **wa** required by the new version is
59 : !
60 : ! 2*m*n + 11m*m + 5*n + 8*m
61 : !
62 : ! the old version required
63 : !
64 : ! 2*m*n + 12m*m + 4*n + 12*m
65 : !
66 : !
67 : ! J. Nocedal Department of Electrical Engineering and
68 : ! Computer Science.
69 : ! Northwestern University. Evanston, IL. USA
70 : !
71 : !
72 : ! J.L Morales Departamento de Matematicas,
73 : ! Instituto Tecnologico Autonomo de Mexico
74 : ! Mexico D.F. Mexico.
75 : !
76 : ! March 2011
77 : !
78 : !=======================================================================
79 : ! **************************************************************************************************
80 : !> \brief This subroutine partitions the working arrays wa and iwa, and
81 : !> then uses the limited memory BFGS method to solve the bound
82 : !> constrained optimization problem by calling mainlb.
83 : !> (The direct method will be used in the subspace minimization.)
84 : !> \param n n is the dimension of the problem.
85 : !> \param m m is the maximum number of variable metric corrections
86 : !> used to define the limited memory matrix.
87 : !> \param x On entry x is an approximation to the solution.
88 : !> On exit x is the current approximation.
89 : !> \param lower_bound the lower bound on x.
90 : !> \param upper_bound the upper bound on x.
91 : !> \param nbd nbd represents the type of bounds imposed on the
92 : !> variables, and must be specified as follows:
93 : !> nbd(i)=0 if x(i) is unbounded,
94 : !> 1 if x(i) has only a lower bound,
95 : !> 2 if x(i) has both lower and upper bounds, and
96 : !> 3 if x(i) has only an upper bound.
97 : !> \param f On first entry f is unspecified.
98 : !> On final exit f is the value of the function at x.
99 : !> \param g On first entry g is unspecified.
100 : !> On final exit g is the value of the gradient at x.
101 : !> \param factr factr >= 0 is specified by the user. The iteration
102 : !> will stop when
103 : !>
104 : !> (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
105 : !>
106 : !> where epsmch is the machine precision, which is automatically
107 : !> generated by the code. Typical values for factr: 1.d+12 for
108 : !> low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely
109 : !> high accuracy.
110 : !> \param pgtol pgtol >= 0 is specified by the user. The iteration
111 : !> will stop when
112 : !>
113 : !> max{|proj g_i | i = 1, ..., n} <= pgtol
114 : !>
115 : !> where pg_i is the ith component of the projected gradient.
116 : !> \param wa working array
117 : !> \param iwa integer working array
118 : !> \param task is a working string of characters of length 60 indicating
119 : !> the current job when entering and quitting this subroutine.
120 : !> \param iprint iprint is a variable that must be set by the user.
121 : !> It controls the frequency and type of output generated:
122 : !> iprint<0 no output is generated;
123 : !> iprint=0 print only one line at the last iteration;
124 : !> 0<iprint<99 print also f and |proj g| every iprint iterations;
125 : !> iprint=99 print details of every iteration except n-vectors;
126 : !> iprint=100 print also the changes of active set and final x;
127 : !> iprint>100 print details of every iteration including x and g;
128 : !> When iprint > 0, the file iterate.dat will be created to
129 : !> summarize the iteration.
130 : !> \param csave is a working string of characters
131 : !> \param lsave lsave is a working array
132 : !> On exit with 'task' = NEW_X, the following information is available:
133 : !> If lsave(1) = .true. then the initial X has been replaced by
134 : !> its projection in the feasible set
135 : !> If lsave(2) = .true. then the problem is constrained;
136 : !> If lsave(3) = .true. then each variable has upper and lower bounds;
137 : !> \param isave isave is a working array
138 : !> On exit with 'task' = NEW_X, the following information is available:
139 : !> isave(22) = the total number of intervals explored in the
140 : !> search of Cauchy points;
141 : !> isave(26) = the total number of skipped BFGS updates before the current iteration;
142 : !> isave(30) = the number of current iteration;
143 : !> isave(31) = the total number of BFGS updates prior the current iteration;
144 : !> isave(33) = the number of intervals explored in the search of
145 : !> Cauchy point in the current iteration;
146 : !> isave(34) = the total number of function and gradient evaluations;
147 : !> isave(36) = the number of function value or gradient
148 : !> evaluations in the current iteration;
149 : !> if isave(37) = 0 then the subspace argmin is within the box;
150 : !> if isave(37) = 1 then the subspace argmin is beyond the box;
151 : !> isave(38) = the number of free variables in the current iteration;
152 : !> isave(39) = the number of active constraints in the current iteration;
153 : !> n + 1 - isave(40) = the number of variables leaving the set of
154 : !> active constraints in the current iteration;
155 : !> isave(41) = the number of variables entering the set of active
156 : !> constraints in the current iteration.
157 : !> \param dsave dsave is a working array of dimension 29.
158 : !> On exit with 'task' = NEW_X, the following information is available:
159 : !> dsave(1) = current 'theta' in the BFGS matrix;
160 : !> dsave(2) = f(x) in the previous iteration;
161 : !> dsave(3) = factr*epsmch;
162 : !> dsave(4) = 2-norm of the line search direction vector;
163 : !> dsave(5) = the machine precision epsmch generated by the code;
164 : !> dsave(7) = the accumulated time spent on searching for Cauchy points;
165 : !> dsave(8) = the accumulated time spent on subspace minimization;
166 : !> dsave(9) = the accumulated time spent on line search;
167 : !> dsave(11) = the slope of the line search function at the current point of line search;
168 : !> dsave(12) = the maximum relative step length imposed in line search;
169 : !> dsave(13) = the infinity norm of the projected gradient;
170 : !> dsave(14) = the relative step length in the line search;
171 : !> dsave(15) = the slope of the line search function at the starting point of the line search;
172 : !> dsave(16) = the square of the 2-norm of the line search direction vector.
173 : !> \param trust_radius ...
174 : !> \param spgr ...
175 : !> \param iwunit User-specified write unit, if not set then WRITE statements
176 : !> write to default_output_unit by default
177 : !> \par History
178 : !> 12.2020 Implementation of Space Group Symmetry [pcazade]
179 : !> \author NEOS, November 1994. (Latest revision June 1996.)
180 : !> Optimization Technology Center.
181 : !> Argonne National Laboratory and Northwestern University.
182 : !> Written by
183 : !> Ciyou Zhu
184 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
185 : ! **************************************************************************************************
186 3181 : SUBROUTINE setulb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, wa, iwa, &
187 : task, iprint, csave, lsave, isave, dsave, trust_radius, spgr, iwunit)
188 :
189 : INTEGER, INTENT(in) :: n, m
190 : REAL(KIND=dp), INTENT(inout) :: x(n)
191 : REAL(KIND=dp) :: lower_bound(n), upper_bound(n)
192 : INTEGER :: nbd(n)
193 : REAL(KIND=dp) :: f, g(n)
194 : REAL(KIND=dp), INTENT(in) :: factr, pgtol
195 : REAL(KIND=dp) :: wa(2*m*n + 5*n + 11*m*m + 8*m)
196 : INTEGER :: iwa(3*n)
197 : CHARACTER(LEN=60) :: task
198 : INTEGER :: iprint
199 : CHARACTER(LEN=60) :: csave
200 : LOGICAL :: lsave(4)
201 : INTEGER :: isave(44)
202 : REAL(KIND=dp) :: dsave(29)
203 : REAL(KIND=dp), INTENT(in) :: trust_radius
204 : TYPE(spgr_type), OPTIONAL, POINTER :: spgr
205 : INTEGER, OPTIONAL :: iwunit
206 :
207 : INTEGER :: i, ld, lr, lsnd, lss, lsy, lt, lwa, lwn, &
208 : lws, lwt, lwy, lxp, lz, wunit
209 :
210 : ! References:
211 : !
212 : ! [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
213 : ! memory algorithm for bound constrained optimization'',
214 : ! SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
215 : !
216 : ! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
217 : ! limited memory FORTRAN code for solving bound constrained
218 : ! optimization problems'', Tech. Report, NAM-11, EECS Department,
219 : ! Northwestern University, 1994.
220 : !
221 : ! (Postscript files of these papers are available via anonymous
222 : ! ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
223 : !
224 : ! * * *
225 :
226 3181 : wunit = default_output_unit
227 3181 : IF (PRESENT(iwunit)) THEN
228 3136 : IF (iwunit > 0) wunit = iwunit
229 : END IF
230 :
231 3181 : IF (task == 'START') THEN
232 45 : CALL cite_reference(Byrd1995)
233 45 : isave(1) = m*n
234 45 : isave(2) = m**2
235 45 : isave(3) = 4*m**2
236 : ! ws m*n
237 45 : isave(4) = 1
238 : ! wy m*n
239 45 : isave(5) = isave(4) + isave(1)
240 : ! wsy m**2
241 45 : isave(6) = isave(5) + isave(1)
242 : ! wss m**2
243 45 : isave(7) = isave(6) + isave(2)
244 : ! wt m**2
245 45 : isave(8) = isave(7) + isave(2)
246 : ! wn 4*m**2
247 45 : isave(9) = isave(8) + isave(2)
248 : ! wsnd 4*m**2
249 45 : isave(10) = isave(9) + isave(3)
250 : ! wz n
251 45 : isave(11) = isave(10) + isave(3)
252 : ! wr n
253 45 : isave(12) = isave(11) + n
254 : ! wd n
255 45 : isave(13) = isave(12) + n
256 : ! wt n
257 45 : isave(14) = isave(13) + n
258 : ! wxp n
259 45 : isave(15) = isave(14) + n
260 : ! wa 8*m
261 45 : isave(16) = isave(15) + n
262 : END IF
263 3181 : lws = isave(4)
264 3181 : lwy = isave(5)
265 3181 : lsy = isave(6)
266 3181 : lss = isave(7)
267 3181 : lwt = isave(8)
268 3181 : lwn = isave(9)
269 3181 : lsnd = isave(10)
270 3181 : lz = isave(11)
271 3181 : lr = isave(12)
272 3181 : ld = isave(13)
273 3181 : lt = isave(14)
274 3181 : lxp = isave(15)
275 3181 : lwa = isave(16)
276 :
277 : !in case we use a trust radius we set the boundaries to be one times the trust radius away from the current positions
278 : !the original implementation only allowed for boundaries that remain constant during the optimization.
279 : !This way of including a trust radius seems to work,
280 : !but the change of the boundaries during optimization might introduce some not yet discovered problems.
281 3181 : IF (trust_radius >= 0) THEN
282 81371 : DO i = 1, n
283 81318 : lower_bound(i) = x(i) - trust_radius
284 81318 : upper_bound(i) = x(i) + trust_radius
285 81371 : nbd(i) = 2
286 : END DO
287 : END IF
288 :
289 : ! passes spgr and wunit to mainlb
290 : CALL mainlb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, &
291 : wa(lws), wa(lwy), wa(lsy), wa(lss), wa(lwt), &
292 : wa(lwn), wa(lsnd), wa(lz), wa(lr), wa(ld), wa(lt), wa(lxp), &
293 : wa(lwa), &
294 : iwa(1), iwa(n + 1), iwa(2*n + 1), task, iprint, &
295 3181 : csave, lsave, isave(22), dsave, spgr, wunit)
296 :
297 3181 : RETURN
298 :
299 : END SUBROUTINE setulb
300 :
301 : ! **************************************************************************************************
302 : !> \brief This subroutine solves bound constrained optimization problems by
303 : !> using the compact formula of the limited memory BFGS updates.
304 : !> \param n n is the number of variables
305 : !> \param m m is the maximum number of variable metric
306 : !> corrections allowed in the limited memory matrix.
307 : !> \param x On entry x is an approximation to the solution.
308 : !> On exit x is the current approximation.
309 : !> \param lower_bound lower_bound is the lower bound of x.
310 : !> \param upper_bound upper_bound is the upper bound of x.
311 : !> \param nbd nbd represents the type of bounds imposed on the
312 : !> variables, and must be specified as follows:
313 : !> nbd(i)=0 if x(i) is unbounded,
314 : !> 1 if x(i) has only a lower bound,
315 : !> 2 if x(i) has both lower and upper bounds,
316 : !> 3 if x(i) has only an upper bound.
317 : !> \param f On first entry f is unspecified.
318 : !> On final exit f is the value of the function at x.
319 : !> \param g On first entry g is unspecified.
320 : !> On final exit g is the value of the gradient at x.
321 : !> \param factr factr >= 0 is specified by the user. The iteration
322 : !> will stop when
323 : !>
324 : !> (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
325 : !>
326 : !> where epsmch is the machine precision, which is automatically
327 : !> generated by the code.
328 : !> \param pgtol pgtol >= 0 is specified by the user. The iteration
329 : !> will stop when
330 : !>
331 : !> max{|proj g_i | i = 1, ..., n} <= pgtol
332 : !>
333 : !> where pg_i is the ith component of the projected gradient.
334 : !> \param ws ws, wy, sy, and wt are working arrays used to store the following
335 : !> information defining the limited memory BFGS matrix:
336 : !> ws stores S, the matrix of s-vectors;
337 : !> \param wy stores Y, the matrix of y-vectors;
338 : !> \param sy stores S'Y;
339 : !> \param ss stores S'S;
340 : !> \param wt stores the Cholesky factorization of (theta*S'S+LD^(-1)L');
341 : !> see eq. (2.26) in [3].
342 : !> \param wn wn is a working array of dimension 2m x 2m
343 : !> used to store the LEL^T factorization of the indefinite matrix
344 : !> K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
345 : !> [L_a -R_z theta*S'AA'S ]
346 : !>
347 : !> where E = [-I 0]
348 : !> [ 0 I]
349 : !> \param snd is a working array of dimension 2m x 2m
350 : !> used to store the lower triangular part of
351 : !> N = [Y' ZZ'Y L_a'+R_z']
352 : !> [L_a +R_z S'AA'S ]
353 : !> \param z z(n),r(n),d(n),t(n), xp(n),wa(8*m) are working arrays
354 : !> z is used at different times to store the Cauchy point and
355 : !> the Newton point.
356 : !> \param r working array
357 : !> \param d working array
358 : !> \param t workign array
359 : !> \param xp xp is a workng array used to safeguard the projected Newton direction
360 : !> \param wa working array
361 : !> \param index In subroutine freev, index is used to store the free and fixed
362 : !> variables at the Generalized Cauchy Point (GCP).
363 : !> \param iwhere iwhere is an integer working array of dimension n used to record
364 : !> the status of the vector x for GCP computation.
365 : !> iwhere(i)=0 or -3 if x(i) is free and has bounds,
366 : !> 1 if x(i) is fixed at l(i), and l(i) .ne. u(i)
367 : !> 2 if x(i) is fixed at u(i), and u(i) .ne. l(i)
368 : !> 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i)
369 : !> -1 if x(i) is always free, i.e., no bounds on it.
370 : !> \param indx2 indx2 is a working array. Within subroutine cauchy, indx2 corresponds to the array iorder.
371 : !> In subroutine freev, a list of variables entering and leaving
372 : !> the free set is stored in indx2, and it is passed on to
373 : !> subroutine formk with this information
374 : !> \param task task is a working string of characters indicating
375 : !> the current job when entering and leaving this subroutine.
376 : !> \param iprint is an variable that must be set by the user.
377 : !> It controls the frequency and type of output generated:
378 : !> iprint<0 no output is generated;
379 : !> iprint=0 print only one line at the last iteration;
380 : !> 0<iprint<99 print also f and |proj g| every iprint iterations;
381 : !> iprint=99 print details of every iteration except n-vectors;
382 : !> iprint=100 print also the changes of active set and final x;
383 : !> iprint>100 print details of every iteration including x and g;
384 : !> When iprint > 0, the file iterate.dat will be created to summarize the iteration.
385 : !> \param csave csave is a working string of characters
386 : !> \param lsave lsave is a logical working array
387 : !> \param isave isave is an integer working array
388 : !> \param dsave is a double precision working array
389 : !> \param spgr ...
390 : !> \param iwunit User-specified write unit, if not set then WRITE statements
391 : !> write to default_output_unit by default
392 : !> \par History
393 : !> 12.2020 Implementation of Space Group Symmetry [pcazade]
394 : !> \author NEOS, November 1994. (Latest revision June 1996.)
395 : !> Optimization Technology Center.
396 : !> Argonne National Laboratory and Northwestern University.
397 : !> Written by
398 : !> Ciyou Zhu
399 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
400 : ! **************************************************************************************************
401 3181 : SUBROUTINE mainlb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, ws, wy, &
402 3181 : sy, ss, wt, wn, snd, z, r, d, t, xp, wa, &
403 3181 : index, iwhere, indx2, task, &
404 : iprint, csave, lsave, isave, dsave, spgr, iwunit)
405 : INTEGER, INTENT(in) :: n, m
406 : REAL(KIND=dp), INTENT(inout) :: x(n)
407 : REAL(KIND=dp), INTENT(in) :: lower_bound(n), upper_bound(n)
408 : INTEGER :: nbd(n)
409 : REAL(KIND=dp) :: f, g(n), factr, pgtol, ws(n, m), wy(n, m), sy(m, m), ss(m, m), wt(m, m), &
410 : wn(2*m, 2*m), snd(2*m, 2*m), z(n), r(n), d(n), t(n), xp(n), wa(8*m)
411 : INTEGER :: INDEX(n), iwhere(n), indx2(n)
412 : CHARACTER(LEN=60) :: task
413 : INTEGER :: iprint
414 : CHARACTER(LEN=60) :: csave
415 : LOGICAL :: lsave(4)
416 : INTEGER :: isave(23)
417 : REAL(KIND=dp) :: dsave(29)
418 : TYPE(spgr_type), OPTIONAL, POINTER :: spgr
419 : INTEGER, OPTIONAL :: iwunit
420 :
421 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
422 :
423 : CHARACTER(LEN=3) :: word
424 : INTEGER :: col, head, i, iback, ifun, ileave, info, &
425 : itail, iter, itfile, iupdat, iword, k, &
426 : nact, nenter, nfgv, nfree, nintol, &
427 : nseg, nskip, wunit
428 : LOGICAL :: boxed, constrained, first, &
429 : keep_space_group, updatd, wrk, &
430 : x_projected
431 : REAL(KIND=dp) :: cachyt, cpu1, cpu2, ddot, ddum, dnorm, dr, dtd, epsmch, fold, g_inf_norm, &
432 : gd, gdold, lnscht, rr, sbtime, step_max, stp, theta, time, time1, time2, tol, xstep
433 :
434 : ! References:
435 : !
436 : ! [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
437 : ! memory algorithm for bound constrained optimization'',
438 : ! SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
439 : !
440 : ! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
441 : ! Subroutines for Large Scale Bound Constrained Optimization''
442 : ! Tech. Report, NAM-11, EECS Department, Northwestern University,
443 : ! 1994.
444 : !
445 : ! [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of
446 : ! Quasi-Newton Matrices and their use in Limited Memory Methods'',
447 : ! Mathematical Programming 63 (1994), no. 4, pp. 129-156.
448 : !
449 : ! (Postscript files of these papers are available via anonymous
450 : ! ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
451 : !
452 : ! * * *
453 :
454 3181 : wunit = default_output_unit
455 3181 : IF (PRESENT(iwunit)) THEN
456 3181 : IF (iwunit > 0) wunit = iwunit
457 : END IF
458 :
459 3181 : keep_space_group = .FALSE.
460 3181 : IF (PRESENT(spgr)) THEN
461 3136 : IF (ASSOCIATED(spgr)) keep_space_group = spgr%keep_space_group
462 : END IF
463 :
464 3181 : IF (task == 'START') THEN
465 :
466 45 : epsmch = EPSILON(one)
467 :
468 45 : CALL timer(time1)
469 :
470 : ! Initialize counters and scalars when task='START'.
471 :
472 : ! for the limited memory BFGS matrices:
473 45 : col = 0
474 45 : head = 1
475 45 : theta = one
476 45 : iupdat = 0
477 45 : updatd = .FALSE.
478 45 : iback = 0
479 45 : itail = 0
480 45 : iword = 0
481 45 : nact = 0
482 45 : ileave = 0
483 45 : nenter = 0
484 45 : fold = zero
485 45 : dnorm = zero
486 45 : cpu1 = zero
487 45 : gd = zero
488 45 : step_max = zero
489 45 : g_inf_norm = zero
490 45 : stp = zero
491 45 : gdold = zero
492 45 : dtd = zero
493 :
494 : ! for operation counts:
495 45 : iter = 0
496 45 : nfgv = 0
497 45 : nseg = 0
498 45 : nintol = 0
499 45 : nskip = 0
500 45 : nfree = n
501 45 : ifun = 0
502 : ! for stopping tolerance:
503 45 : tol = factr*epsmch
504 :
505 : ! for measuring running time:
506 45 : cachyt = 0
507 45 : sbtime = 0
508 45 : lnscht = 0
509 :
510 : ! 'word' records the status of subspace solutions.
511 45 : word = '---'
512 :
513 : ! 'info' records the termination information.
514 45 : info = 0
515 :
516 45 : itfile = 8
517 45 : IF (iprint >= 1) THEN
518 : ! open a summary file 'iterate.dat'
519 45 : CALL open_file(file_name='iterate.dat', unit_number=itfile, file_action='WRITE', file_status='UNKNOWN')
520 : END IF
521 :
522 : ! Check the input arguments for errors.
523 :
524 45 : CALL errclb(n, m, factr, lower_bound, upper_bound, nbd, task, info, k)
525 45 : IF (task(1:5) == 'ERROR') THEN
526 : CALL prn3lb(n, x, f, task, iprint, info, itfile, &
527 : iter, nfgv, nintol, nskip, nact, g_inf_norm, &
528 : zero, nseg, word, iback, stp, xstep, k, &
529 0 : cachyt, sbtime, lnscht, wunit)
530 0 : RETURN
531 : END IF
532 :
533 45 : CALL prn1lb(n, m, lower_bound, upper_bound, x, iprint, itfile, epsmch, wunit)
534 :
535 : ! Initialize iwhere & project x onto the feasible set.
536 :
537 45 : CALL active(n, lower_bound, upper_bound, nbd, x, iwhere, iprint, x_projected, constrained, boxed, wunit)
538 : ! applies rotation matrices to coordinates
539 45 : IF (keep_space_group) THEN
540 0 : CALL spgr_apply_rotations_coord(spgr, x)
541 : END IF
542 :
543 : ! The end of the initialization.
544 45 : task = 'FG_START'
545 : ! return to the driver to calculate f and g; reenter at 111.
546 : CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
547 : iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
548 45 : cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
549 45 : RETURN
550 : ELSE
551 : ! applies rotation matrices to coordinates
552 3136 : IF (keep_space_group) THEN
553 2 : CALL spgr_apply_rotations_coord(spgr, x)
554 2 : CALL spgr_apply_rotations_force(spgr, g)
555 : END IF
556 :
557 : ! restore local variables.
558 :
559 3136 : x_projected = lsave(1)
560 3136 : constrained = lsave(2)
561 3136 : boxed = lsave(3)
562 3136 : updatd = lsave(4)
563 :
564 3136 : nintol = isave(1)
565 3136 : itfile = isave(3)
566 3136 : iback = isave(4)
567 3136 : nskip = isave(5)
568 3136 : head = isave(6)
569 3136 : col = isave(7)
570 3136 : itail = isave(8)
571 3136 : iter = isave(9)
572 3136 : iupdat = isave(10)
573 3136 : nseg = isave(12)
574 3136 : nfgv = isave(13)
575 3136 : info = isave(14)
576 3136 : ifun = isave(15)
577 3136 : iword = isave(16)
578 3136 : nfree = isave(17)
579 3136 : nact = isave(18)
580 3136 : ileave = isave(19)
581 3136 : nenter = isave(20)
582 :
583 3136 : theta = dsave(1)
584 3136 : fold = dsave(2)
585 3136 : tol = dsave(3)
586 3136 : dnorm = dsave(4)
587 3136 : epsmch = dsave(5)
588 3136 : cpu1 = dsave(6)
589 3136 : cachyt = dsave(7)
590 3136 : sbtime = dsave(8)
591 3136 : lnscht = dsave(9)
592 3136 : time1 = dsave(10)
593 3136 : gd = dsave(11)
594 3136 : step_max = dsave(12)
595 3136 : g_inf_norm = dsave(13)
596 3136 : stp = dsave(14)
597 3136 : gdold = dsave(15)
598 3136 : dtd = dsave(16)
599 :
600 : ! After returning from the driver go to the point where execution
601 : ! is to resume.
602 :
603 3136 : IF (task(1:4) == 'STOP') THEN
604 0 : IF (task(7:9) == 'CPU') THEN
605 : ! restore the previous iterate.
606 0 : CALL dcopy(n, t, 1, x, 1)
607 0 : CALL dcopy(n, r, 1, g, 1)
608 0 : f = fold
609 : END IF
610 0 : CALL timer(time2)
611 0 : time = time2 - time1
612 : CALL prn3lb(n, x, f, task, iprint, info, itfile, &
613 : iter, nfgv, nintol, nskip, nact, g_inf_norm, &
614 : time, nseg, word, iback, stp, xstep, k, &
615 0 : cachyt, sbtime, lnscht, wunit)
616 : CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
617 : iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
618 0 : cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
619 0 : RETURN
620 : END IF
621 : END IF
622 :
623 3136 : IF (.NOT. (task(1:5) == 'FG_LN' .OR. task(1:5) == 'NEW_X')) THEN
624 :
625 : ! Compute f0 and g0.
626 44 : nfgv = 1
627 :
628 : ! Compute the infinity norm of the (-) projected gradient.
629 :
630 44 : CALL projgr(n, lower_bound, upper_bound, nbd, x, g, g_inf_norm)
631 :
632 44 : IF (iprint >= 1) THEN
633 44 : WRITE (wunit, 1002) iter, f, g_inf_norm
634 44 : WRITE (itfile, 1003) iter, nfgv, g_inf_norm, f
635 : END IF
636 44 : IF (g_inf_norm <= pgtol) THEN
637 : ! terminate the algorithm.
638 0 : task = 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
639 0 : CALL timer(time2)
640 0 : time = time2 - time1
641 : CALL prn3lb(n, x, f, task, iprint, info, itfile, &
642 : iter, nfgv, nintol, nskip, nact, g_inf_norm, &
643 : time, nseg, word, iback, stp, xstep, k, &
644 0 : cachyt, sbtime, lnscht, wunit)
645 : CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
646 : iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
647 0 : cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
648 0 : RETURN
649 : END IF
650 : END IF
651 :
652 : first = .TRUE.
653 : DO WHILE (.TRUE.)
654 4576 : IF (.NOT. first .OR. .NOT. (task(1:5) == 'FG_LN' .OR. task(1:5) == 'NEW_X')) THEN
655 1484 : IF (iprint >= 99) WRITE (wunit, 1001) iter + 1
656 1484 : iword = -1
657 : !
658 1484 : IF (.NOT. constrained .AND. col > 0) THEN
659 : ! skip the search for GCP.
660 1423 : CALL dcopy(n, x, 1, z, 1)
661 1423 : wrk = updatd
662 1423 : nseg = 0
663 : ELSE
664 :
665 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
666 : !
667 : ! Compute the Generalized Cauchy Point (GCP).
668 : !
669 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
670 :
671 61 : CALL timer(cpu1)
672 : CALL cauchy(n, x, lower_bound, upper_bound, nbd, g, indx2, iwhere, t, d, z, &
673 : m, wy, ws, sy, wt, theta, col, head, &
674 : wa(1), wa(2*m + 1), wa(4*m + 1), wa(6*m + 1), nseg, &
675 61 : iprint, g_inf_norm, info, epsmch, wunit)
676 : ! applies rotation matrices to coordinates
677 61 : IF (keep_space_group) THEN
678 1 : CALL spgr_apply_rotations_coord(spgr, z)
679 : END IF
680 61 : IF (info /= 0) THEN
681 : ! singular triangular system detected; refresh the lbfgs memory.
682 0 : IF (iprint >= 1) WRITE (wunit, 1005)
683 0 : info = 0
684 0 : col = 0
685 0 : head = 1
686 0 : theta = one
687 0 : iupdat = 0
688 0 : updatd = .FALSE.
689 0 : CALL timer(cpu2)
690 0 : cachyt = cachyt + cpu2 - cpu1
691 0 : first = .FALSE.
692 0 : CYCLE
693 : END IF
694 61 : CALL timer(cpu2)
695 61 : cachyt = cachyt + cpu2 - cpu1
696 61 : nintol = nintol + nseg
697 :
698 : ! Count the entering and leaving variables for iter > 0;
699 : ! find the index set of free and active variables at the GCP.
700 :
701 : CALL freev(n, nfree, index, nenter, ileave, indx2, &
702 61 : iwhere, wrk, updatd, constrained, iprint, iter, wunit)
703 61 : nact = n - nfree
704 :
705 : END IF
706 :
707 : ! If there are no free variables or B=theta*I, then
708 : ! skip the subspace minimization.
709 :
710 1484 : IF (.NOT. (nfree == 0 .OR. col == 0)) THEN
711 :
712 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
713 : !
714 : ! Subspace minimization.
715 : !
716 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
717 :
718 1436 : CALL timer(cpu1)
719 :
720 : ! Form the LEL^T factorization of the indefinite
721 : ! matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
722 : ! [L_a -R_z theta*S'AA'S ]
723 : ! where E = [-I 0]
724 : ! [ 0 I]
725 :
726 1436 : IF (wrk) CALL formk(n, nfree, index, nenter, ileave, indx2, iupdat, &
727 1436 : updatd, wn, snd, m, ws, wy, sy, theta, col, head, info)
728 1436 : IF (info /= 0) THEN
729 : ! nonpositive definiteness in Cholesky factorization;
730 : ! refresh the lbfgs memory and restart the iteration.
731 0 : IF (iprint >= 1) WRITE (wunit, 1006)
732 0 : info = 0
733 0 : col = 0
734 0 : head = 1
735 0 : theta = one
736 0 : iupdat = 0
737 0 : updatd = .FALSE.
738 0 : CALL timer(cpu2)
739 0 : sbtime = sbtime + cpu2 - cpu1
740 0 : first = .FALSE.
741 0 : CYCLE
742 : END IF
743 :
744 : ! compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x)
745 : ! from 'cauchy').
746 : CALL cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index, &
747 1436 : theta, col, head, nfree, constrained, info)
748 : ! applies rotation matrices to coordinates
749 1436 : IF (keep_space_group) THEN
750 0 : CALL spgr_apply_rotations_force(spgr, r)
751 : END IF
752 1436 : IF (info == 0) THEN
753 :
754 : ! call the direct method.
755 :
756 : CALL subsm(n, m, nfree, index, lower_bound, upper_bound, nbd, z, r, xp, ws, wy, &
757 1436 : theta, x, g, col, head, iword, wa, wn, iprint, info, wunit)
758 : ! applies rotation matrices to coordinates
759 1436 : IF (keep_space_group) THEN
760 0 : CALL spgr_apply_rotations_coord(spgr, z)
761 0 : CALL spgr_apply_rotations_force(spgr, r)
762 : END IF
763 : END IF
764 1436 : IF (info /= 0) THEN
765 : ! singular triangular system detected;
766 : ! refresh the lbfgs memory and restart the iteration.
767 0 : IF (iprint >= 1) WRITE (wunit, 1005)
768 0 : info = 0
769 0 : col = 0
770 0 : head = 1
771 0 : theta = one
772 0 : iupdat = 0
773 0 : updatd = .FALSE.
774 0 : CALL timer(cpu2)
775 0 : sbtime = sbtime + cpu2 - cpu1
776 0 : first = .FALSE.
777 0 : CYCLE
778 : END IF
779 :
780 1436 : CALL timer(cpu2)
781 1436 : sbtime = sbtime + cpu2 - cpu1
782 : END IF
783 :
784 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
785 : !
786 : ! Line search and optimality tests.
787 : !
788 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
789 :
790 : ! Generate the search direction d:=z-x.
791 : ! applies rotation matrices to coordinates
792 1484 : IF (keep_space_group) THEN
793 1 : CALL spgr_apply_rotations_coord(spgr, x)
794 1 : CALL spgr_apply_rotations_coord(spgr, z)
795 : END IF
796 866144 : DO i = 1, n
797 866144 : d(i) = z(i) - x(i)
798 : END DO
799 1484 : CALL timer(cpu1)
800 : END IF
801 4576 : IF (.NOT. first .OR. .NOT. (task(1:5) == 'NEW_X')) THEN
802 : ! applies rotation matrices to coordinates
803 3135 : IF (keep_space_group) THEN
804 2 : CALL spgr_apply_rotations_coord(spgr, x)
805 2 : CALL spgr_apply_rotations_coord(spgr, z)
806 2 : CALL spgr_apply_rotations_force(spgr, d)
807 2 : CALL spgr_apply_rotations_force(spgr, g)
808 2 : CALL spgr_apply_rotations_force(spgr, r)
809 : END IF
810 : CALL lnsrlb(n, lower_bound, upper_bound, nbd, x, f, fold, gd, gdold, g, d, r, t, z, stp, dnorm, &
811 : dtd, xstep, step_max, iter, ifun, iback, nfgv, info, task, &
812 3135 : boxed, constrained, csave, isave(22), dsave(17), wunit)
813 : ! applies rotation matrices to coordinates
814 3135 : IF (keep_space_group) THEN
815 2 : CALL spgr_apply_rotations_coord(spgr, x)
816 2 : CALL spgr_apply_rotations_force(spgr, g)
817 : END IF
818 3135 : IF (info /= 0 .OR. iback >= 20) THEN
819 : ! restore the previous iterate.
820 0 : CALL dcopy(n, t, 1, x, 1)
821 0 : CALL dcopy(n, r, 1, g, 1)
822 0 : f = fold
823 0 : IF (col == 0) THEN
824 : ! abnormal termination.
825 0 : IF (info == 0) THEN
826 0 : info = -9
827 : ! restore the actual number of f and g evaluations etc.
828 0 : nfgv = nfgv - 1
829 0 : ifun = ifun - 1
830 0 : iback = iback - 1
831 : END IF
832 0 : task = 'ABNORMAL_TERMINATION_IN_LNSRCH'
833 0 : iter = iter + 1
834 0 : CALL timer(time2)
835 0 : time = time2 - time1
836 : CALL prn3lb(n, x, f, task, iprint, info, itfile, &
837 : iter, nfgv, nintol, nskip, nact, g_inf_norm, &
838 : time, nseg, word, iback, stp, xstep, k, &
839 0 : cachyt, sbtime, lnscht, wunit)
840 : CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
841 : iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
842 0 : cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
843 0 : RETURN
844 : ELSE
845 : ! refresh the lbfgs memory and restart the iteration.
846 0 : IF (iprint >= 1) WRITE (wunit, 1008)
847 0 : IF (info == 0) nfgv = nfgv - 1
848 0 : info = 0
849 0 : col = 0
850 0 : head = 1
851 0 : theta = one
852 0 : iupdat = 0
853 0 : updatd = .FALSE.
854 0 : task = 'RESTART_FROM_LNSRCH'
855 0 : CALL timer(cpu2)
856 0 : lnscht = lnscht + cpu2 - cpu1
857 0 : first = .FALSE.
858 0 : CYCLE
859 : END IF
860 3135 : ELSE IF (task(1:5) == 'FG_LN') THEN
861 : ! return to the driver for calculating f and g; reenter at 666.
862 : CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
863 : iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
864 1651 : cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
865 1651 : RETURN
866 : ELSE
867 : ! calculate and print out the quantities related to the new X.
868 1484 : CALL timer(cpu2)
869 1484 : lnscht = lnscht + cpu2 - cpu1
870 1484 : iter = iter + 1
871 :
872 : ! Compute the infinity norm of the projected (-)gradient.
873 :
874 1484 : CALL projgr(n, lower_bound, upper_bound, nbd, x, g, g_inf_norm)
875 :
876 : ! Print iteration information.
877 :
878 : CALL prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact, &
879 1484 : g_inf_norm, nseg, word, iword, iback, stp, xstep, wunit)
880 : CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
881 : iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
882 1484 : cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
883 2925 : RETURN
884 : END IF
885 : END IF
886 :
887 : ! Test for termination.
888 :
889 1441 : IF (g_inf_norm <= pgtol) THEN
890 : ! terminate the algorithm.
891 0 : task = 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
892 0 : CALL timer(time2)
893 0 : time = time2 - time1
894 : CALL prn3lb(n, x, f, task, iprint, info, itfile, &
895 : iter, nfgv, nintol, nskip, nact, g_inf_norm, &
896 : time, nseg, word, iback, stp, xstep, k, &
897 0 : cachyt, sbtime, lnscht, wunit)
898 : CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
899 : iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
900 0 : cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
901 0 : RETURN
902 : END IF
903 :
904 1441 : ddum = MAX(ABS(fold), ABS(f), one)
905 1441 : IF ((fold - f) <= tol*ddum) THEN
906 : ! terminate the algorithm.
907 1 : task = 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
908 1 : IF (iback >= 10) info = -5
909 : ! i.e., to issue a warning if iback>10 in the line search.
910 1 : CALL timer(time2)
911 1 : time = time2 - time1
912 : CALL prn3lb(n, x, f, task, iprint, info, itfile, &
913 : iter, nfgv, nintol, nskip, nact, g_inf_norm, &
914 : time, nseg, word, iback, stp, xstep, k, &
915 1 : cachyt, sbtime, lnscht, wunit)
916 : CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
917 : iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
918 1 : cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
919 1 : RETURN
920 : END IF
921 :
922 : ! Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's.
923 1440 : IF (keep_space_group) THEN
924 0 : CALL spgr_apply_rotations_force(spgr, g)
925 0 : CALL spgr_apply_rotations_force(spgr, r)
926 : END IF
927 837525 : DO i = 1, n
928 837525 : r(i) = g(i) - r(i)
929 : END DO
930 1440 : rr = ddot(n, r, 1, r, 1)
931 1440 : IF (stp == one) THEN
932 1306 : dr = gd - gdold
933 1306 : ddum = -gdold
934 : ELSE
935 134 : dr = (gd - gdold)*stp
936 134 : CALL dscal(n, stp, d, 1)
937 134 : ddum = -gdold*stp
938 : END IF
939 :
940 1440 : IF (dr <= epsmch*ddum) THEN
941 : ! skip the L-BFGS update.
942 4 : nskip = nskip + 1
943 4 : updatd = .FALSE.
944 4 : IF (iprint >= 1) WRITE (wunit, 1004) dr, ddum
945 : first = .FALSE.
946 : CYCLE
947 : END IF
948 :
949 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
950 : !
951 : ! Update the L-BFGS matrix.
952 : !
953 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
954 :
955 1436 : updatd = .TRUE.
956 1436 : iupdat = iupdat + 1
957 :
958 : ! Update matrices WS and WY and form the middle matrix in B.
959 :
960 : CALL matupd(n, m, ws, wy, sy, ss, d, r, itail, &
961 1436 : iupdat, col, head, theta, rr, dr, stp, dtd)
962 :
963 : ! Form the upper half of the pds T = theta*SS + L*D^(-1)*L';
964 : ! Store T in the upper triangular of the array wt;
965 : ! Cholesky factorize T to J*J' with
966 : ! J' stored in the upper triangular of wt.
967 :
968 1436 : CALL formt(m, wt, sy, ss, col, theta, info)
969 :
970 1436 : IF (info /= 0) THEN
971 : ! nonpositive definiteness in Cholesky factorization;
972 : ! refresh the lbfgs memory and restart the iteration.
973 0 : IF (iprint >= 1) WRITE (wunit, 1007)
974 0 : info = 0
975 0 : col = 0
976 0 : head = 1
977 0 : theta = one
978 0 : iupdat = 0
979 0 : updatd = .FALSE.
980 : END IF
981 :
982 : ! Now the inverse of the middle matrix in B is
983 :
984 : ! [ D^(1/2) O ] [ -D^(1/2) D^(-1/2)*L' ]
985 : ! [ -L*D^(-1/2) J ] [ 0 J' ]
986 :
987 : first = .FALSE.
988 : END DO
989 :
990 : 1001 FORMAT(//, ' L-BFGS| ITERATION ', i5)
991 : 1002 FORMAT &
992 : (/, ' L-BFGS| At iterate', i5, 4x, 'f= ', 1p, d12.5, 4x, '|proj g|= ', 1p, d12.5)
993 : 1003 FORMAT(2(1x, i4), 5x, '-', 5x, '-', 3x, '-', 5x, '-', 5x, '-', 8x, '-', 3x, &
994 : 1p, 2(1x, d10.3))
995 : 1004 FORMAT(' L-BFGS| ys=', 1p, e10.3, ' -gs=', 1p, e10.3, ' BFGS update SKIPPED')
996 : 1005 FORMAT(/, &
997 : ' L-BFGS| Singular triangular system detected;', /, &
998 : ' L-BFGS| refresh the lbfgs memory and restart the iteration.')
999 : 1006 FORMAT(/, &
1000 : ' L-BFGS| Nonpositive definiteness in Cholesky factorization in formk;', /, &
1001 : ' L-BFGS| refresh the lbfgs memory and restart the iteration.')
1002 : 1007 FORMAT(/, &
1003 : ' L-BFGS| Nonpositive definiteness in Cholesky factorization in formt;', /, &
1004 : ' L-BFGS| refresh the lbfgs memory and restart the iteration.')
1005 : 1008 FORMAT(/, &
1006 : ' L-BFGS| Bad direction in the line search;', /, &
1007 : ' L-BFGS| refresh the lbfgs memory and restart the iteration.')
1008 :
1009 : RETURN
1010 :
1011 : END SUBROUTINE mainlb
1012 :
1013 : ! **************************************************************************************************
1014 : !> \brief This subroutine initializes iwhere and projects the initial x to the feasible set if necessary.
1015 : !> \param n ...
1016 : !> \param lower_bound the lower bound on x.
1017 : !> \param upper_bound the upper bound on x.
1018 : !> \param nbd ...
1019 : !> \param x ...
1020 : !> \param iwhere iwhere(i)=-1 if x(i) has no bounds
1021 : !> 3 if l(i)=u(i)
1022 : !> 0 otherwise.
1023 : !> In cauchy, iwhere is given finer gradations.
1024 : !> \param iprint ...
1025 : !> \param x_projected ...
1026 : !> \param constrained ...
1027 : !> \param boxed ...
1028 : !> \param iwunit User-specified write unit, if not set then WRITE statements
1029 : !> write to default_output_unit by default
1030 : !> \author NEOS, November 1994. (Latest revision June 1996.)
1031 : !> Optimization Technology Center.
1032 : !> Argonne National Laboratory and Northwestern University.
1033 : !> Written by
1034 : !> Ciyou Zhu
1035 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
1036 : ! **************************************************************************************************
1037 45 : SUBROUTINE active(n, lower_bound, upper_bound, nbd, x, iwhere, iprint, &
1038 : x_projected, constrained, boxed, iwunit)
1039 :
1040 : INTEGER, INTENT(in) :: n
1041 : REAL(KIND=dp), INTENT(in) :: lower_bound(n), upper_bound(n)
1042 : INTEGER :: nbd(n)
1043 : REAL(KIND=dp) :: x(n)
1044 : INTEGER, INTENT(out) :: iwhere(n)
1045 : INTEGER :: iprint
1046 : LOGICAL :: x_projected, constrained, boxed
1047 : INTEGER, OPTIONAL :: iwunit
1048 :
1049 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
1050 :
1051 : INTEGER :: i, nbdd, wunit
1052 :
1053 45 : wunit = default_output_unit
1054 45 : IF (PRESENT(iwunit)) THEN
1055 45 : IF (iwunit > 0) wunit = iwunit
1056 : END IF
1057 :
1058 : ! ************
1059 : ! Initialize nbdd, x_projected, constrained and boxed.
1060 :
1061 45 : nbdd = 0
1062 45 : x_projected = .FALSE.
1063 45 : constrained = .FALSE.
1064 45 : boxed = .TRUE.
1065 :
1066 : ! Project the initial x to the easible set if necessary.
1067 :
1068 28632 : DO i = 1, n
1069 28632 : IF (nbd(i) > 0) THEN
1070 9018 : IF (nbd(i) <= 2 .AND. x(i) <= lower_bound(i)) THEN
1071 0 : IF (x(i) < lower_bound(i)) THEN
1072 0 : x_projected = .TRUE.
1073 0 : x(i) = lower_bound(i)
1074 : END IF
1075 0 : nbdd = nbdd + 1
1076 9018 : ELSE IF (nbd(i) >= 2 .AND. x(i) >= upper_bound(i)) THEN
1077 0 : IF (x(i) > upper_bound(i)) THEN
1078 0 : x_projected = .TRUE.
1079 0 : x(i) = upper_bound(i)
1080 : END IF
1081 0 : nbdd = nbdd + 1
1082 : END IF
1083 : END IF
1084 : END DO
1085 :
1086 : ! Initialize iwhere and assign values to constrained and boxed.
1087 :
1088 28632 : DO i = 1, n
1089 28587 : IF (nbd(i) /= 2) boxed = .FALSE.
1090 28632 : IF (nbd(i) == 0) THEN
1091 : ! this variable is always free
1092 19569 : iwhere(i) = -1
1093 :
1094 : ! otherwise set x(i)=mid(x(i), u(i), l(i)).
1095 : ELSE
1096 9018 : constrained = .TRUE.
1097 9018 : IF (nbd(i) == 2 .AND. upper_bound(i) - lower_bound(i) <= zero) THEN
1098 : ! this variable is always fixed
1099 0 : iwhere(i) = 3
1100 : ELSE
1101 9018 : iwhere(i) = 0
1102 : END IF
1103 : END IF
1104 : END DO
1105 :
1106 45 : IF (iprint >= 0) THEN
1107 45 : IF (x_projected) WRITE (wunit, 2001)
1108 45 : IF (.NOT. constrained) WRITE (wunit, 3001)
1109 : END IF
1110 :
1111 45 : IF (iprint > 0) WRITE (wunit, 1001) nbdd
1112 :
1113 : 1001 FORMAT(/, ' L-BFGS| At X0 ', i9, ' variables are exactly at the bounds')
1114 : 2001 FORMAT(' L-BFGS| The initial X is infeasible. Restart with its projection.')
1115 : 3001 FORMAT(' L-BFGS| This problem is unconstrained.')
1116 :
1117 45 : RETURN
1118 :
1119 : END SUBROUTINE active
1120 :
1121 : ! **************************************************************************************************
1122 : !> \brief This subroutine computes the product of the 2m x 2m middle matrix
1123 : !> in the compact L-BFGS formula of B and a 2m vector v;
1124 : !> it returns the product in p.
1125 : !> \param m m is the maximum number of variable metric corrections
1126 : !> used to define the limited memory matrix.
1127 : !> \param sy sy specifies the matrix S'Y.
1128 : !> \param wt wt specifies the upper triangular matrix J' which is
1129 : !> the Cholesky factor of (thetaS'S+LD^(-1)L').
1130 : !> \param col col specifies the number of s-vectors (or y-vectors)
1131 : !> stored in the compact L-BFGS formula.
1132 : !> \param v v specifies vector v.
1133 : !> \param p p is the product Mv.
1134 : !> \param info info = 0 for normal return,
1135 : !> = nonzero for abnormal return when the system to be solved by dtrsl is singular.
1136 : !> \author NEOS, November 1994. (Latest revision June 1996.)
1137 : !> Optimization Technology Center.
1138 : !> Argonne National Laboratory and Northwestern University.
1139 : !> Written by
1140 : !> Ciyou Zhu
1141 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
1142 : ! **************************************************************************************************
1143 34 : SUBROUTINE bmv(m, sy, wt, col, v, p, info)
1144 :
1145 : INTEGER :: m
1146 : REAL(KIND=dp) :: sy(m, m), wt(m, m)
1147 : INTEGER :: col
1148 : REAL(KIND=dp), INTENT(in) :: v(2*col)
1149 : REAL(KIND=dp), INTENT(out) :: p(2*col)
1150 : INTEGER, INTENT(out) :: info
1151 :
1152 : INTEGER :: i, i2, k
1153 : REAL(KIND=dp) :: sum
1154 :
1155 34 : IF (col == 0) RETURN
1156 :
1157 : ! PART I: solve [ D^(1/2) O ] [ p1 ] = [ v1 ]
1158 : ! [ -L*D^(-1/2) J ] [ p2 ] [ v2 ].
1159 :
1160 : ! solve Jp2=v2+LD^(-1)v1.
1161 34 : p(col + 1) = v(col + 1)
1162 88 : DO i = 2, col
1163 54 : i2 = col + i
1164 54 : sum = 0.0_dp
1165 160 : DO k = 1, i - 1
1166 160 : sum = sum + sy(i, k)*v(k)/sy(k, k)
1167 : END DO
1168 88 : p(i2) = v(i2) + sum
1169 : END DO
1170 : ! Solve the triangular system
1171 34 : CALL dtrsl(wt, m, col, p(col + 1), 11, info)
1172 34 : IF (info /= 0) RETURN
1173 :
1174 : ! solve D^(1/2)p1=v1.
1175 122 : DO i = 1, col
1176 122 : p(i) = v(i)/SQRT(sy(i, i))
1177 : END DO
1178 :
1179 : ! PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ]
1180 : ! [ 0 J' ] [ p2 ] [ p2 ].
1181 :
1182 : ! solve J^Tp2=p2.
1183 34 : CALL dtrsl(wt, m, col, p(col + 1), 01, info)
1184 34 : IF (info /= 0) RETURN
1185 :
1186 : ! compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2)
1187 : ! =-D^(-1/2)p1+D^(-1)L'p2.
1188 122 : DO i = 1, col
1189 122 : p(i) = -p(i)/SQRT(sy(i, i))
1190 : END DO
1191 122 : DO i = 1, col
1192 88 : sum = 0._dp
1193 194 : DO k = i + 1, col
1194 194 : sum = sum + sy(k, i)*p(col + k)/sy(i, i)
1195 : END DO
1196 122 : p(i) = p(i) + sum
1197 : END DO
1198 :
1199 : RETURN
1200 :
1201 : END SUBROUTINE bmv
1202 :
1203 : ! **************************************************************************************************
1204 : !> \brief For given x, l, u, g (with g_inf_norm > 0), and a limited memory
1205 : !> BFGS matrix B defined in terms of matrices WY, WS, WT, and
1206 : !> scalars head, col, and theta, this subroutine computes the
1207 : !> generalized Cauchy point (GCP), defined as the first local
1208 : !> minimizer of the quadratic
1209 : !>
1210 : !> Q(x + s) = g's + 1/2 s'Bs
1211 : !>
1212 : !> along the projected gradient direction P(x-tg,l,u).
1213 : !> The routine returns the GCP in xcp.
1214 : !> \param n n is the dimension of the problem.
1215 : !> \param x x is the starting point for the GCP computation.
1216 : !> \param lower_bound the lower bound on x.
1217 : !> \param upper_bound the upper bound on x.
1218 : !> \param nbd nbd represents the type of bounds imposed on the
1219 : !> variables, and must be specified as follows:
1220 : !> nbd(i)=0 if x(i) is unbounded,
1221 : !> 1 if x(i) has only a lower bound,
1222 : !> 2 if x(i) has both lower and upper bounds, and
1223 : !> 3 if x(i) has only an upper bound.
1224 : !> \param g g is the gradient of f(x). g must be a nonzero vector.
1225 : !> \param iorder iorder will be used to store the breakpoints in the piecewise
1226 : !> linear path and free variables encountered. On exit,
1227 : !> iorder(1),...,iorder(nleft) are indices of breakpoints
1228 : !> which have not been encountered;
1229 : !> iorder(nleft+1),...,iorder(nbreak) are indices of
1230 : !> encountered breakpoints; and
1231 : !> iorder(nfree),...,iorder(n) are indices of variables which
1232 : !> have no bound constraits along the search direction.
1233 : !> \param iwhere On entry iwhere indicates only the permanently fixed (iwhere=3)
1234 : !> or free (iwhere= -1) components of x.
1235 : !> On exit iwhere records the status of the current x variables.
1236 : !> iwhere(i)=-3 if x(i) is free and has bounds, but is not moved
1237 : !> 0 if x(i) is free and has bounds, and is moved
1238 : !> 1 if x(i) is fixed at l(i), and l(i) .ne. u(i)
1239 : !> 2 if x(i) is fixed at u(i), and u(i) .ne. l(i)
1240 : !> 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i)
1241 : !> -1 if x(i) is always free, i.e., it has no bounds.
1242 : !> \param t t will be used to store the break points.
1243 : !> \param d d is used to store the Cauchy direction P(x-tg)-x.
1244 : !> \param xcp is a double precision array of dimension n used to return the GCP on exit.
1245 : !> \param m m is the maximum number of variable metric corrections used to define the limited memory matrix.
1246 : !> \param wy ws, wy, sy, and wt are double precision arrays.
1247 : !> On entry they store information that defines the limited memory BFGS matrix:
1248 : !> wy(n,m) stores Y, a set of y-vectors;
1249 : !> \param ws ws(n,m) stores S, a set of s-vectors;
1250 : !> \param sy sy(m,m) stores S'Y;
1251 : !> \param wt wt(m,m) stores the Cholesky factorization of (theta*S'S+LD^(-1)L').
1252 : !> \param theta theta is the scaling factor specifying B_0 = theta I.
1253 : !> \param col col is the actual number of variable metric corrections stored so far.
1254 : !> \param head head is the location of the first s-vector (or y-vector in S (or Y)
1255 : !> \param p p will be used to store the vector p = W^(T)d.
1256 : !> \param c c will be used to store the vector c = W^(T)(xcp-x).
1257 : !> \param wbp wbp will be used to store the row of W corresponding to a breakpoint.
1258 : !> \param v v is a double precision working array.
1259 : !> \param nseg On exit nseg records the number of quadratic segments explored in searching for the GCP.
1260 : !> \param iprint iprint is an INTEGER variable that must be set by the user.
1261 : !> It controls the frequency and type of output generated:
1262 : !> iprint<0 no output is generated;
1263 : !> iprint=0 print only one line at the last iteration;
1264 : !> 0<iprint<99 print also f and |proj g| every iprint iterations;
1265 : !> iprint=99 print details of every iteration except n-vectors;
1266 : !> iprint=100 print also the changes of active set and final x;
1267 : !> iprint>100 print details of every iteration including x and g;
1268 : !> When iprint > 0, the file iterate.dat will be created to summarize the iteration.
1269 : !> \param g_inf_norm g_inf_norm is the norm of the projected gradient at x.
1270 : !> \param info On entry info is 0.
1271 : !> On exit info = 0 for normal return,
1272 : !> = nonzero for abnormal return when the the system
1273 : !> used in routine bmv is singular.
1274 : !> \param epsmch ...
1275 : !> \param iwunit User-specified write unit, if not set then WRITE statements
1276 : !> write to default_output_unit by default
1277 : !> \author NEOS, November 1994. (Latest revision June 1996.)
1278 : !> Optimization Technology Center.
1279 : !> Argonne National Laboratory and Northwestern University.
1280 : !> Written by
1281 : !> Ciyou Zhu
1282 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
1283 : ! **************************************************************************************************
1284 61 : SUBROUTINE cauchy(n, x, lower_bound, upper_bound, nbd, g, iorder, iwhere, t, d, xcp, &
1285 61 : m, wy, ws, sy, wt, theta, col, head, p, c, wbp, &
1286 61 : v, nseg, iprint, g_inf_norm, info, epsmch, iwunit)
1287 : INTEGER, INTENT(in) :: n
1288 : REAL(KIND=dp), INTENT(in) :: x(n), lower_bound(n), upper_bound(n)
1289 : INTEGER, INTENT(in) :: nbd(n)
1290 : REAL(KIND=dp), INTENT(in) :: g(n)
1291 : INTEGER :: iorder(n)
1292 : INTEGER, INTENT(inout) :: iwhere(n)
1293 : REAL(KIND=dp) :: t(n), d(n), xcp(n)
1294 : INTEGER, INTENT(in) :: m
1295 : REAL(KIND=dp), INTENT(in) :: sy(m, m), wt(m, m), theta
1296 : INTEGER, INTENT(in) :: col
1297 : REAL(KIND=dp), INTENT(in) :: ws(n, col), wy(n, col)
1298 : INTEGER, INTENT(in) :: head
1299 : REAL(KIND=dp) :: p(2*m), c(2*m), wbp(2*m), v(2*m)
1300 : INTEGER :: nseg, iprint
1301 : REAL(KIND=dp), INTENT(in) :: g_inf_norm
1302 : INTEGER, INTENT(inout) :: info
1303 : REAL(KIND=dp) :: epsmch
1304 : INTEGER, OPTIONAL :: iwunit
1305 :
1306 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
1307 :
1308 : INTEGER :: col2, i, ibkmin, ibp, iter, j, nbreak, &
1309 : nfree, nleft, pointr, wunit
1310 : LOGICAL :: bnded, xlower, xupper
1311 : REAL(KIND=dp) :: bkmin, ddot, dibp, dibp2, dt, dtm, f1, &
1312 : f2, f2_org, neggi, tj, tj0, tl, tsum, &
1313 : tu, wmc, wmp, wmw, zibp
1314 :
1315 61 : wunit = default_output_unit
1316 61 : IF (PRESENT(iwunit)) THEN
1317 61 : IF (iwunit > 0) wunit = iwunit
1318 : END IF
1319 :
1320 : ! References:
1321 : !
1322 : ! [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
1323 : ! memory algorithm for bound constrained optimization'',
1324 : ! SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
1325 : !
1326 : ! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
1327 : ! Subroutines for Large Scale Bound Constrained Optimization''
1328 : ! Tech. Report, NAM-11, EECS Department, Northwestern University,
1329 : ! 1994.
1330 : !
1331 : ! (Postscript files of these papers are available via anonymous
1332 : ! ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
1333 : !
1334 : ! * * *
1335 : ! Check the status of the variables, reset iwhere(i) if necessary;
1336 : ! compute the Cauchy direction d and the breakpoints t; initialize
1337 : ! the derivative f1 and the vector p = W'd (for theta = 1).
1338 :
1339 61 : IF (g_inf_norm <= zero) THEN
1340 0 : IF (iprint >= 0) WRITE (wunit, 7010)
1341 0 : CALL dcopy(n, x, 1, xcp, 1)
1342 0 : RETURN
1343 : END IF
1344 61 : bnded = .TRUE.
1345 61 : nfree = n + 1
1346 61 : nbreak = 0
1347 61 : ibkmin = 0
1348 61 : bkmin = zero
1349 61 : col2 = 2*col
1350 61 : f1 = zero
1351 61 : IF (iprint >= 99) WRITE (wunit, 3010)
1352 :
1353 : ! We set p to zero and build it up as we determine d.
1354 :
1355 137 : DO i = 1, col2
1356 137 : p(i) = zero
1357 : END DO
1358 :
1359 : ! In the following loop we determine for each variable its bound
1360 : ! status and its breakpoint, and update p accordingly.
1361 : ! Smallest breakpoint is identified.
1362 :
1363 55738 : DO i = 1, n
1364 55677 : neggi = -g(i)
1365 55677 : IF (iwhere(i) /= 3 .AND. iwhere(i) /= -1) THEN
1366 : ! if x(i) is not a constant and has bounds,
1367 : ! compute the difference between x(i) and its bounds.
1368 36120 : IF (nbd(i) <= 2) tl = x(i) - lower_bound(i)
1369 36120 : IF (nbd(i) >= 2) tu = upper_bound(i) - x(i)
1370 :
1371 : ! If a variable is close enough to a bound
1372 : ! we treat it as at bound.
1373 36120 : xlower = nbd(i) <= 2 .AND. tl <= zero
1374 36120 : xupper = nbd(i) >= 2 .AND. tu <= zero
1375 :
1376 : ! reset iwhere(i).
1377 36120 : iwhere(i) = 0
1378 36120 : IF (xlower) THEN
1379 0 : IF (neggi <= zero) iwhere(i) = 1
1380 36120 : ELSE IF (xupper) THEN
1381 0 : IF (neggi >= zero) iwhere(i) = 2
1382 : ELSE
1383 36120 : IF (ABS(neggi) <= zero) iwhere(i) = -3
1384 : END IF
1385 : END IF
1386 55677 : pointr = head
1387 55738 : IF (iwhere(i) /= 0 .AND. iwhere(i) /= -1) THEN
1388 15 : d(i) = zero
1389 : ELSE
1390 55662 : d(i) = neggi
1391 55662 : f1 = f1 - neggi*neggi
1392 : ! calculate p := p - W'e_i* (g_i).
1393 114360 : DO j = 1, col
1394 58698 : p(j) = p(j) + wy(i, pointr)*neggi
1395 58698 : p(col + j) = p(col + j) + ws(i, pointr)*neggi
1396 114360 : pointr = MOD(pointr, m) + 1
1397 : END DO
1398 : IF (nbd(i) <= 2 .AND. nbd(i) /= 0 &
1399 55662 : & .AND. neggi < zero) THEN
1400 : ! x(i) + d(i) is bounded; compute t(i).
1401 18110 : nbreak = nbreak + 1
1402 18110 : iorder(nbreak) = i
1403 18110 : t(nbreak) = tl/(-neggi)
1404 18110 : IF (nbreak == 1 .OR. t(nbreak) < bkmin) THEN
1405 64 : bkmin = t(nbreak)
1406 64 : ibkmin = nbreak
1407 : END IF
1408 37552 : ELSE IF (nbd(i) >= 2 .AND. neggi > zero) THEN
1409 : ! x(i) + d(i) is bounded; compute t(i).
1410 17995 : nbreak = nbreak + 1
1411 17995 : iorder(nbreak) = i
1412 17995 : t(nbreak) = tu/neggi
1413 17995 : IF (nbreak == 1 .OR. t(nbreak) < bkmin) THEN
1414 49 : bkmin = t(nbreak)
1415 49 : ibkmin = nbreak
1416 : END IF
1417 : ELSE
1418 : ! x(i) + d(i) is not bounded.
1419 19557 : nfree = nfree - 1
1420 19557 : iorder(nfree) = i
1421 19557 : IF (ABS(neggi) > zero) bnded = .FALSE.
1422 : END IF
1423 : END IF
1424 : END DO
1425 :
1426 : ! The indices of the nonzero components of d are now stored
1427 : ! in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n).
1428 : ! The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin.
1429 :
1430 61 : IF (theta /= one) THEN
1431 : ! complete the initialization of p for theta not= one.
1432 13 : CALL dscal(col, theta, p(col + 1), 1)
1433 : END IF
1434 :
1435 : ! Initialize GCP xcp = x.
1436 :
1437 61 : CALL dcopy(n, x, 1, xcp, 1)
1438 :
1439 61 : IF (nbreak == 0 .AND. nfree == n + 1) THEN
1440 : ! is a zero vector, return with the initial xcp as GCP.
1441 0 : IF (iprint > 100) WRITE (wunit, 1010) (xcp(i), i=1, n)
1442 0 : RETURN
1443 : END IF
1444 :
1445 : ! Initialize c = W'(xcp - x) = 0.
1446 :
1447 137 : DO j = 1, col2
1448 137 : c(j) = zero
1449 : END DO
1450 :
1451 : ! Initialize derivative f2.
1452 :
1453 61 : f2 = -theta*f1
1454 61 : f2_org = f2
1455 61 : IF (col > 0) THEN
1456 13 : CALL bmv(m, sy, wt, col, p, v, info)
1457 13 : IF (info /= 0) RETURN
1458 13 : f2 = f2 - ddot(col2, v, 1, p, 1)
1459 : END IF
1460 61 : dtm = -f1/f2
1461 61 : tsum = zero
1462 61 : nseg = 1
1463 61 : IF (iprint >= 99) THEN
1464 0 : WRITE (wunit, 1011) nbreak
1465 : END IF
1466 :
1467 61 : nleft = nbreak
1468 61 : iter = 1
1469 :
1470 61 : tj = zero
1471 :
1472 : ! If there are no breakpoints, locate the GCP and return.
1473 :
1474 61 : IF (nleft == 0) THEN
1475 41 : IF (iprint >= 99) THEN
1476 0 : WRITE (wunit, 4012)
1477 0 : WRITE (wunit, 4010) nseg, f1, f2
1478 0 : WRITE (wunit, 6010) dtm
1479 : END IF
1480 41 : IF (dtm <= zero) dtm = zero
1481 41 : tsum = tsum + dtm
1482 :
1483 : ! Move free variables (i.e., the ones w/o breakpoints) and
1484 : ! the variables whose breakpoints haven't been reached.
1485 :
1486 41 : CALL daxpy(n, tsum, d, 1, xcp, 1)
1487 : END IF
1488 :
1489 71 : DO WHILE (nleft > 0)
1490 :
1491 : ! Find the next smallest breakpoint;
1492 : ! compute dt = t(nleft) - t(nleft + 1).
1493 :
1494 30 : tj0 = tj
1495 30 : IF (iter == 1) THEN
1496 : ! Since we already have the smallest breakpoint we need not do
1497 : ! heapsort yet. Often only one breakpoint is used and the
1498 : ! cost of heapsort is avoided.
1499 20 : tj = bkmin
1500 20 : ibp = iorder(ibkmin)
1501 : ELSE
1502 10 : IF (iter == 2) THEN
1503 : ! Replace the already used smallest breakpoint with the
1504 : ! breakpoint numbered nbreak > nlast, before heapsort call.
1505 5 : IF (ibkmin /= nbreak) THEN
1506 4 : t(ibkmin) = t(nbreak)
1507 4 : iorder(ibkmin) = iorder(nbreak)
1508 : END IF
1509 : ! Update heap structure of breakpoints
1510 : ! (if iter=2, initialize heap).
1511 : END IF
1512 10 : CALL hpsolb(nleft, t, iorder, iter - 2)
1513 10 : tj = t(nleft)
1514 10 : ibp = iorder(nleft)
1515 : END IF
1516 :
1517 30 : dt = tj - tj0
1518 :
1519 30 : IF (dt /= zero .AND. iprint >= 100) THEN
1520 0 : WRITE (wunit, 4011) nseg, f1, f2
1521 0 : WRITE (wunit, 5010) dt
1522 0 : WRITE (wunit, 6010) dtm
1523 : END IF
1524 :
1525 : ! If a minimizer is within this interval, locate the GCP and return.
1526 :
1527 30 : IF (dtm < dt) THEN
1528 20 : IF (iprint >= 99) THEN
1529 0 : WRITE (wunit, 4012)
1530 0 : WRITE (wunit, 4010) nseg, f1, f2
1531 0 : WRITE (wunit, 6010) dtm
1532 : END IF
1533 20 : IF (dtm <= zero) dtm = zero
1534 20 : tsum = tsum + dtm
1535 :
1536 : ! Move free variables (i.e., the ones w/o breakpoints) and
1537 : ! the variables whose breakpoints haven't been reached.
1538 :
1539 20 : CALL daxpy(n, tsum, d, 1, xcp, 1)
1540 20 : EXIT
1541 : END IF
1542 :
1543 : ! Otherwise fix one variable and
1544 : ! reset the corresponding component of d to zero.
1545 :
1546 10 : tsum = tsum + dt
1547 10 : nleft = nleft - 1
1548 10 : iter = iter + 1
1549 10 : dibp = d(ibp)
1550 10 : d(ibp) = zero
1551 10 : IF (dibp > zero) THEN
1552 2 : zibp = upper_bound(ibp) - x(ibp)
1553 2 : xcp(ibp) = upper_bound(ibp)
1554 2 : iwhere(ibp) = 2
1555 : ELSE
1556 8 : zibp = lower_bound(ibp) - x(ibp)
1557 8 : xcp(ibp) = lower_bound(ibp)
1558 8 : iwhere(ibp) = 1
1559 : END IF
1560 10 : IF (iprint >= 100) WRITE (wunit, 8010) ibp
1561 10 : IF (nleft == 0 .AND. nbreak == n) THEN
1562 : ! all n variables are fixed,
1563 : ! return with xcp as GCP.
1564 0 : dtm = dt
1565 0 : EXIT
1566 : END IF
1567 :
1568 : ! Update the derivative information.
1569 :
1570 10 : nseg = nseg + 1
1571 10 : dibp2 = dibp**2
1572 :
1573 : ! Update f1 and f2.
1574 :
1575 : ! temporarily set f1 and f2 for col=0.
1576 10 : f1 = f1 + dt*f2 + dibp2 - theta*dibp*zibp
1577 10 : f2 = f2 - theta*dibp2
1578 :
1579 10 : IF (col > 0) THEN
1580 : ! update c = c + dt*p.
1581 8 : CALL daxpy(col2, dt, p, 1, c, 1)
1582 :
1583 : ! choose wbp,
1584 : ! the row of W corresponding to the breakpoint encountered.
1585 8 : pointr = head
1586 20 : DO j = 1, col
1587 12 : wbp(j) = wy(ibp, pointr)
1588 12 : wbp(col + j) = theta*ws(ibp, pointr)
1589 20 : pointr = MOD(pointr, m) + 1
1590 : END DO
1591 :
1592 : ! compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'.
1593 8 : CALL bmv(m, sy, wt, col, wbp, v, info)
1594 8 : IF (info /= 0) RETURN
1595 8 : wmc = ddot(col2, c, 1, v, 1)
1596 8 : wmp = ddot(col2, p, 1, v, 1)
1597 8 : wmw = ddot(col2, wbp, 1, v, 1)
1598 :
1599 : ! update p = p - dibp*wbp.
1600 8 : CALL daxpy(col2, -dibp, wbp, 1, p, 1)
1601 :
1602 : ! complete updating f1 and f2 while col > 0.
1603 8 : f1 = f1 + dibp*wmc
1604 8 : f2 = f2 + 2.0_dp*dibp*wmp - dibp2*wmw
1605 : END IF
1606 :
1607 10 : f2 = MAX(epsmch*f2_org, f2)
1608 51 : IF (nleft > 0) THEN
1609 10 : dtm = -f1/f2
1610 : CYCLE
1611 : ! to repeat the loop for unsearched intervals.
1612 : ELSE
1613 0 : IF (bnded) THEN
1614 0 : f1 = zero
1615 0 : f2 = zero
1616 0 : dtm = zero
1617 : ELSE
1618 0 : dtm = -f1/f2
1619 : END IF
1620 0 : IF (iprint >= 99) THEN
1621 0 : WRITE (wunit, 4012)
1622 0 : WRITE (wunit, 4010) nseg, f1, f2
1623 0 : WRITE (wunit, 6010) dtm
1624 : END IF
1625 0 : IF (dtm <= zero) dtm = zero
1626 0 : tsum = tsum + dtm
1627 :
1628 : ! Move free variables (i.e., the ones w/o breakpoints) and
1629 : ! the variables whose breakpoints haven't been reached.
1630 :
1631 0 : CALL daxpy(n, tsum, d, 1, xcp, 1)
1632 0 : EXIT
1633 : END IF
1634 : END DO
1635 :
1636 : ! Update c = c + dtm*p = W'(x^c - x)
1637 : ! which will be used in computing r = Z'(B(x^c - x) + g).
1638 :
1639 61 : IF (col > 0) CALL daxpy(col2, dtm, p, 1, c, 1)
1640 61 : IF (iprint > 100) WRITE (wunit, 1010) (xcp(i), i=1, n)
1641 61 : IF (iprint >= 99) WRITE (wunit, 2010)
1642 :
1643 : 1010 FORMAT(' L-BFGS| Cauchy X = ', /, (4x, 1p, 6(1x, d11.4)))
1644 : 1011 FORMAT(/, ' L-BFGS| There are ', i12, ' breakpoints ')
1645 : 2010 FORMAT(/, ' L-BFGS| ---------------- exit CAUCHY-----------------')
1646 : 3010 FORMAT(/, ' L-BFGS| ---------------- enter CAUCHY ---------------')
1647 : 4010 FORMAT(' L-BFGS| Piece ', i3, ' --f1, f2 at start point ', 1p, 2(1x, d11.4))
1648 : 4011 FORMAT(/, ' L-BFGS| Piece ', i3, ' --f1, f2 at start point ', &
1649 : 1p, 2(1x, d11.4))
1650 : 4012 FORMAT(/, ' L-BFGS| GCP found in this segment')
1651 : 5010 FORMAT(' L-BFGS| Distance to the next break point = ', 1p, d11.4)
1652 : 6010 FORMAT(' L-BFGS| Distance to the stationary point = ', 1p, d11.4)
1653 : 7010 FORMAT(' L-BFGS| Subgnorm = 0. GCP = X.')
1654 : 8010 FORMAT(' L-BFGS| Variable ', i12, ' is fixed.')
1655 :
1656 : RETURN
1657 :
1658 : END SUBROUTINE cauchy
1659 :
1660 : ! **************************************************************************************************
1661 : !> \brief This subroutine computes r=-Z'B(xcp-xk)-Z'g by using
1662 : !> wa(2m+1)=W'(xcp-x) from subroutine cauchy.
1663 : !> \param n ...
1664 : !> \param m ...
1665 : !> \param x ...
1666 : !> \param g ...
1667 : !> \param ws ...
1668 : !> \param wy ...
1669 : !> \param sy ...
1670 : !> \param wt ...
1671 : !> \param z ...
1672 : !> \param r ...
1673 : !> \param wa ...
1674 : !> \param index ...
1675 : !> \param theta ...
1676 : !> \param col ...
1677 : !> \param head ...
1678 : !> \param nfree ...
1679 : !> \param constrained ...
1680 : !> \param info ...
1681 : !> \author NEOS, November 1994. (Latest revision June 1996.)
1682 : !> Optimization Technology Center.
1683 : !> Argonne National Laboratory and Northwestern University.
1684 : !> Written by
1685 : !> Ciyou Zhu
1686 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
1687 : ! **************************************************************************************************
1688 1436 : SUBROUTINE cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index, &
1689 : theta, col, head, nfree, constrained, info)
1690 :
1691 : INTEGER, INTENT(in) :: n, m
1692 : REAL(KIND=dp), INTENT(in) :: x(n), g(n), ws(n, m), wy(n, m), &
1693 : sy(m, m), wt(m, m), z(n)
1694 : REAL(KIND=dp), INTENT(out) :: r(n), wa(4*m)
1695 : INTEGER, INTENT(in) :: INDEX(n)
1696 : REAL(KIND=dp), INTENT(in) :: theta
1697 : INTEGER, INTENT(in) :: col, head, nfree
1698 : LOGICAL, INTENT(in) :: constrained
1699 : INTEGER :: info
1700 :
1701 : INTEGER :: i, j, k, pointr
1702 : REAL(KIND=dp) :: a1, a2
1703 :
1704 1436 : IF (.NOT. constrained .AND. col > 0) THEN
1705 810406 : DO i = 1, n
1706 810406 : r(i) = -g(i)
1707 : END DO
1708 : ELSE
1709 27083 : DO i = 1, nfree
1710 27070 : k = INDEX(i)
1711 27083 : r(i) = -theta*(z(k) - x(k)) - g(k)
1712 : END DO
1713 13 : CALL bmv(m, sy, wt, col, wa(2*m + 1), wa(1), info)
1714 13 : IF (info /= 0) THEN
1715 0 : info = -8
1716 0 : RETURN
1717 : END IF
1718 13 : pointr = head
1719 51 : DO j = 1, col
1720 38 : a1 = wa(j)
1721 38 : a2 = theta*wa(col + j)
1722 58754 : DO i = 1, nfree
1723 58716 : k = INDEX(i)
1724 58754 : r(i) = r(i) + wy(k, pointr)*a1 + ws(k, pointr)*a2
1725 : END DO
1726 51 : pointr = MOD(pointr, m) + 1
1727 : END DO
1728 : END IF
1729 :
1730 : RETURN
1731 :
1732 : END SUBROUTINE cmprlb
1733 :
1734 : ! **************************************************************************************************
1735 : !> \brief This subroutine checks the validity of the input data.
1736 : !> \param n ...
1737 : !> \param m ...
1738 : !> \param factr ...
1739 : !> \param lower_bound the lower bound on x.
1740 : !> \param upper_bound the upper bound on x.
1741 : !> \param nbd ...
1742 : !> \param task ...
1743 : !> \param info ...
1744 : !> \param k ...
1745 : !> \author NEOS, November 1994. (Latest revision June 1996.)
1746 : !> Optimization Technology Center.
1747 : !> Argonne National Laboratory and Northwestern University.
1748 : !> Written by
1749 : !> Ciyou Zhu
1750 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
1751 : ! **************************************************************************************************
1752 45 : SUBROUTINE errclb(n, m, factr, lower_bound, upper_bound, nbd, task, info, k)
1753 :
1754 : INTEGER, INTENT(in) :: n, m
1755 : REAL(KIND=dp), INTENT(in) :: factr, lower_bound(n), upper_bound(n)
1756 : INTEGER :: nbd(n)
1757 : CHARACTER(LEN=60) :: task
1758 : INTEGER :: info, k
1759 :
1760 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
1761 :
1762 : INTEGER :: i
1763 :
1764 : ! Check the input arguments for errors.
1765 :
1766 45 : IF (n <= 0) task = 'ERROR: N <= 0'
1767 45 : IF (m <= 0) task = 'ERROR: M <= 0'
1768 45 : IF (factr < zero) task = 'ERROR: FACTR < 0'
1769 :
1770 : ! Check the validity of the arrays nbd(i), u(i), and l(i).
1771 :
1772 28632 : DO i = 1, n
1773 28587 : IF (nbd(i) < 0 .OR. nbd(i) > 3) THEN
1774 : ! return
1775 0 : task = 'ERROR: INVALID NBD'
1776 0 : info = -6
1777 0 : k = i
1778 : END IF
1779 28632 : IF (nbd(i) == 2) THEN
1780 9018 : IF (lower_bound(i) > upper_bound(i)) THEN
1781 : ! return
1782 0 : task = 'ERROR: NO FEASIBLE SOLUTION'
1783 0 : info = -7
1784 0 : k = i
1785 : END IF
1786 : END IF
1787 : END DO
1788 :
1789 45 : RETURN
1790 :
1791 : END SUBROUTINE errclb
1792 :
1793 : ! **************************************************************************************************
1794 : !> \brief This subroutine forms the LEL^T factorization of the indefinite
1795 : !> matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
1796 : !> [L_a -R_z theta*S'AA'S ]
1797 : !> where E = [-I 0]
1798 : !> [ 0 I]
1799 : !> The matrix K can be shown to be equal to the matrix M^[-1]N
1800 : !> occurring in section 5.1 of [1], as well as to the matrix
1801 : !> Mbar^[-1] Nbar in section 5.3.
1802 : !> \param n n is the dimension of the problem.
1803 : !> \param nsub nsub is the number of subspace variables in free set.
1804 : !> \param ind ind specifies the indices of subspace variables.
1805 : !> \param nenter nenter is the number of variables entering the free set.
1806 : !> \param ileave indx2(ileave),...,indx2(n) are the variables leaving the free set.
1807 : !> \param indx2 indx2(1),...,indx2(nenter) are the variables entering the free set,
1808 : !> while indx2(ileave),...,indx2(n) are the variables leaving the free set.
1809 : !> \param iupdat iupdat is the total number of BFGS updates made so far.
1810 : !> \param updatd 'updatd' is true if the L-BFGS matrix is updatd.
1811 : !> \param wn the upper triangle of wn stores the LEL^T factorization
1812 : !> of the 2*col x 2*col indefinite matrix
1813 : !> [-D -Y'ZZ'Y/theta L_a'-R_z' ]
1814 : !> [L_a -R_z theta*S'AA'S ]
1815 : !> \param wn1 On entry wn1 stores the lower triangular part of
1816 : !> [Y' ZZ'Y L_a'+R_z']
1817 : !> [L_a+R_z S'AA'S ]
1818 : !> in the previous iteration.
1819 : !> On exit wn1 stores the corresponding updated matrices.
1820 : !> The purpose of wn1 is just to store these inner products
1821 : !> so they can be easily updated and inserted into wn.
1822 : !> \param m m is the maximum number of variable metric corrections
1823 : !> used to define the limited memory matrix.
1824 : !> \param ws ws(n,m) stores S, a set of s-vectors;
1825 : !> \param wy wy(n,m) stores Y, a set of y-vectors;
1826 : !> \param sy sy(m,m) stores S'Y;
1827 : !> \param theta is the scaling factor specifying B_0 = theta I;
1828 : !> \param col is the number of variable metric corrections stored;
1829 : !> \param head is the location of the 1st s- (or y-) vector in S (or Y).
1830 : !> \param info info = 0 for normal return;
1831 : !> = -1 when the 1st Cholesky factorization failed;
1832 : !> = -2 when the 2st Cholesky factorization failed.
1833 : !> \author NEOS, November 1994. (Latest revision June 1996.)
1834 : !> Optimization Technology Center.
1835 : !> Argonne National Laboratory and Northwestern University.
1836 : !> Written by
1837 : !> Ciyou Zhu
1838 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
1839 : ! **************************************************************************************************
1840 1436 : SUBROUTINE formk(n, nsub, ind, nenter, ileave, indx2, iupdat, &
1841 1436 : updatd, wn, wn1, m, ws, wy, sy, theta, col, &
1842 : head, info)
1843 :
1844 : INTEGER, INTENT(in) :: n, nsub, ind(n), nenter, ileave, &
1845 : indx2(n), iupdat
1846 : LOGICAL :: updatd
1847 : INTEGER, INTENT(in) :: m
1848 : REAL(KIND=dp) :: wn1(2*m, 2*m)
1849 : REAL(KIND=dp), INTENT(out) :: wn(2*m, 2*m)
1850 : REAL(KIND=dp), INTENT(in) :: ws(n, m), wy(n, m), sy(m, m), theta
1851 : INTEGER, INTENT(in) :: col, head
1852 : INTEGER, INTENT(out) :: info
1853 :
1854 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
1855 :
1856 : INTEGER :: col2, dbegin, dend, i, ipntr, is, is1, &
1857 : iy, jpntr, js, js1, jy, k, k1, m2, &
1858 : pbegin, pend, upcl
1859 : REAL(KIND=dp) :: ddot, temp1, temp2, temp3, temp4
1860 :
1861 : ! References:
1862 : ! [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
1863 : ! memory algorithm for bound constrained optimization'',
1864 : ! SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
1865 : !
1866 : ! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
1867 : ! limited memory FORTRAN code for solving bound constrained
1868 : ! optimization problems'', Tech. Report, NAM-11, EECS Department,
1869 : ! Northwestern University, 1994.
1870 : !
1871 : ! (Postscript files of these papers are available via anonymous
1872 : ! ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
1873 : !
1874 : ! * * *
1875 : ! Form the lower triangular part of
1876 : ! WN1 = [Y' ZZ'Y L_a'+R_z']
1877 : ! [L_a+R_z S'AA'S ]
1878 : ! where L_a is the strictly lower triangular part of S'AA'Y
1879 : ! R_z is the upper triangular part of S'ZZ'Y.
1880 :
1881 1436 : IF (updatd) THEN
1882 1436 : IF (iupdat > m) THEN
1883 : ! shift old part of WN1.
1884 6245 : DO jy = 1, m - 1
1885 4996 : js = m + jy
1886 4996 : CALL dcopy(m - jy, wn1(jy + 1, jy + 1), 1, wn1(jy, jy), 1)
1887 4996 : CALL dcopy(m - jy, wn1(js + 1, js + 1), 1, wn1(js, js), 1)
1888 6245 : CALL dcopy(m - 1, wn1(m + 2, jy + 1), 1, wn1(m + 1, jy), 1)
1889 : END DO
1890 : END IF
1891 :
1892 : ! put new rows in blocks (1,1), (2,1) and (2,2).
1893 1436 : pbegin = 1
1894 1436 : pend = nsub
1895 1436 : dbegin = nsub + 1
1896 1436 : dend = n
1897 1436 : iy = col
1898 1436 : is = m + col
1899 1436 : ipntr = head + col - 1
1900 1436 : IF (ipntr > m) ipntr = ipntr - m
1901 1436 : jpntr = head
1902 8656 : DO jy = 1, col
1903 7220 : js = m + jy
1904 7220 : temp1 = zero
1905 7220 : temp2 = zero
1906 7220 : temp3 = zero
1907 : ! compute element jy of row 'col' of Y'ZZ'Y
1908 3981683 : DO k = pbegin, pend
1909 3974463 : k1 = ind(k)
1910 3981683 : temp1 = temp1 + wy(k1, ipntr)*wy(k1, jpntr)
1911 : END DO
1912 : ! compute elements jy of row 'col' of L_a and S'AA'S
1913 7232 : DO k = dbegin, dend
1914 12 : k1 = ind(k)
1915 12 : temp2 = temp2 + ws(k1, ipntr)*ws(k1, jpntr)
1916 7232 : temp3 = temp3 + ws(k1, ipntr)*wy(k1, jpntr)
1917 : END DO
1918 7220 : wn1(iy, jy) = temp1
1919 7220 : wn1(is, js) = temp2
1920 7220 : wn1(is, jy) = temp3
1921 8656 : jpntr = MOD(jpntr, m) + 1
1922 : END DO
1923 :
1924 : ! put new column in block (2,1).
1925 1436 : jy = col
1926 1436 : jpntr = head + col - 1
1927 1436 : IF (jpntr > m) jpntr = jpntr - m
1928 : ipntr = head
1929 8656 : DO i = 1, col
1930 7220 : is = m + i
1931 7220 : temp3 = zero
1932 : ! compute element i of column 'col' of R_z
1933 3981683 : DO k = pbegin, pend
1934 3974463 : k1 = ind(k)
1935 3981683 : temp3 = temp3 + ws(k1, ipntr)*wy(k1, jpntr)
1936 : END DO
1937 7220 : ipntr = MOD(ipntr, m) + 1
1938 8656 : wn1(is, jy) = temp3
1939 : END DO
1940 1436 : upcl = col - 1
1941 : ELSE
1942 0 : upcl = col
1943 : END IF
1944 :
1945 : ! modify the old parts in blocks (1,1) and (2,2) due to changes
1946 : ! in the set of free variables.
1947 1436 : ipntr = head
1948 7220 : DO iy = 1, upcl
1949 5784 : is = m + iy
1950 5784 : jpntr = head
1951 22807 : DO jy = 1, iy
1952 17023 : js = m + jy
1953 17023 : temp1 = zero
1954 17023 : temp2 = zero
1955 17023 : temp3 = zero
1956 17023 : temp4 = zero
1957 17035 : DO k = 1, nenter
1958 12 : k1 = indx2(k)
1959 12 : temp1 = temp1 + wy(k1, ipntr)*wy(k1, jpntr)
1960 17035 : temp2 = temp2 + ws(k1, ipntr)*ws(k1, jpntr)
1961 : END DO
1962 17023 : DO k = ileave, n
1963 0 : k1 = indx2(k)
1964 0 : temp3 = temp3 + wy(k1, ipntr)*wy(k1, jpntr)
1965 17023 : temp4 = temp4 + ws(k1, ipntr)*ws(k1, jpntr)
1966 : END DO
1967 17023 : wn1(iy, jy) = wn1(iy, jy) + temp1 - temp3
1968 17023 : wn1(is, js) = wn1(is, js) - temp2 + temp4
1969 22807 : jpntr = MOD(jpntr, m) + 1
1970 : END DO
1971 7220 : ipntr = MOD(ipntr, m) + 1
1972 : END DO
1973 :
1974 : ! modify the old parts in block (2,1).
1975 1436 : ipntr = head
1976 7220 : DO is = m + 1, m + upcl
1977 : jpntr = head
1978 34046 : DO jy = 1, upcl
1979 28262 : temp1 = zero
1980 28262 : temp3 = zero
1981 28278 : DO k = 1, nenter
1982 16 : k1 = indx2(k)
1983 28278 : temp1 = temp1 + ws(k1, ipntr)*wy(k1, jpntr)
1984 : END DO
1985 28262 : DO k = ileave, n
1986 0 : k1 = indx2(k)
1987 28262 : temp3 = temp3 + ws(k1, ipntr)*wy(k1, jpntr)
1988 : END DO
1989 28262 : IF (is <= jy + m) THEN
1990 17023 : wn1(is, jy) = wn1(is, jy) + temp1 - temp3
1991 : ELSE
1992 11239 : wn1(is, jy) = wn1(is, jy) - temp1 + temp3
1993 : END IF
1994 34046 : jpntr = MOD(jpntr, m) + 1
1995 : END DO
1996 7220 : ipntr = MOD(ipntr, m) + 1
1997 : END DO
1998 :
1999 : ! Form the upper triangle of WN = [D+Y' ZZ'Y/theta -L_a'+R_z' ]
2000 : ! [-L_a +R_z S'AA'S*theta]
2001 :
2002 1436 : m2 = 2*m
2003 8656 : DO iy = 1, col
2004 7220 : is = col + iy
2005 7220 : is1 = m + iy
2006 31463 : DO jy = 1, iy
2007 24243 : js = col + jy
2008 24243 : js1 = m + jy
2009 24243 : wn(jy, iy) = wn1(iy, jy)/theta
2010 31463 : wn(js, is) = wn1(is1, js1)*theta
2011 : END DO
2012 24243 : DO jy = 1, iy - 1
2013 24243 : wn(jy, is) = -wn1(is1, jy)
2014 : END DO
2015 31463 : DO jy = iy, col
2016 31463 : wn(jy, is) = wn1(is1, jy)
2017 : END DO
2018 8656 : wn(iy, iy) = wn(iy, iy) + sy(iy, iy)
2019 : END DO
2020 :
2021 : ! Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')]
2022 : ! [(-L_a +R_z)L'^-1 S'AA'S*theta ]
2023 :
2024 : ! first Cholesky factor (1,1) block of wn to get LL'
2025 : ! with L' stored in the upper triangle of wn.
2026 1436 : CALL dpofa(wn, m2, col, info)
2027 1436 : IF (info /= 0) THEN
2028 0 : info = -1
2029 0 : RETURN
2030 : END IF
2031 : ! then form L^-1(-L_a'+R_z') in the (1,2) block.
2032 1436 : col2 = 2*col
2033 8656 : DO js = col + 1, col2
2034 8656 : CALL dtrsl(wn, m2, col, wn(1, js), 11, info)
2035 : END DO
2036 :
2037 : ! Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the
2038 : ! upper triangle of (2,2) block of wn.
2039 :
2040 8656 : DO is = col + 1, col2
2041 32899 : DO js = is, col2
2042 31463 : wn(is, js) = wn(is, js) + ddot(col, wn(1, is), 1, wn(1, js), 1)
2043 : END DO
2044 : END DO
2045 :
2046 : ! Cholesky factorization of (2,2) block of wn.
2047 :
2048 1436 : CALL dpofa(wn(col + 1, col + 1), m2, col, info)
2049 1436 : IF (info /= 0) THEN
2050 0 : info = -2
2051 0 : RETURN
2052 : END IF
2053 :
2054 : RETURN
2055 :
2056 : END SUBROUTINE formk
2057 :
2058 : ! **************************************************************************************************
2059 : !> \brief This subroutine forms the upper half of the pos. def. and symm.
2060 : !> T = theta*SS + L*D^(-1)*L', stores T in the upper triangle
2061 : !> of the array wt, and performs the Cholesky factorization of T
2062 : !> to produce J*J', with J' stored in the upper triangle of wt.
2063 : !> \param m ...
2064 : !> \param wt ...
2065 : !> \param sy ...
2066 : !> \param ss ...
2067 : !> \param col ...
2068 : !> \param theta ...
2069 : !> \param info ...
2070 : !> \author NEOS, November 1994. (Latest revision June 1996.)
2071 : !> Optimization Technology Center.
2072 : !> Argonne National Laboratory and Northwestern University.
2073 : !> Written by
2074 : !> Ciyou Zhu
2075 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
2076 : ! **************************************************************************************************
2077 1436 : SUBROUTINE formt(m, wt, sy, ss, col, theta, info)
2078 :
2079 : INTEGER :: m
2080 : REAL(KIND=dp) :: wt(m, m), sy(m, m), ss(m, m)
2081 : INTEGER :: col
2082 : REAL(KIND=dp) :: theta
2083 : INTEGER :: info
2084 :
2085 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
2086 :
2087 : INTEGER :: i, j, k, k1
2088 : REAL(KIND=dp) :: ddum
2089 :
2090 : ! Form the upper half of T = theta*SS + L*D^(-1)*L',
2091 : ! store T in the upper triangle of the array wt.
2092 :
2093 8656 : DO j = 1, col
2094 8656 : wt(1, j) = theta*ss(1, j)
2095 : END DO
2096 7220 : DO i = 2, col
2097 24243 : DO j = i, col
2098 17023 : k1 = MIN(i, j) - 1
2099 17023 : ddum = zero
2100 66619 : DO k = 1, k1
2101 66619 : ddum = ddum + sy(i, k)*sy(j, k)/sy(k, k)
2102 : END DO
2103 22807 : wt(i, j) = ddum + theta*ss(i, j)
2104 : END DO
2105 : END DO
2106 :
2107 : ! Cholesky factorize T to J*J' with
2108 : ! J' stored in the upper triangle of wt.
2109 :
2110 1436 : CALL dpofa(wt, m, col, info)
2111 1436 : IF (info /= 0) THEN
2112 0 : info = -3
2113 : END IF
2114 :
2115 1436 : RETURN
2116 :
2117 : END SUBROUTINE formt
2118 :
2119 : ! **************************************************************************************************
2120 : !> \brief This subroutine counts the entering and leaving variables when
2121 : !> iter > 0, and finds the index set of free and active variables
2122 : !> at the GCP.
2123 : !> \param n ...
2124 : !> \param nfree ...
2125 : !> \param index for i=1,...,nfree, index(i) are the indices of free variables
2126 : !> for i=nfree+1,...,n, index(i) are the indices of bound variables
2127 : !> On entry after the first iteration, index gives
2128 : !> the free variables at the previous iteration.
2129 : !> On exit it gives the free variables based on the determination
2130 : !> in cauchy using the array iwhere.
2131 : !> \param nenter ...
2132 : !> \param ileave ...
2133 : !> \param indx2 On exit with iter>0, indx2 indicates which variables
2134 : !> have changed status since the previous iteration.
2135 : !> For i= 1,...,nenter, indx2(i) have changed from bound to free.
2136 : !> For i= ileave+1,...,n, indx2(i) have changed from free to bound.
2137 : !> \param iwhere ...
2138 : !> \param wrk ...
2139 : !> \param updatd ...
2140 : !> \param constrained A variable indicating whether bounds are present
2141 : !> \param iprint ...
2142 : !> \param iter ...
2143 : !> \param iwunit User-specified write unit, if not set then WRITE statements
2144 : !> write to default_output_unit by default
2145 : !> \author NEOS, November 1994. (Latest revision June 1996.)
2146 : !> Optimization Technology Center.
2147 : !> Argonne National Laboratory and Northwestern University.
2148 : !> Written by
2149 : !> Ciyou Zhu
2150 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
2151 : ! **************************************************************************************************
2152 61 : SUBROUTINE freev(n, nfree, index, nenter, ileave, indx2, &
2153 61 : iwhere, wrk, updatd, constrained, iprint, iter, iwunit)
2154 :
2155 : INTEGER :: n, nfree
2156 : INTEGER, INTENT(inout) :: INDEX(n)
2157 : INTEGER :: nenter, ileave
2158 : INTEGER, INTENT(out) :: indx2(n)
2159 : INTEGER :: iwhere(n)
2160 : LOGICAL :: wrk, updatd, constrained
2161 : INTEGER :: iprint, iter
2162 : INTEGER, OPTIONAL :: iwunit
2163 :
2164 : INTEGER :: i, iact, k, wunit
2165 :
2166 61 : wunit = default_output_unit
2167 61 : IF (PRESENT(iwunit)) THEN
2168 61 : IF (iwunit > 0) wunit = iwunit
2169 : END IF
2170 :
2171 61 : nenter = 0
2172 61 : ileave = n + 1
2173 61 : IF (iter > 0 .AND. constrained) THEN
2174 : ! count the entering and leaving variables.
2175 27109 : DO i = 1, nfree
2176 27092 : k = INDEX(i)
2177 :
2178 27109 : IF (iwhere(k) > 0) THEN
2179 2 : ileave = ileave - 1
2180 2 : indx2(ileave) = k
2181 2 : IF (iprint >= 100) WRITE (wunit, 1030) k
2182 : END IF
2183 : END DO
2184 27 : DO i = 1 + nfree, n
2185 10 : k = INDEX(i)
2186 27 : IF (iwhere(k) <= 0) THEN
2187 4 : nenter = nenter + 1
2188 4 : indx2(nenter) = k
2189 4 : IF (iprint >= 100) WRITE (wunit, 2030) k
2190 : END IF
2191 : END DO
2192 17 : IF (iprint >= 99) WRITE (wunit, 3030) n + 1 - ileave, nenter
2193 : END IF
2194 61 : wrk = (ileave < n + 1) .OR. (nenter > 0) .OR. updatd
2195 :
2196 : ! Find the index set of free and active variables at the GCP.
2197 :
2198 61 : nfree = 0
2199 61 : iact = n + 1
2200 55738 : DO i = 1, n
2201 55738 : IF (iwhere(i) <= 0) THEN
2202 55667 : nfree = nfree + 1
2203 55667 : INDEX(nfree) = i
2204 : ELSE
2205 10 : iact = iact - 1
2206 10 : INDEX(iact) = i
2207 : END IF
2208 : END DO
2209 61 : IF (iprint >= 99) WRITE (wunit, 4030) nfree, iter + 1
2210 :
2211 : 1030 FORMAT(' L-BFGS| Variable ', i12, ' leaves the set of free variables')
2212 : 2030 FORMAT(' L-BFGS| Variable ', i12, ' enters the set of free variables')
2213 : 3030 FORMAT(' L-BFGS| ', i12, ' variables leave; ', i12, ' variables enter')
2214 : 4030 FORMAT(' L-BFGS| ', i12, ' variables are free at GCP ', i12)
2215 :
2216 61 : RETURN
2217 :
2218 : END SUBROUTINE freev
2219 :
2220 : ! **************************************************************************************************
2221 : !> \brief This subroutine sorts out the least element of t, and puts the
2222 : !> remaining elements of t in a heap.
2223 : !> \param n n is the dimension of the arrays t and iorder.
2224 : !> \param t On entry t stores the elements to be sorted,
2225 : !> On exit t(n) stores the least elements of t, and t(1) to t(n-1)
2226 : !> stores the remaining elements in the form of a heap.
2227 : !> \param iorder On entry iorder(i) is the index of t(i).
2228 : !> On exit iorder(i) is still the index of t(i), but iorder may be
2229 : !> permuted in accordance with t.
2230 : !> \param iheap iheap should be set as follows:
2231 : !> iheap .eq. 0 if t(1) to t(n) is not in the form of a heap,
2232 : !> iheap .ne. 0 if otherwise.
2233 : !> \author NEOS, November 1994. (Latest revision June 1996.)
2234 : !> Optimization Technology Center.
2235 : !> Argonne National Laboratory and Northwestern University.
2236 : !> Written by
2237 : !> Ciyou Zhu
2238 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
2239 : ! **************************************************************************************************
2240 10 : SUBROUTINE hpsolb(n, t, iorder, iheap)
2241 : INTEGER, INTENT(in) :: n
2242 : REAL(KIND=dp), INTENT(inout) :: t(n)
2243 : INTEGER, INTENT(inout) :: iorder(n)
2244 : INTEGER, INTENT(in) :: iheap
2245 :
2246 : INTEGER :: i, indxin, indxou, j, k
2247 : REAL(KIND=dp) :: ddum, out
2248 :
2249 : !
2250 : ! References:
2251 : ! Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT.
2252 : !
2253 : ! * * *
2254 :
2255 10 : IF (iheap == 0) THEN
2256 :
2257 : ! Rearrange the elements t(1) to t(n) to form a heap.
2258 :
2259 13516 : DO k = 2, n
2260 13511 : ddum = t(k)
2261 13511 : indxin = iorder(k)
2262 :
2263 : ! Add ddum to the heap.
2264 13511 : i = k
2265 32488 : DO WHILE (i > 1)
2266 32458 : j = i/2
2267 32488 : IF (ddum < t(j)) THEN
2268 18977 : t(i) = t(j)
2269 18977 : iorder(i) = iorder(j)
2270 18977 : i = j
2271 : ELSE
2272 : EXIT
2273 : END IF
2274 : END DO
2275 13511 : t(i) = ddum
2276 13516 : iorder(i) = indxin
2277 : END DO
2278 : END IF
2279 :
2280 : ! Assign to 'out' the value of t(1), the least member of the heap,
2281 : ! and rearrange the remaining members to form a heap as
2282 : ! elements 1 to n-1 of t.
2283 :
2284 10 : IF (n > 1) THEN
2285 10 : i = 1
2286 10 : out = t(1)
2287 10 : indxou = iorder(1)
2288 10 : ddum = t(n)
2289 10 : indxin = iorder(n)
2290 :
2291 : ! Restore the heap
2292 10 : j = 2*i
2293 85 : DO WHILE (j <= n - 1)
2294 77 : IF (t(j + 1) < t(j)) j = j + 1
2295 77 : IF (t(j) < ddum) THEN
2296 75 : t(i) = t(j)
2297 75 : iorder(i) = iorder(j)
2298 75 : i = j
2299 : ELSE
2300 : EXIT
2301 : END IF
2302 77 : j = 2*i
2303 : END DO
2304 10 : t(i) = ddum
2305 10 : iorder(i) = indxin
2306 :
2307 : ! Put the least member in t(n).
2308 :
2309 10 : t(n) = out
2310 10 : iorder(n) = indxou
2311 : END IF
2312 :
2313 10 : RETURN
2314 :
2315 : END SUBROUTINE hpsolb
2316 :
2317 : ! **************************************************************************************************
2318 : !> \brief This subroutine calls subroutine dcsrch from the Minpack2 library
2319 : !> to perform the line search. Subroutine dscrch is safeguarded so
2320 : !> that all trial points lie within the feasible region.
2321 : !> \param n ...
2322 : !> \param lower_bound the lower bound on x.
2323 : !> \param upper_bound the upper bound on x.
2324 : !> \param nbd ...
2325 : !> \param x ...
2326 : !> \param f ...
2327 : !> \param fold ...
2328 : !> \param gd ...
2329 : !> \param gdold ...
2330 : !> \param g ...
2331 : !> \param d ...
2332 : !> \param r ...
2333 : !> \param t ...
2334 : !> \param z ...
2335 : !> \param stp ...
2336 : !> \param dnorm ...
2337 : !> \param dtd ...
2338 : !> \param xstep ...
2339 : !> \param step_max ...
2340 : !> \param iter ...
2341 : !> \param ifun ...
2342 : !> \param iback ...
2343 : !> \param nfgv ...
2344 : !> \param info ...
2345 : !> \param task ...
2346 : !> \param boxed ...
2347 : !> \param constrained ...
2348 : !> \param csave ...
2349 : !> \param isave ...
2350 : !> \param dsave ...
2351 : !> \param iwunit User-specified write unit, if not set then WRITE statements
2352 : !> write to default_output_unit by default
2353 : !> \author NEOS, November 1994. (Latest revision June 1996.)
2354 : !> Optimization Technology Center.
2355 : !> Argonne National Laboratory and Northwestern University.
2356 : !> Written by
2357 : !> Ciyou Zhu
2358 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
2359 : ! **************************************************************************************************
2360 3135 : SUBROUTINE lnsrlb(n, lower_bound, upper_bound, nbd, x, f, fold, gd, gdold, g, d, r, t, &
2361 3135 : z, stp, dnorm, dtd, xstep, step_max, iter, ifun, &
2362 : iback, nfgv, info, task, boxed, constrained, csave, &
2363 : isave, dsave, iwunit)
2364 :
2365 : INTEGER, INTENT(in) :: n
2366 : REAL(KIND=dp), INTENT(in) :: lower_bound(n), upper_bound(n)
2367 : INTEGER :: nbd(n)
2368 : REAL(KIND=dp) :: x(n), f, fold, gd, gdold, g(n), d(n), &
2369 : r(n), t(n), z(n), stp, dnorm, dtd, &
2370 : xstep, step_max
2371 : INTEGER :: iter, ifun, iback, nfgv, info
2372 : CHARACTER(LEN=60) :: task
2373 : LOGICAL :: boxed, constrained
2374 : CHARACTER(LEN=60) :: csave
2375 : INTEGER :: isave(2)
2376 : REAL(KIND=dp) :: dsave(13)
2377 : INTEGER, OPTIONAL :: iwunit
2378 :
2379 : REAL(KIND=dp), PARAMETER :: big = 1.0E10_dp, ftol = 1.0E-3_dp, &
2380 : gtol = 0.9_dp, one = 1.0_dp, &
2381 : xtol = 0.1_dp, zero = 0.0_dp
2382 :
2383 : INTEGER :: i, wunit
2384 : REAL(KIND=dp) :: a1, a2, ddot
2385 :
2386 3135 : wunit = default_output_unit
2387 3135 : IF (PRESENT(iwunit)) THEN
2388 3135 : IF (iwunit > 0) wunit = iwunit
2389 : END IF
2390 :
2391 3135 : IF (.NOT. (task(1:5) == 'FG_LN')) THEN
2392 :
2393 1484 : dtd = ddot(n, d, 1, d, 1)
2394 1484 : dnorm = SQRT(dtd)
2395 :
2396 : ! Determine the maximum step length.
2397 :
2398 1484 : step_max = big
2399 1484 : IF (constrained) THEN
2400 20 : IF (iter == 0) THEN
2401 3 : step_max = one
2402 : ELSE
2403 27119 : DO i = 1, n
2404 27102 : a1 = d(i)
2405 27119 : IF (nbd(i) /= 0) THEN
2406 27102 : IF (a1 < zero .AND. nbd(i) <= 2) THEN
2407 8594 : a2 = lower_bound(i) - x(i)
2408 8594 : IF (a2 >= zero) THEN
2409 0 : step_max = zero
2410 8594 : ELSE IF (a1*step_max < a2) THEN
2411 13 : step_max = a2/a1
2412 : END IF
2413 18508 : ELSE IF (a1 > zero .AND. nbd(i) >= 2) THEN
2414 7872 : a2 = upper_bound(i) - x(i)
2415 7872 : IF (a2 <= zero) THEN
2416 0 : step_max = zero
2417 7872 : ELSE IF (a1*step_max > a2) THEN
2418 12 : step_max = a2/a1
2419 : END IF
2420 : END IF
2421 : END IF
2422 : END DO
2423 : END IF
2424 : END IF
2425 :
2426 1484 : IF (iter == 0 .AND. .NOT. boxed) THEN
2427 41 : stp = MIN(one/dnorm, step_max)
2428 : ELSE
2429 1443 : stp = one
2430 : END IF
2431 :
2432 1484 : CALL dcopy(n, x, 1, t, 1)
2433 1484 : CALL dcopy(n, g, 1, r, 1)
2434 1484 : fold = f
2435 1484 : ifun = 0
2436 1484 : iback = 0
2437 1484 : csave = 'START'
2438 : END IF
2439 3135 : gd = ddot(n, g, 1, d, 1)
2440 3135 : IF (ifun == 0) THEN
2441 1484 : gdold = gd
2442 1484 : IF (gd >= zero) THEN
2443 : ! the directional derivative >=0.
2444 : ! Line search is impossible.
2445 0 : WRITE (wunit, 1020) gd
2446 0 : info = -4
2447 0 : RETURN
2448 : END IF
2449 : END IF
2450 :
2451 3135 : CALL dcsrch(f, gd, stp, ftol, gtol, xtol, zero, step_max, csave, isave, dsave)
2452 :
2453 3135 : xstep = stp*dnorm
2454 3135 : IF (csave(1:4) /= 'CONV' .AND. csave(1:4) /= 'WARN') THEN
2455 1651 : task = 'FG_LNSRCH'
2456 1651 : ifun = ifun + 1
2457 1651 : nfgv = nfgv + 1
2458 1651 : iback = ifun - 1
2459 1651 : IF (stp == one) THEN
2460 1443 : CALL dcopy(n, z, 1, x, 1)
2461 : ELSE
2462 167926 : DO i = 1, n
2463 167926 : x(i) = stp*d(i) + t(i)
2464 : END DO
2465 : END IF
2466 : ELSE
2467 1484 : task = 'NEW_X'
2468 : END IF
2469 :
2470 : 1020 FORMAT(' L-BFGS| ascent direction in projection gd = ', d12.5)
2471 :
2472 : RETURN
2473 :
2474 : END SUBROUTINE lnsrlb
2475 :
2476 : ! **************************************************************************************************
2477 : !> \brief This subroutine updates matrices WS and WY, and forms the middle matrix in B.
2478 : !> \param n ...
2479 : !> \param m ...
2480 : !> \param ws ...
2481 : !> \param wy ...
2482 : !> \param sy ...
2483 : !> \param ss ...
2484 : !> \param d ...
2485 : !> \param r ...
2486 : !> \param itail ...
2487 : !> \param iupdat ...
2488 : !> \param col ...
2489 : !> \param head ...
2490 : !> \param theta ...
2491 : !> \param rr ...
2492 : !> \param dr ...
2493 : !> \param stp ...
2494 : !> \param dtd ...
2495 : !> \author NEOS, November 1994. (Latest revision June 1996.)
2496 : !> Optimization Technology Center.
2497 : !> Argonne National Laboratory and Northwestern University.
2498 : !> Written by
2499 : !> Ciyou Zhu
2500 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
2501 : ! **************************************************************************************************
2502 1436 : SUBROUTINE matupd(n, m, ws, wy, sy, ss, d, r, itail, &
2503 : iupdat, col, head, theta, rr, dr, stp, dtd)
2504 :
2505 : INTEGER :: n, m
2506 : REAL(KIND=dp) :: ws(n, m), wy(n, m), sy(m, m), ss(m, m), &
2507 : d(n), r(n)
2508 : INTEGER :: itail, iupdat, col, head
2509 : REAL(KIND=dp) :: theta, rr, dr, stp, dtd
2510 :
2511 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp
2512 :
2513 : INTEGER :: j, pointr
2514 : REAL(KIND=dp) :: ddot
2515 :
2516 : ! ************
2517 : ! Set pointers for matrices WS and WY.
2518 :
2519 1436 : IF (iupdat <= m) THEN
2520 187 : col = iupdat
2521 187 : itail = MOD(head + iupdat - 2, m) + 1
2522 : ELSE
2523 1249 : itail = MOD(itail, m) + 1
2524 1249 : head = MOD(head, m) + 1
2525 : END IF
2526 :
2527 : ! Update matrices WS and WY.
2528 :
2529 1436 : CALL dcopy(n, d, 1, ws(1, itail), 1)
2530 1436 : CALL dcopy(n, r, 1, wy(1, itail), 1)
2531 :
2532 : ! Set theta=yy/ys.
2533 :
2534 1436 : theta = rr/dr
2535 :
2536 : ! Form the middle matrix in B.
2537 :
2538 : ! update the upper triangle of SS,
2539 : ! and the lower triangle of SY:
2540 1436 : IF (iupdat > m) THEN
2541 : ! move old information
2542 6245 : DO j = 1, col - 1
2543 4996 : CALL dcopy(j, ss(2, j + 1), 1, ss(1, j), 1)
2544 6245 : CALL dcopy(col - j, sy(j + 1, j + 1), 1, sy(j, j), 1)
2545 : END DO
2546 : END IF
2547 : ! add new information: the last row of SY
2548 : ! and the last column of SS:
2549 1436 : pointr = head
2550 7220 : DO j = 1, col - 1
2551 5784 : sy(col, j) = ddot(n, d, 1, wy(1, pointr), 1)
2552 5784 : ss(j, col) = ddot(n, ws(1, pointr), 1, d, 1)
2553 7220 : pointr = MOD(pointr, m) + 1
2554 : END DO
2555 1436 : IF (stp == one) THEN
2556 1305 : ss(col, col) = dtd
2557 : ELSE
2558 131 : ss(col, col) = stp*stp*dtd
2559 : END IF
2560 1436 : sy(col, col) = dr
2561 :
2562 1436 : RETURN
2563 :
2564 : END SUBROUTINE matupd
2565 :
2566 : ! **************************************************************************************************
2567 : !> \brief This subroutine prints the input data, initial point, upper and
2568 : !> lower bounds of each variable, machine precision, as well as
2569 : !> the headings of the output.
2570 : !>
2571 : !> \param n ...
2572 : !> \param m ...
2573 : !> \param lower_bound the lower bound on x.
2574 : !> \param upper_bound the upper bound on x.
2575 : !> \param x ...
2576 : !> \param iprint ...
2577 : !> \param itfile ...
2578 : !> \param epsmch ...
2579 : !> \param iwunit User-specified write unit, if not set then WRITE statements
2580 : !> write to default_output_unit by default
2581 : !> \author NEOS, November 1994. (Latest revision June 1996.)
2582 : !> Optimization Technology Center.
2583 : !> Argonne National Laboratory and Northwestern University.
2584 : !> Written by
2585 : !> Ciyou Zhu
2586 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
2587 : ! **************************************************************************************************
2588 45 : SUBROUTINE prn1lb(n, m, lower_bound, upper_bound, x, iprint, itfile, epsmch, iwunit)
2589 :
2590 : INTEGER, INTENT(in) :: n, m
2591 : REAL(KIND=dp), INTENT(in) :: lower_bound(n), upper_bound(n), x(n)
2592 : INTEGER :: iprint, itfile
2593 : REAL(KIND=dp) :: epsmch
2594 : INTEGER, OPTIONAL :: iwunit
2595 :
2596 : INTEGER :: i, wunit
2597 :
2598 45 : wunit = default_output_unit
2599 45 : IF (PRESENT(iwunit)) THEN
2600 45 : IF (iwunit > 0) wunit = iwunit
2601 : END IF
2602 :
2603 45 : IF (iprint >= 0) THEN
2604 45 : WRITE (wunit, 7001) epsmch
2605 45 : WRITE (wunit, 7002) n, m
2606 45 : IF (iprint >= 1) THEN
2607 45 : WRITE (itfile, 2001) epsmch
2608 45 : WRITE (itfile, 7003) n, m
2609 45 : WRITE (itfile, 9001)
2610 45 : IF (iprint > 100) THEN
2611 0 : WRITE (wunit, 1004) ' L-BFGS| L =', (lower_bound(i), i=1, n)
2612 0 : WRITE (wunit, 1004) ' L-BFGS| X0 =', (x(i), i=1, n)
2613 0 : WRITE (wunit, 1004) ' L-BFGS| U =', (upper_bound(i), i=1, n)
2614 : END IF
2615 : END IF
2616 : END IF
2617 :
2618 : 1004 FORMAT(/, a13, 1p, /, (4x, 1p, 6(1x, d11.4)))
2619 : 2001 FORMAT('RUNNING THE L-BFGS-B CODE', /, /, &
2620 : 'it = iteration number', /, &
2621 : 'nf = number of function evaluations', /, &
2622 : 'nseg = number of segments explored during the Cauchy search', /, &
2623 : 'nact = number of active bounds at the generalized Cauchy point' &
2624 : , /, &
2625 : 'sub = manner in which the subspace minimization terminated:' &
2626 : , /, ' con = converged, bnd = a bound was reached', /, &
2627 : 'itls = number of iterations performed in the line search', /, &
2628 : 'stepl = step length used', /, &
2629 : 'tstep = norm of the displacement (total step)', /, &
2630 : 'projg = norm of the projected gradient', /, &
2631 : 'f = function value', /, /, &
2632 : ' * * *', /, /, &
2633 : 'Machine precision =', 1p, d10.3)
2634 : 7001 FORMAT(/, ' L-BFGS| RUNNING THE L-BFGS-B CODE', /, &
2635 : ' L-BFGS| Machine precision =', 1p, d10.3)
2636 : 7002 FORMAT(/, ' L-BFGS| N = ', i12, ' M = ', i12)
2637 : 7003 FORMAT(' N = ', i12, ' M = ', i12)
2638 : 9001 FORMAT(/, 3x, 'it', 3x, 'nf', 2x, 'nseg', 2x, 'nact', 2x, 'sub', 2x, 'itls', &
2639 : 2x, 'stepl', 4x, 'tstep', 5x, 'projg', 8x, 'f')
2640 :
2641 45 : RETURN
2642 :
2643 : END SUBROUTINE prn1lb
2644 :
2645 : ! **************************************************************************************************
2646 : !> \brief This subroutine prints out new information after a successful line search.
2647 : !> \param n ...
2648 : !> \param x ...
2649 : !> \param f ...
2650 : !> \param g ...
2651 : !> \param iprint ...
2652 : !> \param itfile ...
2653 : !> \param iter ...
2654 : !> \param nfgv ...
2655 : !> \param nact ...
2656 : !> \param g_inf_norm ...
2657 : !> \param nseg ...
2658 : !> \param word ...
2659 : !> \param iword ...
2660 : !> \param iback ...
2661 : !> \param stp ...
2662 : !> \param xstep ...
2663 : !> \param iwunit User-specified write unit, if not set then WRITE statements
2664 : !> write to default_output_unit by default
2665 : !> \author NEOS, November 1994. (Latest revision June 1996.)
2666 : !> Optimization Technology Center.
2667 : !> Argonne National Laboratory and Northwestern University.
2668 : !> Written by
2669 : !> Ciyou Zhu
2670 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
2671 : ! **************************************************************************************************
2672 1484 : SUBROUTINE prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact, &
2673 : g_inf_norm, nseg, word, iword, iback, stp, xstep, iwunit)
2674 :
2675 : INTEGER, INTENT(in) :: n
2676 : REAL(KIND=dp), INTENT(in) :: x(n), f, g(n)
2677 : INTEGER, INTENT(in) :: iprint, itfile, iter, nfgv, nact
2678 : REAL(KIND=dp), INTENT(in) :: g_inf_norm
2679 : INTEGER, INTENT(in) :: nseg
2680 : CHARACTER(LEN=3) :: word
2681 : INTEGER :: iword, iback
2682 : REAL(KIND=dp) :: stp, xstep
2683 : INTEGER, OPTIONAL :: iwunit
2684 :
2685 : INTEGER :: i, imod, wunit
2686 :
2687 1484 : wunit = default_output_unit
2688 1484 : IF (PRESENT(iwunit)) THEN
2689 1484 : IF (iwunit > 0) wunit = iwunit
2690 : END IF
2691 :
2692 : ! 'word' records the status of subspace solutions.
2693 :
2694 1484 : IF (iword == 0) THEN
2695 : ! the subspace minimization converged.
2696 1436 : word = 'con'
2697 48 : ELSE IF (iword == 1) THEN
2698 : ! the subspace minimization stopped at a bound.
2699 0 : word = 'bnd'
2700 48 : ELSE IF (iword == 5) THEN
2701 : ! the truncated Newton step has been used.
2702 0 : word = 'TNT'
2703 : ELSE
2704 48 : word = '---'
2705 : END IF
2706 1484 : IF (iprint >= 99) THEN
2707 0 : WRITE (wunit, 2002) iback, xstep
2708 0 : WRITE (wunit, 2001) iter, f, g_inf_norm
2709 0 : IF (iprint > 100) THEN
2710 0 : WRITE (wunit, 1004) ' L-BFGS| X =', (x(i), i=1, n)
2711 0 : WRITE (wunit, 1004) ' L-BFGS| G =', (g(i), i=1, n)
2712 : END IF
2713 1484 : ELSE IF (iprint > 0) THEN
2714 1484 : imod = MOD(iter, iprint)
2715 1484 : IF (imod == 0) WRITE (wunit, 2001) iter, f, g_inf_norm
2716 : END IF
2717 1484 : IF (iprint >= 1) WRITE (itfile, 3001) &
2718 1484 : iter, nfgv, nseg, nact, word, iback, stp, xstep, g_inf_norm, f
2719 :
2720 : 1004 FORMAT(/, a12, 1p, /, (4x, 1p, 6(1x, d11.4)))
2721 : 2001 FORMAT &
2722 : (/, ' L-BFGS| At iterate', i5, 4x, 'f= ', 1p, d12.5, 4x, '|proj g|= ', 1p, d12.5)
2723 : 2002 FORMAT(/, ' L-BFGS| LINE SEARCH ', i12, ' times; norm of step = ', 1p, d24.15)
2724 : 3001 FORMAT(2(1x, i4), 2(1x, i5), 2x, a3, 1x, i4, 1p, 2(2x, d7.1), 1p, 2(1x, d10.3))
2725 :
2726 1484 : RETURN
2727 :
2728 : END SUBROUTINE prn2lb
2729 :
2730 : ! **************************************************************************************************
2731 : !> \brief This subroutine prints out information when either a built-in
2732 : !> convergence test is satisfied or when an error message is
2733 : !> generated.
2734 : !> \param n ...
2735 : !> \param x ...
2736 : !> \param f ...
2737 : !> \param task ...
2738 : !> \param iprint ...
2739 : !> \param info ...
2740 : !> \param itfile ...
2741 : !> \param iter ...
2742 : !> \param nfgv ...
2743 : !> \param nintol ...
2744 : !> \param nskip ...
2745 : !> \param nact ...
2746 : !> \param g_inf_norm ...
2747 : !> \param time ...
2748 : !> \param nseg ...
2749 : !> \param word ...
2750 : !> \param iback ...
2751 : !> \param stp ...
2752 : !> \param xstep ...
2753 : !> \param k ...
2754 : !> \param cachyt ...
2755 : !> \param sbtime ...
2756 : !> \param lnscht ...
2757 : !> \param iwunit User-specified write unit, if not set then WRITE statements
2758 : !> write to default_output_unit by default
2759 : !> \author NEOS, November 1994. (Latest revision June 1996.)
2760 : !> Optimization Technology Center.
2761 : !> Argonne National Laboratory and Northwestern University.
2762 : !> Written by
2763 : !> Ciyou Zhu
2764 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
2765 : ! **************************************************************************************************
2766 1 : SUBROUTINE prn3lb(n, x, f, task, iprint, info, itfile, &
2767 : iter, nfgv, nintol, nskip, nact, g_inf_norm, &
2768 : time, nseg, word, iback, stp, xstep, k, &
2769 : cachyt, sbtime, lnscht, iwunit)
2770 :
2771 : INTEGER, INTENT(in) :: n
2772 : REAL(KIND=dp), INTENT(in) :: x(n), f
2773 : CHARACTER(LEN=60), INTENT(in) :: task
2774 : INTEGER, INTENT(in) :: iprint, info, itfile, iter, nfgv, &
2775 : nintol, nskip, nact
2776 : REAL(KIND=dp), INTENT(in) :: g_inf_norm, time
2777 : INTEGER, INTENT(in) :: nseg
2778 : CHARACTER(LEN=3) :: word
2779 : INTEGER :: iback
2780 : REAL(KIND=dp) :: stp, xstep
2781 : INTEGER :: k
2782 : REAL(KIND=dp) :: cachyt, sbtime, lnscht
2783 : INTEGER, OPTIONAL :: iwunit
2784 :
2785 : INTEGER :: i, wunit
2786 :
2787 1 : wunit = default_output_unit
2788 1 : IF (PRESENT(iwunit)) THEN
2789 1 : IF (iwunit > 0) wunit = iwunit
2790 : END IF
2791 :
2792 1 : IF (iprint >= 0 .AND. .NOT. (task(1:5) == 'ERROR')) THEN
2793 1 : WRITE (wunit, 3003)
2794 1 : WRITE (wunit, 3004)
2795 1 : WRITE (wunit, 3005) n, iter, nfgv, nintol, nskip, nact, g_inf_norm, f
2796 1 : IF (iprint >= 100) THEN
2797 0 : WRITE (wunit, 1004) ' L-BFGS| X =', (x(i), i=1, n)
2798 : END IF
2799 1 : IF (iprint >= 1) WRITE (wunit, 3006) f
2800 : END IF
2801 1 : IF (iprint >= 0) THEN
2802 :
2803 1 : WRITE (wunit, 3001)
2804 1 : WRITE (wunit, 3009) task
2805 1 : IF (info /= 0) THEN
2806 0 : IF (info == -1) WRITE (wunit, 9011)
2807 0 : IF (info == -2) WRITE (wunit, 9012)
2808 0 : IF (info == -3) WRITE (wunit, 9013)
2809 0 : IF (info == -4) WRITE (wunit, 9014)
2810 0 : IF (info == -5) WRITE (wunit, 9015)
2811 0 : IF (info == -6) WRITE (wunit, 9016) k
2812 0 : IF (info == -7) WRITE (wunit, 9017) k, k
2813 0 : IF (info == -8) WRITE (wunit, 9018)
2814 0 : IF (info == -9) WRITE (wunit, 9019)
2815 : END IF
2816 1 : IF (iprint >= 1) WRITE (wunit, 3007) cachyt, sbtime, lnscht
2817 1 : WRITE (wunit, 3008) time
2818 1 : WRITE (wunit, 3001)
2819 :
2820 1 : IF (iprint >= 1) THEN
2821 1 : IF (info == -4 .OR. info == -9) THEN
2822 : WRITE (itfile, 3002) &
2823 0 : iter, nfgv, nseg, nact, word, iback, stp, xstep
2824 : END IF
2825 1 : WRITE (itfile, 4009) task
2826 1 : IF (info /= 0) THEN
2827 0 : IF (info == -1) WRITE (itfile, 9011)
2828 0 : IF (info == -2) WRITE (itfile, 9012)
2829 0 : IF (info == -3) WRITE (itfile, 9013)
2830 0 : IF (info == -4) WRITE (itfile, 9014)
2831 0 : IF (info == -5) WRITE (itfile, 9015)
2832 0 : IF (info == -8) WRITE (itfile, 9018)
2833 0 : IF (info == -9) WRITE (itfile, 9019)
2834 : END IF
2835 1 : WRITE (itfile, 3008) time
2836 : END IF
2837 : END IF
2838 :
2839 : 1004 FORMAT(/, a12, 1p, /, (4x, 1p, 6(1x, d11.4)))
2840 : 3001 FORMAT(/, ' L-BFGS| ---------------- Information ----------------')
2841 : 3002 FORMAT(2(1x, i4), 2(1x, i5), 2x, a3, 1x, i4, 1p, 2(2x, d7.1), 6x, '-', 10x, '-')
2842 : 3003 FORMAT(/, &
2843 : ' L-BFGS| * * *', /, /, &
2844 : ' L-BFGS| Tit = total number of iterations', /, &
2845 : ' L-BFGS| Tnf = total number of function evaluations', /, &
2846 : ' L-BFGS| Tnint = total number of segments explored during', &
2847 : ' L-BFGS| Cauchy searches', /, &
2848 : ' L-BFGS| Skip = number of BFGS updates skipped', /, &
2849 : ' L-BFGS| Nact = number of active bounds at final generalized', &
2850 : ' L-BFGS| Cauchy point', /, &
2851 : ' L-BFGS| Projg = norm of the final projected gradient', /, &
2852 : ' L-BFGS| F = final function value', /, /, &
2853 : ' L-BFGS| * * *')
2854 : 3004 FORMAT(/, ' L-BFGS| ', 3x, 'N', 4x, 'Tit', 5x, 'Tnf', 2x, 'Tnint', 2x, &
2855 : 'Skip', 2x, 'Nact', 5x, 'Projg', 8x, 'F')
2856 : 3005 FORMAT(' L-BFGS| ', i5, 2(1x, i6), (1x, i6), (2x, i4), (1x, i5), 1p, 2(2x, d10.3))
2857 : 3006 FORMAT(' L-BFGS| F =', d12.5)
2858 : 3007 FORMAT(/, &
2859 : ' L-BFGS| Cauchy time', 1p, e10.3, ' seconds.', / &
2860 : ' L-BFGS| Subspace minimization time', 1p, e10.3, ' seconds.', / &
2861 : ' L-BFGS| Line search time', 1p, e10.3, ' seconds.')
2862 : 3008 FORMAT(/, ' Total User time', 1p, e10.3, ' seconds.',/)
2863 : 3009 FORMAT(/, ' L-BFGS| ', a60)
2864 : 4009 FORMAT(/, a60)
2865 : 9011 FORMAT(/, &
2866 : ' Matrix in 1st Cholesky factorization in formk is not Pos. Def.')
2867 : 9012 FORMAT(/, &
2868 : ' Matrix in 2st Cholesky factorization in formk is not Pos. Def.')
2869 : 9013 FORMAT(/, &
2870 : ' Matrix in the Cholesky factorization in formt is not Pos. Def.')
2871 : 9014 FORMAT(/, &
2872 : ' Derivative >= 0, backtracking line search impossible.', /, &
2873 : ' Previous x, f and g restored.', /, &
2874 : ' Possible causes: 1 error in function or gradient evaluation;', /, &
2875 : ' 2 rounding errors dominate computation.')
2876 : 9015 FORMAT(/, &
2877 : ' Warning: more than 10 function and gradient', /, &
2878 : ' evaluations in the last line search. Termination', /, &
2879 : ' may possibly be caused by a bad search direction.')
2880 : 9016 FORMAT(' Input nbd(', i12, ') is invalid.')
2881 : 9017 FORMAT(' l(', i12, ') > u(', i12, '). No feasible solution.')
2882 : 9018 FORMAT(/, ' The triangular system is singular.')
2883 : 9019 FORMAT(/, &
2884 : ' Line search cannot locate an adequate point after 20 function', /, &
2885 : ' and gradient evaluations. Previous x, f and g restored.', /, &
2886 : ' Possible causes: 1 error in function or gradient evaluation;', /, &
2887 : ' 2 rounding error dominate computation.')
2888 :
2889 1 : RETURN
2890 :
2891 : END SUBROUTINE prn3lb
2892 :
2893 : ! **************************************************************************************************
2894 : !> \brief This subroutine computes the infinity norm of the projected gradient.
2895 : !> \param n ...
2896 : !> \param lower_bound the lower bound on x.
2897 : !> \param upper_bound the upper bound on x.
2898 : !> \param nbd ...
2899 : !> \param x ...
2900 : !> \param g ...
2901 : !> \param g_inf_norm ...
2902 : !> \author NEOS, November 1994. (Latest revision June 1996.)
2903 : !> Optimization Technology Center.
2904 : !> Argonne National Laboratory and Northwestern University.
2905 : !> Written by
2906 : !> Ciyou Zhu
2907 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
2908 : ! **************************************************************************************************
2909 1528 : SUBROUTINE projgr(n, lower_bound, upper_bound, nbd, x, g, g_inf_norm)
2910 :
2911 : INTEGER, INTENT(in) :: n
2912 : REAL(KIND=dp), INTENT(in) :: lower_bound(n), upper_bound(n)
2913 : INTEGER, INTENT(in) :: nbd(n)
2914 : REAL(KIND=dp), INTENT(in) :: x(n), g(n)
2915 : REAL(KIND=dp) :: g_inf_norm
2916 :
2917 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
2918 :
2919 : INTEGER :: i
2920 : REAL(KIND=dp) :: gi
2921 :
2922 1528 : g_inf_norm = zero
2923 894763 : DO i = 1, n
2924 893235 : gi = g(i)
2925 893235 : IF (nbd(i) /= 0) THEN
2926 45138 : IF (gi < zero) THEN
2927 22652 : IF (nbd(i) >= 2) gi = MAX((x(i) - upper_bound(i)), gi)
2928 : ELSE
2929 22486 : IF (nbd(i) <= 2) gi = MIN((x(i) - lower_bound(i)), gi)
2930 : END IF
2931 : END IF
2932 894763 : g_inf_norm = MAX(g_inf_norm, ABS(gi))
2933 : END DO
2934 :
2935 1528 : RETURN
2936 :
2937 : END SUBROUTINE projgr
2938 :
2939 : ! **************************************************************************************************
2940 : !> \brief This routine contains the major changes in the updated version.
2941 : !> The changes are described in the accompanying paper
2942 : !>
2943 : !> Jose Luis Morales, Jorge Nocedal
2944 : !> "Remark On Algorithm 788: L-BFGS-B: Fortran Subroutines for Large
2945 : !> Bound Constrained Optimization". Decemmber 27, 2010.
2946 : !>
2947 : !> J.L. Morales Departamento de Matematicas,
2948 : !> Instituto Tecnologico Autonomo de Mexico
2949 : !> Mexico D.F.
2950 : !>
2951 : !> J, Nocedal Department of Electrical Engineering and
2952 : !> Computer Science.
2953 : !> Northwestern University. Evanston, IL. USA
2954 : !>
2955 : !> January 17, 2011
2956 : !>
2957 : !> *****************************************************************
2958 : !>
2959 : !> Given xcp, l, u, r, an index set that specifies
2960 : !> the active set at xcp, and an l-BFGS matrix B
2961 : !> (in terms of WY, WS, SY, WT, head, col, and theta),
2962 : !> this subroutine computes an approximate solution
2963 : !> of the subspace problem
2964 : !>
2965 : !> (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)
2966 : !>
2967 : !> subject to l<=x<=u
2968 : !> x_i=xcp_i for all i in A(xcp)
2969 : !>
2970 : !> along the subspace unconstrained Newton direction
2971 : !>
2972 : !> d = -(Z'BZ)^(-1) r.
2973 : !>
2974 : !> The formula for the Newton direction, given the L-BFGS matrix
2975 : !> and the Sherman-Morrison formula, is
2976 : !>
2977 : !> d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.
2978 : !>
2979 : !> where
2980 : !> K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
2981 : !> [L_a -R_z theta*S'AA'S ]
2982 : !>
2983 : !> Note that this procedure for computing d differs
2984 : !> from that described in [1]. One can show that the matrix K is
2985 : !> equal to the matrix M^[-1]N in that paper.
2986 : !> \param n n is the dimension of the problem.
2987 : !> \param m m is the maximum number of variable metric corrections
2988 : !> used to define the limited memory matrix.
2989 : !> \param nsub nsub is the number of free variables.
2990 : !> \param ind ind specifies the coordinate indices of free variables.
2991 : !> \param lower_bound the lower bound on x.
2992 : !> \param upper_bound the upper bound on x.
2993 : !> \param nbd nbd represents the type of bounds imposed on the
2994 : !> variables, and must be specified as follows:
2995 : !> nbd(i)=0 if x(i) is unbounded,
2996 : !> 1 if x(i) has only a lower bound,
2997 : !> 2 if x(i) has both lower and upper bounds, and
2998 : !> 3 if x(i) has only an upper bound.
2999 : !> \param x x is a double precision array of dimension n.
3000 : !> On entry x specifies the Cauchy point xcp.
3001 : !> On exit x(i) is the minimizer of Q over the subspace of free variables.
3002 : !> \param d On entry d is the reduced gradient of Q at xcp.
3003 : !> On exit d is the Newton direction of Q.
3004 : !> \param xp xp is a double precision array of dimension n.
3005 : !> used to safeguard the projected Newton direction
3006 : !> \param ws ws and wy are double precision arrays;
3007 : !> On entry they store the information defining the limited memory BFGS matrix:
3008 : !> ws(n,m) stores S, a set of s-vectors;
3009 : !> \param wy wy(n,m) stores Y, a set of y-vectors;
3010 : !> \param theta theta is the scaling factor specifying B_0 = theta I;
3011 : !> \param xx xx holds the current iterate
3012 : !> \param gg gg holds the gradient at the current iterate
3013 : !> \param col is the number of variable metric corrections stored;
3014 : !> \param head head is the location of the 1st s- (or y-) vector in S (or Y).
3015 : !> \param iword iword specifies the status of the subspace solution.
3016 : !> iword = 0 if the solution is in the box,
3017 : !> 1 if some bound is encountered.
3018 : !> \param wv wv is a working array
3019 : !> \param wn the upper triangle of wn stores the LEL^T factorization
3020 : !> of the indefinite matrix
3021 : !>
3022 : !> K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
3023 : !> [L_a -R_z theta*S'AA'S ]
3024 : !> where E = [-I 0]
3025 : !> [ 0 I]
3026 : !> \param iprint iprint is an INTEGER variable that must be set by the user.
3027 : !> It controls the frequency and type of output generated:
3028 : !> iprint<0 no output is generated;
3029 : !> iprint=0 print only one line at the last iteration;
3030 : !> 0<iprint<99 print also f and |proj g| every iprint iterations;
3031 : !> iprint=99 print details of every iteration except n-vectors;
3032 : !> iprint=100 print also the changes of active set and final x;
3033 : !> iprint>100 print details of every iteration including x and g;
3034 : !> When iprint > 0, the file iterate.dat will be created to summarize the iteration.
3035 : !> \param info info = 0 for normal return,
3036 : !> = nonzero for abnormal return when the matrix K is ill-conditioned.
3037 : !> \param iwunit User-specified write unit, if not set then WRITE statements
3038 : !> write to default_output_unit by default
3039 : !> \author NEOS, November 1994. (Latest revision June 1996.)
3040 : !> Optimization Technology Center.
3041 : !> Argonne National Laboratory and Northwestern University.
3042 : !> Written by
3043 : !> Ciyou Zhu
3044 : !> in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
3045 : ! **************************************************************************************************
3046 1436 : SUBROUTINE subsm(n, m, nsub, ind, lower_bound, upper_bound, nbd, x, d, xp, ws, wy, &
3047 1436 : theta, xx, gg, &
3048 1436 : col, head, iword, wv, wn, iprint, info, iwunit)
3049 : INTEGER, INTENT(in) :: n, m, nsub, ind(nsub)
3050 : REAL(KIND=dp), INTENT(in) :: lower_bound(n), upper_bound(n)
3051 : INTEGER, INTENT(in) :: nbd(n)
3052 : REAL(KIND=dp), INTENT(inout) :: x(n), d(n)
3053 : REAL(KIND=dp) :: xp(n)
3054 : REAL(KIND=dp), INTENT(in) :: ws(n, m), wy(n, m), theta, xx(n), gg(n)
3055 : INTEGER, INTENT(in) :: col, head
3056 : INTEGER, INTENT(out) :: iword
3057 : REAL(KIND=dp) :: wv(2*m)
3058 : REAL(KIND=dp), INTENT(in) :: wn(2*m, 2*m)
3059 : INTEGER :: iprint
3060 : INTEGER, INTENT(out) :: info
3061 : INTEGER, OPTIONAL :: iwunit
3062 :
3063 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
3064 :
3065 : INTEGER :: col2, i, ibd, j, js, jy, k, m2, pointr, &
3066 : wunit
3067 : REAL(KIND=dp) :: alpha, dd_p, dk, temp1, temp2, xk
3068 :
3069 : ! References:
3070 : !
3071 : ! [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
3072 : ! memory algorithm for bound constrained optimization'',
3073 : ! SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
3074 : !
3075 : !
3076 : !
3077 : ! * * *
3078 : !
3079 :
3080 1436 : wunit = default_output_unit
3081 1436 : IF (PRESENT(iwunit)) THEN
3082 1436 : IF (iwunit > 0) wunit = iwunit
3083 : END IF
3084 :
3085 1436 : IF (nsub <= 0) RETURN
3086 1436 : IF (iprint >= 99) WRITE (wunit, 4001)
3087 :
3088 : ! Compute wv = W'Zd.
3089 :
3090 1436 : pointr = head
3091 8656 : DO i = 1, col
3092 : temp1 = zero
3093 : temp2 = zero
3094 3981683 : DO j = 1, nsub
3095 3974463 : k = ind(j)
3096 3974463 : temp1 = temp1 + wy(k, pointr)*d(j)
3097 3981683 : temp2 = temp2 + ws(k, pointr)*d(j)
3098 : END DO
3099 7220 : wv(i) = temp1
3100 7220 : wv(col + i) = theta*temp2
3101 8656 : pointr = MOD(pointr, m) + 1
3102 : END DO
3103 :
3104 : ! Compute wv:=K^(-1)wv.
3105 :
3106 1436 : m2 = 2*m
3107 1436 : col2 = 2*col
3108 1436 : CALL dtrsl(wn, m2, col2, wv, 11, info)
3109 1436 : IF (info /= 0) RETURN
3110 8656 : DO i = 1, col
3111 8656 : wv(i) = -wv(i)
3112 : END DO
3113 1436 : CALL dtrsl(wn, m2, col2, wv, 01, info)
3114 1436 : IF (info /= 0) RETURN
3115 :
3116 : ! Compute d = (1/theta)d + (1/theta**2)Z'W wv.
3117 :
3118 : pointr = head
3119 8656 : DO jy = 1, col
3120 7220 : js = col + jy
3121 3981683 : DO i = 1, nsub
3122 3974463 : k = ind(i)
3123 : d(i) = d(i) + wy(k, pointr)*wv(jy)/theta &
3124 3981683 : & + ws(k, pointr)*wv(js)
3125 : END DO
3126 8656 : pointr = MOD(pointr, m) + 1
3127 : END DO
3128 :
3129 1436 : CALL dscal(nsub, one/theta, d, 1)
3130 : !
3131 : !-----------------------------------------------------------------
3132 : ! Let us try the projection, d is the Newton direction
3133 :
3134 1436 : iword = 0
3135 :
3136 1436 : CALL dcopy(n, x, 1, xp, 1)
3137 : !
3138 837489 : DO i = 1, nsub
3139 836053 : k = ind(i)
3140 836053 : dk = d(i)
3141 836053 : xk = x(k)
3142 837489 : IF (nbd(k) /= 0) THEN
3143 : !
3144 : ! lower bounds only
3145 27070 : IF (nbd(k) == 1) THEN
3146 0 : x(k) = MAX(lower_bound(k), xk + dk)
3147 0 : IF (x(k) == lower_bound(k)) iword = 1
3148 : ELSE
3149 : !
3150 : ! upper and lower bounds
3151 27070 : IF (nbd(k) == 2) THEN
3152 27070 : xk = MAX(lower_bound(k), xk + dk)
3153 27070 : x(k) = MIN(upper_bound(k), xk)
3154 27070 : IF (x(k) == lower_bound(k) .OR. x(k) == upper_bound(k)) iword = 1
3155 : ELSE
3156 : !
3157 : ! upper bounds only
3158 0 : IF (nbd(k) == 3) THEN
3159 0 : x(k) = MIN(upper_bound(k), xk + dk)
3160 0 : IF (x(k) == upper_bound(k)) iword = 1
3161 : END IF
3162 : END IF
3163 : END IF
3164 : !
3165 : ! free variables
3166 : ELSE
3167 808983 : x(k) = xk + dk
3168 : END IF
3169 : END DO
3170 : !
3171 1436 : IF (.NOT. (iword == 0)) THEN
3172 : !
3173 : ! check sign of the directional derivative
3174 : !
3175 : dd_p = zero
3176 0 : DO i = 1, n
3177 0 : dd_p = dd_p + (x(i) - xx(i))*gg(i)
3178 : END DO
3179 0 : IF (dd_p > zero) THEN
3180 0 : CALL dcopy(n, xp, 1, x, 1)
3181 0 : IF (iprint > 0) WRITE (wunit, 4002)
3182 0 : IF (iprint > 0) WRITE (wunit, 4003)
3183 0 : alpha = one
3184 0 : temp1 = alpha
3185 0 : ibd = 0
3186 0 : DO i = 1, nsub
3187 0 : k = ind(i)
3188 0 : dk = d(i)
3189 0 : IF (nbd(k) /= 0) THEN
3190 0 : IF (dk < zero .AND. nbd(k) <= 2) THEN
3191 0 : temp2 = lower_bound(k) - x(k)
3192 0 : IF (temp2 >= zero) THEN
3193 : temp1 = zero
3194 0 : ELSE IF (dk*alpha < temp2) THEN
3195 0 : temp1 = temp2/dk
3196 : END IF
3197 0 : ELSE IF (dk > zero .AND. nbd(k) >= 2) THEN
3198 0 : temp2 = upper_bound(k) - x(k)
3199 0 : IF (temp2 <= zero) THEN
3200 : temp1 = zero
3201 0 : ELSE IF (dk*alpha > temp2) THEN
3202 0 : temp1 = temp2/dk
3203 : END IF
3204 : END IF
3205 0 : IF (temp1 < alpha) THEN
3206 0 : alpha = temp1
3207 0 : ibd = i
3208 : END IF
3209 : END IF
3210 : END DO
3211 :
3212 0 : IF (alpha < one) THEN
3213 0 : dk = d(ibd)
3214 0 : k = ind(ibd)
3215 0 : IF (dk > zero) THEN
3216 0 : x(k) = upper_bound(k)
3217 0 : d(ibd) = zero
3218 0 : ELSE IF (dk < zero) THEN
3219 0 : x(k) = lower_bound(k)
3220 0 : d(ibd) = zero
3221 : END IF
3222 : END IF
3223 0 : DO i = 1, nsub
3224 0 : k = ind(i)
3225 0 : x(k) = x(k) + alpha*d(i)
3226 : END DO
3227 : END IF
3228 : END IF
3229 :
3230 1436 : IF (iprint >= 99) WRITE (wunit, 4004)
3231 :
3232 : 4001 FORMAT(/, ' L-BFGS| ---------------- enter SUBSM ----------------',/)
3233 : 4002 FORMAT(' L-BFGS| Positive dir derivative in projection ')
3234 : 4003 FORMAT(' L-BFGS| Using the backtracking step ')
3235 : 4004 FORMAT(/, ' L-BFGS| ---------------- exit SUBSM -----------------',/)
3236 :
3237 : RETURN
3238 :
3239 : END SUBROUTINE subsm
3240 :
3241 : ! **************************************************************************************************
3242 : !> \brief This subroutine finds a step that satisfies a sufficient
3243 : !> decrease condition and a curvature condition.
3244 : !>
3245 : !> Each call of the subroutine updates an interval with
3246 : !> endpoints stx and sty. The interval is initially chosen
3247 : !> so that it contains a minimizer of the modified function
3248 : !>
3249 : !> psi(stp) = f(stp) - f(0) - ftol*stp*f'(0).
3250 : !>
3251 : !> If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
3252 : !> interval is chosen so that it contains a minimizer of f.
3253 : !>
3254 : !> The algorithm is designed to find a step that satisfies
3255 : !> the sufficient decrease condition
3256 : !>
3257 : !> f(stp) <= f(0) + ftol*stp*f'(0),
3258 : !>
3259 : !> and the curvature condition
3260 : !>
3261 : !> abs(f'(stp)) <= gtol*abs(f'(0)).
3262 : !>
3263 : !> If ftol is less than gtol and if, for example, the function
3264 : !> is bounded below, then there is always a step which satisfies
3265 : !> both conditions.
3266 : !>
3267 : !> If no step can be found that satisfies both conditions, then
3268 : !> the algorithm stops with a warning. In this case stp only
3269 : !> satisfies the sufficient decrease condition.
3270 : !>
3271 : !> A typical invocation of dcsrch has the following outline:
3272 : !>
3273 : !> task = 'START'
3274 : !> DO WHILE (.TRUE.)
3275 : !> call dcsrch( ... )
3276 : !> if (task .eq. 'FG') then
3277 : !> Evaluate the function and the gradient at stp
3278 : !> else
3279 : !> exit
3280 : !> end if
3281 : !> END DO
3282 : !> \param f On initial entry f is the value of the function at 0.
3283 : !> On subsequent entries f is the value of the
3284 : !> function at stp.
3285 : !> On exit f is the value of the function at stp.
3286 : !> \param g On initial entry g is the derivative of the function at 0.
3287 : !> On subsequent entries g is the derivative of the
3288 : !> function at stp.
3289 : !> On exit g is the derivative of the function at stp.
3290 : !> \param stp On entry stp is the current estimate of a satisfactory
3291 : !> step. On initial entry, a positive initial estimate
3292 : !> must be provided.
3293 : !> On exit stp is the current estimate of a satisfactory step
3294 : !> if task = 'FG'. If task = 'CONV' then stp satisfies
3295 : !> the sufficient decrease and curvature condition.
3296 : !> \param ftol ftol specifies a nonnegative tolerance for the
3297 : !> sufficient decrease condition.
3298 : !> \param gtol gtol specifies a nonnegative tolerance for the
3299 : !> curvature condition.
3300 : !> \param xtol xtol specifies a nonnegative relative tolerance
3301 : !> for an acceptable step. The subroutine exits with a
3302 : !> warning if the relative difference between sty and stx
3303 : !> is less than xtol.
3304 : !> \param stpmin stpmin is a nonnegative lower bound for the step.
3305 : !> \param stpmax stpmax is a nonnegative upper bound for the step.
3306 : !> \param task task is a character variable of length at least 60.
3307 : !> On initial entry task must be set to 'START'.
3308 : !> On exit task indicates the required action:
3309 : !>
3310 : !> If task(1:2) = 'FG' then evaluate the function and
3311 : !> derivative at stp and call dcsrch again.
3312 : !>
3313 : !> If task(1:4) = 'CONV' then the search is successful.
3314 : !>
3315 : !> If task(1:4) = 'WARN' then the subroutine is not able
3316 : !> to satisfy the convergence conditions. The exit value of
3317 : !> stp contains the best point found during the search.
3318 : !>
3319 : !> If task(1:5) = 'ERROR' then there is an error in the
3320 : !> input arguments.
3321 : !>
3322 : !> On exit with convergence, a warning or an error, the
3323 : !> variable task contains additional information.
3324 : !> \param isave is work array
3325 : !> \param dsave is a work array
3326 : ! **************************************************************************************************
3327 3135 : SUBROUTINE dcsrch(f, g, stp, ftol, gtol, xtol, stpmin, stpmax, &
3328 : task, isave, dsave)
3329 : REAL(KIND=dp) :: f, g
3330 : REAL(KIND=dp), INTENT(inout) :: stp
3331 : REAL(KIND=dp) :: ftol, gtol, xtol, stpmin, stpmax
3332 : CHARACTER(LEN=*) :: task
3333 : INTEGER :: isave(2)
3334 : REAL(KIND=dp) :: dsave(13)
3335 :
3336 : REAL(KIND=dp), PARAMETER :: p5 = 0.5_dp, p66 = 0.66_dp, &
3337 : xtrapl = 1.1_dp, xtrapu = 4.0_dp, &
3338 : zero = 0.0_dp
3339 :
3340 : INTEGER :: stage
3341 : LOGICAL :: brackt
3342 : REAL(KIND=dp) :: finit, fm, ftest, fx, fxm, fy, fym, &
3343 : ginit, gm, gtest, gx, gxm, gy, gym, &
3344 : stmax, stmin, stx, sty, width, width1
3345 :
3346 : !
3347 : ! NOTE: The user must no alter work arrays between calls.
3348 : !
3349 : !
3350 : ! MINPACK-1 Project. June 1983.
3351 : ! Argonne National Laboratory.
3352 : ! Jorge J. More' and David J. Thuente.
3353 : !
3354 : ! MINPACK-2 Project. October 1993.
3355 : ! Argonne National Laboratory and University of Minnesota.
3356 : ! Brett M. Averick, Richard G. Carter, and Jorge J. More'.
3357 : !
3358 : ! **********
3359 : ! Initialization block.
3360 :
3361 3135 : IF (task(1:5) == 'START') THEN
3362 :
3363 : ! Check the input arguments for errors.
3364 :
3365 1484 : IF (stp < stpmin) task = 'ERROR: STP < STPMIN'
3366 1484 : IF (stp > stpmax) task = 'ERROR: STP > STPMAX'
3367 1484 : IF (g >= zero) task = 'ERROR: INITIAL G >= ZERO'
3368 1484 : IF (ftol < zero) task = 'ERROR: FTOL < ZERO'
3369 1484 : IF (gtol < zero) task = 'ERROR: GTOL < ZERO'
3370 1484 : IF (xtol < zero) task = 'ERROR: XTOL < ZERO'
3371 1484 : IF (stpmin < zero) task = 'ERROR: STPMIN < ZERO'
3372 1484 : IF (stpmax < stpmin) task = 'ERROR: STPMAX < STPMIN'
3373 :
3374 : ! Exit if there are errors on input.
3375 :
3376 1484 : IF (task(1:5) == 'ERROR') RETURN
3377 :
3378 : ! Initialize local variables.
3379 :
3380 1484 : brackt = .FALSE.
3381 1484 : stage = 1
3382 1484 : finit = f
3383 1484 : ginit = g
3384 1484 : gtest = ftol*ginit
3385 1484 : width = stpmax - stpmin
3386 1484 : width1 = width/p5
3387 :
3388 : ! The variables stx, fx, gx contain the values of the step,
3389 : ! function, and derivative at the best step.
3390 : ! The variables sty, fy, gy contain the value of the step,
3391 : ! function, and derivative at sty.
3392 : ! The variables stp, f, g contain the values of the step,
3393 : ! function, and derivative at stp.
3394 :
3395 1484 : stx = zero
3396 1484 : fx = finit
3397 1484 : gx = ginit
3398 1484 : sty = zero
3399 1484 : fy = finit
3400 1484 : gy = ginit
3401 1484 : stmin = zero
3402 1484 : stmax = stp + xtrapu*stp
3403 1484 : task = 'FG'
3404 :
3405 : ELSE
3406 :
3407 : ! Restore local variables.
3408 :
3409 1651 : IF (isave(1) == 1) THEN
3410 154 : brackt = .TRUE.
3411 : ELSE
3412 1497 : brackt = .FALSE.
3413 : END IF
3414 1651 : stage = isave(2)
3415 1651 : ginit = dsave(1)
3416 1651 : gtest = dsave(2)
3417 1651 : gx = dsave(3)
3418 1651 : gy = dsave(4)
3419 1651 : finit = dsave(5)
3420 1651 : fx = dsave(6)
3421 1651 : fy = dsave(7)
3422 1651 : stx = dsave(8)
3423 1651 : sty = dsave(9)
3424 1651 : stmin = dsave(10)
3425 1651 : stmax = dsave(11)
3426 1651 : width = dsave(12)
3427 1651 : width1 = dsave(13)
3428 :
3429 : ! If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
3430 : ! algorithm enters the second stage.
3431 :
3432 1651 : ftest = finit + stp*gtest
3433 1651 : IF (stage == 1 .AND. f <= ftest .AND. g >= zero) THEN
3434 329 : stage = 2
3435 : END IF
3436 :
3437 : ! Test for warnings.
3438 :
3439 1651 : IF (brackt .AND. (stp <= stmin .OR. stp >= stmax)) THEN
3440 3 : task = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS'
3441 : END IF
3442 1651 : IF (brackt .AND. stmax - stmin <= xtol*stmax) THEN
3443 2 : task = 'WARNING: XTOL TEST SATISFIED'
3444 : END IF
3445 1651 : IF (stp == stpmax .AND. f <= ftest .AND. g <= gtest) THEN
3446 11 : task = 'WARNING: STP = STPMAX'
3447 : END IF
3448 1651 : IF (stp == stpmin .AND. (f > ftest .OR. g >= gtest)) THEN
3449 1 : task = 'WARNING: STP = STPMIN'
3450 : END IF
3451 :
3452 : ! Test for convergence.
3453 :
3454 1651 : IF (f <= ftest .AND. ABS(g) <= gtol*(-ginit)) THEN
3455 1476 : task = 'CONVERGENCE'
3456 : END IF
3457 :
3458 : ! Test for termination.
3459 :
3460 1651 : IF (.NOT. (task(1:4) == 'WARN' .OR. task(1:4) == 'CONV')) THEN
3461 :
3462 : ! A modified function is used to predict the step during the
3463 : ! first stage if a lower function value has been obtained but
3464 : ! the decrease is not sufficient.
3465 :
3466 167 : IF (stage == 1 .AND. f <= fx .AND. f > ftest) THEN
3467 :
3468 : ! Define the modified function and derivative values.
3469 :
3470 0 : fm = f - stp*gtest
3471 0 : fxm = fx - stx*gtest
3472 0 : fym = fy - sty*gtest
3473 0 : gm = g - gtest
3474 0 : gxm = gx - gtest
3475 0 : gym = gy - gtest
3476 :
3477 : ! Call dcstep to update stx, sty, and to compute the new step.
3478 :
3479 : CALL dcstep(stx, fxm, gxm, sty, fym, gym, stp, fm, gm, &
3480 0 : brackt, stmin, stmax)
3481 :
3482 : ! Reset the function and derivative values for f.
3483 :
3484 0 : fx = fxm + stx*gtest
3485 0 : fy = fym + sty*gtest
3486 0 : gx = gxm + gtest
3487 0 : gy = gym + gtest
3488 :
3489 : ELSE
3490 :
3491 : ! Call dcstep to update stx, sty, and to compute the new step.
3492 :
3493 : CALL dcstep(stx, fx, gx, sty, fy, gy, stp, f, g, &
3494 167 : brackt, stmin, stmax)
3495 :
3496 : END IF
3497 :
3498 : ! Decide if a bisection step is needed.
3499 :
3500 167 : IF (brackt) THEN
3501 154 : IF (ABS(sty - stx) >= p66*width1) stp = stx + p5*(sty - stx)
3502 154 : width1 = width
3503 154 : width = ABS(sty - stx)
3504 : END IF
3505 :
3506 : ! Set the minimum and maximum steps allowed for stp.
3507 :
3508 167 : IF (brackt) THEN
3509 154 : stmin = MIN(stx, sty)
3510 154 : stmax = MAX(stx, sty)
3511 : ELSE
3512 13 : stmin = stp + xtrapl*(stp - stx)
3513 13 : stmax = stp + xtrapu*(stp - stx)
3514 : END IF
3515 :
3516 : ! Force the step to be within the bounds stpmax and stpmin.
3517 :
3518 167 : stp = MAX(stp, stpmin)
3519 167 : stp = MIN(stp, stpmax)
3520 :
3521 : ! If further progress is not possible, let stp be the best
3522 : ! point obtained during the search.
3523 :
3524 : IF (brackt .AND. (stp <= stmin .OR. stp >= stmax) &
3525 167 : .OR. (brackt .AND. stmax - stmin <= xtol*stmax)) stp = stx
3526 :
3527 : ! Obtain another function and derivative.
3528 :
3529 167 : task = 'FG'
3530 :
3531 : END IF
3532 : END IF
3533 :
3534 : ! Save local variables.
3535 :
3536 3135 : IF (brackt) THEN
3537 273 : isave(1) = 1
3538 : ELSE
3539 2862 : isave(1) = 0
3540 : END IF
3541 3135 : isave(2) = stage
3542 3135 : dsave(1) = ginit
3543 3135 : dsave(2) = gtest
3544 3135 : dsave(3) = gx
3545 3135 : dsave(4) = gy
3546 3135 : dsave(5) = finit
3547 3135 : dsave(6) = fx
3548 3135 : dsave(7) = fy
3549 3135 : dsave(8) = stx
3550 3135 : dsave(9) = sty
3551 3135 : dsave(10) = stmin
3552 3135 : dsave(11) = stmax
3553 3135 : dsave(12) = width
3554 3135 : dsave(13) = width1
3555 :
3556 3135 : RETURN
3557 3135 : END SUBROUTINE dcsrch
3558 :
3559 : ! **************************************************************************************************
3560 : !> \brief This subroutine computes a safeguarded step for a search
3561 : !> procedure and updates an interval that contains a step that
3562 : !> satisfies a sufficient decrease and a curvature condition.
3563 : !>
3564 : !> The parameter stx contains the step with the least function
3565 : !> value. If brackt is set to .true. then a minimizer has
3566 : !> been bracketed in an interval with endpoints stx and sty.
3567 : !> The parameter stp contains the current step.
3568 : !> The subroutine assumes that if brackt is set to .true. then
3569 : !>
3570 : !> min(stx,sty) < stp < max(stx,sty),
3571 : !>
3572 : !> and that the derivative at stx is negative in the direction
3573 : !> of the step.
3574 : !> \param stx On entry stx is the best step obtained so far and is an
3575 : !> endpoint of the interval that contains the minimizer.
3576 : !> On exit stx is the updated best step.
3577 : !> \param fx fx is the function at stx.
3578 : !> \param dx On entry dx is the derivative of the function at
3579 : !> stx. The derivative must be negative in the direction of
3580 : !> the step, that is, dx and stp - stx must have opposite
3581 : !> signs.
3582 : !> On exit dx is the derivative of the function at stx.
3583 : !> \param sty On entry sty is the second endpoint of the interval that
3584 : !> contains the minimizer.
3585 : !> On exit sty is the updated endpoint of the interval that
3586 : !> contains the minimizer.
3587 : !> \param fy fy is the function at sty.
3588 : !> \param dy On entry dy is the derivative of the function at sty.
3589 : !> On exit dy is the derivative of the function at the exit sty.
3590 : !> \param stp On entry stp is the current step. If brackt is set to .true.
3591 : !> then on input stp must be between stx and sty.
3592 : !> On exit stp is a new trial step.
3593 : !> \param fp fp is the function at stp
3594 : !> \param dp_loc dp_loc is the the derivative of the function at stp.
3595 : !> \param brackt On entry brackt specifies if a minimizer has been bracketed.
3596 : !> Initially brackt must be set to .false.
3597 : !> On exit brackt specifies if a minimizer has been bracketed.
3598 : !> When a minimizer is bracketed brackt is set to .true.
3599 : !> \param stpmin stpmin is a lower bound for the step.
3600 : !> \param stpmax stpmax is an upper bound for the step.
3601 : ! **************************************************************************************************
3602 167 : SUBROUTINE dcstep(stx, fx, dx, sty, fy, dy, stp, fp, dp_loc, brackt, &
3603 : stpmin, stpmax)
3604 : REAL(KIND=dp), INTENT(inout) :: stx, fx, dx, sty, fy, dy, stp
3605 : REAL(KIND=dp), INTENT(in) :: fp, dp_loc
3606 : LOGICAL, INTENT(inout) :: brackt
3607 : REAL(KIND=dp), INTENT(in) :: stpmin, stpmax
3608 :
3609 : REAL(KIND=dp), PARAMETER :: p66 = 0.66_dp, three = 3.0_dp, &
3610 : two = 2.0_dp, zero = 0.0_dp
3611 :
3612 : REAL(KIND=dp) :: gamma, p, q, r, s, sgnd, stpc, stpf, &
3613 : stpq, theta
3614 :
3615 : !
3616 : ! MINPACK-1 Project. June 1983
3617 : ! Argonne National Laboratory.
3618 : ! Jorge J. More' and David J. Thuente.
3619 : !
3620 : ! MINPACK-2 Project. October 1993.
3621 : ! Argonne National Laboratory and University of Minnesota.
3622 : ! Brett M. Averick and Jorge J. More'.
3623 : !
3624 : ! **********
3625 :
3626 167 : sgnd = dp_loc*SIGN(1.0_dp, dx)
3627 :
3628 : ! First case: A higher function value. The minimum is bracketed.
3629 : ! If the cubic step is closer to stx than the quadratic step, the
3630 : ! cubic step is taken, otherwise the average of the cubic and
3631 : ! quadratic steps is taken.
3632 :
3633 167 : IF (fp > fx) THEN
3634 123 : theta = three*(fx - fp)/(stp - stx) + dx + dp_loc
3635 123 : s = MAX(ABS(theta), ABS(dx), ABS(dp_loc))
3636 123 : gamma = s*SQRT((theta/s)**2 - (dx/s)*(dp_loc/s))
3637 123 : IF (stp < stx) gamma = -gamma
3638 123 : p = (gamma - dx) + theta
3639 123 : q = ((gamma - dx) + gamma) + dp_loc
3640 123 : r = p/q
3641 123 : stpc = stx + r*(stp - stx)
3642 : stpq = stx + ((dx/((fx - fp)/(stp - stx) + dx))/two)* &
3643 123 : & (stp - stx)
3644 123 : IF (ABS(stpc - stx) < ABS(stpq - stx)) THEN
3645 : stpf = stpc
3646 : ELSE
3647 55 : stpf = stpc + (stpq - stpc)/two
3648 : END IF
3649 123 : brackt = .TRUE.
3650 :
3651 : ! Second case: A lower function value and derivatives of opposite
3652 : ! sign. The minimum is bracketed. If the cubic step is farther from
3653 : ! stp than the secant step, the cubic step is taken, otherwise the
3654 : ! secant step is taken.
3655 :
3656 44 : ELSE IF (sgnd < zero) THEN
3657 15 : theta = three*(fx - fp)/(stp - stx) + dx + dp_loc
3658 15 : s = MAX(ABS(theta), ABS(dx), ABS(dp_loc))
3659 15 : gamma = s*SQRT((theta/s)**2 - (dx/s)*(dp_loc/s))
3660 15 : IF (stp > stx) gamma = -gamma
3661 15 : p = (gamma - dp_loc) + theta
3662 15 : q = ((gamma - dp_loc) + gamma) + dx
3663 15 : r = p/q
3664 15 : stpc = stp + r*(stx - stp)
3665 15 : stpq = stp + (dp_loc/(dp_loc - dx))*(stx - stp)
3666 15 : IF (ABS(stpc - stp) > ABS(stpq - stp)) THEN
3667 : stpf = stpc
3668 : ELSE
3669 13 : stpf = stpq
3670 : END IF
3671 15 : brackt = .TRUE.
3672 :
3673 : ! Third case: A lower function value, derivatives of the same sign,
3674 : ! and the magnitude of the derivative decreases.
3675 :
3676 29 : ELSE IF (ABS(dp_loc) < ABS(dx)) THEN
3677 :
3678 : ! The cubic step is computed only if the cubic tends to infinity
3679 : ! in the direction of the step or if the minimum of the cubic
3680 : ! is beyond stp. Otherwise the cubic step is defined to be the
3681 : ! secant step.
3682 :
3683 8 : theta = three*(fx - fp)/(stp - stx) + dx + dp_loc
3684 8 : s = MAX(ABS(theta), ABS(dx), ABS(dp_loc))
3685 :
3686 : ! The case gamma = 0 only arises if the cubic does not tend
3687 : ! to infinity in the direction of the step.
3688 :
3689 8 : gamma = s*SQRT(MAX(zero, (theta/s)**2 - (dx/s)*(dp_loc/s)))
3690 8 : IF (stp > stx) gamma = -gamma
3691 8 : p = (gamma - dp_loc) + theta
3692 8 : q = (gamma + (dx - dp_loc)) + gamma
3693 8 : r = p/q
3694 8 : IF (r < zero .AND. gamma /= zero) THEN
3695 8 : stpc = stp + r*(stx - stp)
3696 0 : ELSE IF (stp > stx) THEN
3697 0 : stpc = stpmax
3698 : ELSE
3699 0 : stpc = stpmin
3700 : END IF
3701 8 : stpq = stp + (dp_loc/(dp_loc - dx))*(stx - stp)
3702 :
3703 8 : IF (brackt) THEN
3704 :
3705 : ! A minimizer has been bracketed. If the cubic step is
3706 : ! closer to stp than the secant step, the cubic step is
3707 : ! taken, otherwise the secant step is taken.
3708 :
3709 4 : IF (ABS(stpc - stp) < ABS(stpq - stp)) THEN
3710 : stpf = stpc
3711 : ELSE
3712 0 : stpf = stpq
3713 : END IF
3714 4 : IF (stp > stx) THEN
3715 4 : stpf = MIN(stp + p66*(sty - stp), stpf)
3716 : ELSE
3717 0 : stpf = MAX(stp + p66*(sty - stp), stpf)
3718 : END IF
3719 : ELSE
3720 :
3721 : ! A minimizer has not been bracketed. If the cubic step is
3722 : ! farther from stp than the secant step, the cubic step is
3723 : ! taken, otherwise the secant step is taken.
3724 :
3725 4 : IF (ABS(stpc - stp) > ABS(stpq - stp)) THEN
3726 : stpf = stpc
3727 : ELSE
3728 3 : stpf = stpq
3729 : END IF
3730 4 : stpf = MIN(stpmax, stpf)
3731 4 : stpf = MAX(stpmin, stpf)
3732 : END IF
3733 :
3734 : ! Fourth case: A lower function value, derivatives of the same sign,
3735 : ! and the magnitude of the derivative does not decrease. If the
3736 : ! minimum is not bracketed, the step is either stpmin or stpmax,
3737 : ! otherwise the cubic step is taken.
3738 :
3739 : ELSE
3740 21 : IF (brackt) THEN
3741 12 : theta = three*(fp - fy)/(sty - stp) + dy + dp_loc
3742 12 : s = MAX(ABS(theta), ABS(dy), ABS(dp_loc))
3743 12 : gamma = s*SQRT((theta/s)**2 - (dy/s)*(dp_loc/s))
3744 12 : IF (stp > sty) gamma = -gamma
3745 12 : p = (gamma - dp_loc) + theta
3746 12 : q = ((gamma - dp_loc) + gamma) + dy
3747 12 : r = p/q
3748 12 : stpc = stp + r*(sty - stp)
3749 12 : stpf = stpc
3750 9 : ELSE IF (stp > stx) THEN
3751 9 : stpf = stpmax
3752 : ELSE
3753 0 : stpf = stpmin
3754 : END IF
3755 : END IF
3756 :
3757 : ! Update the interval which contains a minimizer.
3758 :
3759 167 : IF (fp > fx) THEN
3760 123 : sty = stp
3761 123 : fy = fp
3762 123 : dy = dp_loc
3763 : ELSE
3764 44 : IF (sgnd < zero) THEN
3765 15 : sty = stx
3766 15 : fy = fx
3767 15 : dy = dx
3768 : END IF
3769 44 : stx = stp
3770 44 : fx = fp
3771 44 : dx = dp_loc
3772 : END IF
3773 :
3774 : ! Compute the new step.
3775 :
3776 167 : stp = stpf
3777 :
3778 167 : RETURN
3779 : END SUBROUTINE dcstep
3780 :
3781 : !MK LINPACK
3782 :
3783 : ! **************************************************************************************************
3784 : !> \brief factors a double precision symmetric positive definite
3785 : !> matrix.
3786 : !>
3787 : !> dpofa is usually called by dpoco, but it can be called
3788 : !> directly with a saving in time if rcond is not needed.
3789 : !> (time for dpoco) = (1 + 18/n)*(time for dpofa) .
3790 : !> \param a the symmetric matrix to be factored. only the
3791 : !> diagonal and upper triangle are used.
3792 : !> on return
3793 : !> an upper triangular matrix r so that a = trans(r)*r
3794 : !> where trans(r) is the transpose.
3795 : !> the strict lower triangle is unaltered.
3796 : !> if info .ne. 0 , the factorization is not complete.
3797 : !> \param lda the leading dimension of the array a .
3798 : !> \param n the order of the matrix a .
3799 : !> \param info = 0 for normal return.
3800 : !> = k signals an error condition. the leading minor
3801 : !> of order k is not positive definite.
3802 : ! **************************************************************************************************
3803 4308 : SUBROUTINE dpofa(a, lda, n, info)
3804 : INTEGER, INTENT(in) :: lda
3805 : REAL(KIND=dp) :: a(lda, *)
3806 : INTEGER, INTENT(in) :: n
3807 : INTEGER :: info
3808 :
3809 : INTEGER :: j, jm1, k
3810 : REAL(KIND=dp) :: ddot, s, t
3811 :
3812 : !
3813 : ! linpack. this version dated 08/14/78 .
3814 : ! cleve moler, university of new mexico, argonne national lab.
3815 : !
3816 : ! begin block with ...exits to 40
3817 : !
3818 : !
3819 :
3820 25968 : DO j = 1, n
3821 21660 : info = j
3822 21660 : s = 0.0_dp
3823 21660 : jm1 = j - 1
3824 21660 : IF (.NOT. (jm1 < 1)) THEN
3825 68421 : DO k = 1, jm1
3826 51069 : t = a(k, j) - ddot(k - 1, a(1, k), 1, a(1, j), 1)
3827 51069 : t = t/a(k, k)
3828 51069 : a(k, j) = t
3829 68421 : s = s + t*t
3830 : END DO
3831 : END IF
3832 21660 : s = a(j, j) - s
3833 : ! ......exit
3834 21660 : IF (s <= 0.0_dp) EXIT
3835 21660 : a(j, j) = SQRT(s)
3836 25968 : info = 0
3837 : END DO
3838 4308 : RETURN
3839 : END SUBROUTINE dpofa
3840 :
3841 : ! **************************************************************************************************
3842 : !> \brief dtrsl solves systems of the form
3843 : !>
3844 : !> t * x = b
3845 : !> or
3846 : !> trans(t) * x = b
3847 : !>
3848 : !> where t is a triangular matrix of order n. here trans(t)
3849 : !> denotes the transpose of the matrix t.
3850 : !> \param t t contains the matrix of the system. the zero
3851 : !> elements of the matrix are not referenced, and
3852 : !> the corresponding elements of the array can be
3853 : !> used to store other information.
3854 : !> \param ldt ldt is the leading dimension of the array t.
3855 : !> \param n n is the order of the system.
3856 : !> \param b contains the right hand side of the system.
3857 : !> on return
3858 : !> b contains the solution, if info .eq. 0.
3859 : !> otherwise b is unaltered.
3860 : !> \param job job specifies what kind of system is to be solved.
3861 : !> if job is
3862 : !> 00 solve t*x=b, t lower triangular,
3863 : !> 01 solve t*x=b, t upper triangular,
3864 : !> 10 solve trans(t)*x=b, t lower triangular,
3865 : !> 11 solve trans(t)*x=b, t upper triangular.
3866 : !> \param info on return
3867 : !> info contains zero if the system is nonsingular.
3868 : !> otherwise info contains the index of
3869 : !> the first zero diagonal element of t.
3870 : ! **************************************************************************************************
3871 10160 : SUBROUTINE dtrsl(t, ldt, n, b, job, info)
3872 : INTEGER, INTENT(in) :: ldt
3873 : REAL(KIND=dp), INTENT(in) :: t(ldt, *)
3874 : INTEGER, INTENT(in) :: n
3875 : REAL(KIND=dp), INTENT(inout) :: b(*)
3876 : INTEGER, INTENT(in) :: job
3877 : INTEGER, INTENT(out) :: info
3878 :
3879 : INTEGER :: CASE, j, jj
3880 : REAL(KIND=dp) :: ddot, temp
3881 :
3882 : ! linpack. this version dated 08/14/78 .
3883 : ! g. w. stewart, university of maryland, argonne national lab.
3884 : !
3885 : ! begin block permitting ...exits to 150
3886 : !
3887 : ! check for zero diagonal elements.
3888 : !
3889 :
3890 80482 : DO info = 1, n
3891 : ! ......exit
3892 80482 : IF (t(info, info) == 0.0_dp) RETURN
3893 : END DO
3894 10160 : info = 0
3895 : !
3896 : ! determine the task and go to it.
3897 : !
3898 10160 : CASE = 1
3899 10160 : IF (MOD(job, 10) /= 0) CASE = 2
3900 10160 : IF (MOD(job, 100)/10 /= 0) CASE = CASE + 2
3901 :
3902 0 : SELECT CASE (CASE)
3903 : CASE (1)
3904 : !
3905 : ! solve t*x=b for t lower triangular
3906 : !
3907 0 : b(1) = b(1)/t(1, 1)
3908 0 : IF (n > 1) THEN
3909 0 : DO j = 2, n
3910 0 : temp = -b(j - 1)
3911 0 : CALL daxpy(n - j + 1, temp, t(j, j - 1), 1, b(j), 1)
3912 0 : b(j) = b(j)/t(j, j)
3913 : END DO
3914 : END IF
3915 : CASE (2)
3916 : !
3917 : ! solve t*x=b for t upper triangular.
3918 : !
3919 1470 : b(n) = b(n)/t(n, n)
3920 1470 : IF (n > 1) THEN
3921 14518 : DO jj = 2, n
3922 13058 : j = n - jj + 1
3923 13058 : temp = -b(j + 1)
3924 13058 : CALL daxpy(j, temp, t(1, j + 1), 1, b(1), 1)
3925 14518 : b(j) = b(j)/t(j, j)
3926 : END DO
3927 : END IF
3928 : CASE (3)
3929 : !
3930 : ! solve trans(t)*x=b for t lower triangular.
3931 : !
3932 0 : b(n) = b(n)/t(n, n)
3933 0 : IF (n > 1) THEN
3934 0 : DO jj = 2, n
3935 0 : j = n - jj + 1
3936 0 : b(j) = b(j) - ddot(jj - 1, t(j + 1, j), 1, b(j + 1), 1)
3937 0 : b(j) = b(j)/t(j, j)
3938 : END DO
3939 : END IF
3940 : CASE (4)
3941 : !
3942 : ! solve trans(t)*x=b for t upper triangular.
3943 : !
3944 8690 : b(1) = b(1)/t(1, 1)
3945 8690 : IF (.NOT. (n < 2)) THEN
3946 55741 : DO j = 2, n
3947 47104 : b(j) = b(j) - ddot(j - 1, t(1, j), 1, b(1), 1)
3948 55741 : b(j) = b(j)/t(j, j)
3949 : END DO
3950 : END IF
3951 : CASE DEFAULT
3952 10160 : CPABORT("unexpected case")
3953 : END SELECT
3954 :
3955 : RETURN
3956 : END SUBROUTINE dtrsl
3957 :
3958 : !MK Timer
3959 :
3960 : ! **************************************************************************************************
3961 : !> \brief This routine computes cpu time in double precision; it makes use o
3962 : !> the intrinsic f90 cpu_time therefore a conversion type is
3963 : !> needed.
3964 : !> \param ttime ...
3965 : ! **************************************************************************************************
3966 6008 : SUBROUTINE timer(ttime)
3967 : REAL(KIND=dp) :: ttime
3968 :
3969 : !
3970 : ! REAL temp
3971 : !
3972 : ! J.L Morales Departamento de Matematicas,
3973 : ! Instituto Tecnologico Autonomo de Mexico
3974 : ! Mexico D.F.
3975 : !
3976 : ! J.L Nocedal Department of Electrical Engineering and
3977 : ! Computer Science.
3978 : ! Northwestern University. Evanston, IL. USA
3979 : !
3980 : ! January 21, 2011
3981 : !
3982 : !MK temp = sngl(ttime)
3983 : !MK CALL cpu_time(temp)
3984 : !MK ttime = REAL(temp, KIND=dp)
3985 :
3986 6008 : ttime = m_walltime()
3987 :
3988 6008 : END SUBROUTINE timer
3989 :
3990 : ! **************************************************************************************************
3991 : !> \brief Saves the lcoal variables, long term this should be replaces by a lbfgs type
3992 : !> \param lsave lsave is a working array
3993 : !> On exit with 'task' = NEW_X, the following information is available:
3994 : !> If lsave(1) = .true. then the initial X has been replaced by
3995 : !> its projection in the feasible set
3996 : !> If lsave(2) = .true. then the problem is constrained;
3997 : !> If lsave(3) = .true. then each variable has upper and lower bounds;
3998 : !> \param isave isave is a working array
3999 : !> On exit with 'task' = NEW_X, the following information is available:
4000 : !> isave(22) = the total number of intervals explored in the
4001 : !> search of Cauchy points;
4002 : !> isave(26) = the total number of skipped BFGS updates before the current iteration;
4003 : !> isave(30) = the number of current iteration;
4004 : !> isave(31) = the total number of BFGS updates prior the current iteration;
4005 : !> isave(33) = the number of intervals explored in the search of
4006 : !> Cauchy point in the current iteration;
4007 : !> isave(34) = the total number of function and gradient evaluations;
4008 : !> isave(36) = the number of function value or gradient
4009 : !> evaluations in the current iteration;
4010 : !> if isave(37) = 0 then the subspace argmin is within the box;
4011 : !> if isave(37) = 1 then the subspace argmin is beyond the box;
4012 : !> isave(38) = the number of free variables in the current iteration;
4013 : !> isave(39) = the number of active constraints in the current iteration;
4014 : !> n + 1 - isave(40) = the number of variables leaving the set of
4015 : !> active constraints in the current iteration;
4016 : !> isave(41) = the number of variables entering the set of active
4017 : !> constraints in the current iteration.
4018 : !> \param dsave dsave is a working array of dimension 29.
4019 : !> On exit with 'task' = NEW_X, the following information is available:
4020 : !> dsave(1) = current 'theta' in the BFGS matrix;
4021 : !> dsave(2) = f(x) in the previous iteration;
4022 : !> dsave(3) = factr*epsmch;
4023 : !> dsave(4) = 2-norm of the line search direction vector;
4024 : !> dsave(5) = the machine precision epsmch generated by the code;
4025 : !> dsave(7) = the accumulated time spent on searching for Cauchy points;
4026 : !> dsave(8) = the accumulated time spent on subspace minimization;
4027 : !> dsave(9) = the accumulated time spent on line search;
4028 : !> dsave(11) = the slope of the line search function at the current point of line search;
4029 : !> dsave(12) = the maximum relative step length imposed in line search;
4030 : !> dsave(13) = the infinity norm of the projected gradient;
4031 : !> dsave(14) = the relative step length in the line search;
4032 : !> dsave(15) = the slope of the line search function at the starting point of the line search;
4033 : !> dsave(16) = the square of the 2-norm of the line search direction vector.
4034 : !> \param x_projected ...
4035 : !> \param constrained ...
4036 : !> \param boxed ...
4037 : !> \param updatd ...
4038 : !> \param nintol ...
4039 : !> \param itfile ...
4040 : !> \param iback ...
4041 : !> \param nskip ...
4042 : !> \param head ...
4043 : !> \param col ...
4044 : !> \param itail ...
4045 : !> \param iter ...
4046 : !> \param iupdat ...
4047 : !> \param nseg ...
4048 : !> \param nfgv ...
4049 : !> \param info ...
4050 : !> \param ifun ...
4051 : !> \param iword ...
4052 : !> \param nfree ...
4053 : !> \param nact ...
4054 : !> \param ileave ...
4055 : !> \param nenter ...
4056 : !> \param theta ...
4057 : !> \param fold ...
4058 : !> \param tol ...
4059 : !> \param dnorm ...
4060 : !> \param epsmch ...
4061 : !> \param cpu1 ...
4062 : !> \param cachyt ...
4063 : !> \param sbtime ...
4064 : !> \param lnscht ...
4065 : !> \param time1 ...
4066 : !> \param gd ...
4067 : !> \param step_max ...
4068 : !> \param g_inf_norm ...
4069 : !> \param stp ...
4070 : !> \param gdold ...
4071 : !> \param dtd ...
4072 : !> \author Samuel Andermatt (01.15)
4073 : ! **************************************************************************************************
4074 :
4075 3181 : SUBROUTINE save_local(lsave,isave,dsave,x_projected,constrained,boxed,updatd,nintol,itfile,iback,nskip,head,col,itail,&
4076 : iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, cpu1, &
4077 : cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
4078 : LOGICAL, INTENT(out) :: lsave(4)
4079 : INTEGER, INTENT(out) :: isave(23)
4080 : REAL(KIND=dp), INTENT(out) :: dsave(29)
4081 : LOGICAL, INTENT(in) :: x_projected, constrained, boxed, updatd
4082 : INTEGER, INTENT(in) :: nintol, itfile, iback, nskip, head, col, &
4083 : itail, iter, iupdat, nseg, nfgv, info, &
4084 : ifun, iword, nfree, nact, ileave, &
4085 : nenter
4086 : REAL(KIND=dp), INTENT(in) :: theta, fold, tol, dnorm, epsmch, cpu1, &
4087 : cachyt, sbtime, lnscht, time1, gd, &
4088 : step_max, g_inf_norm, stp, gdold, dtd
4089 :
4090 3181 : lsave(1) = x_projected
4091 3181 : lsave(2) = constrained
4092 3181 : lsave(3) = boxed
4093 3181 : lsave(4) = updatd
4094 :
4095 3181 : isave(1) = nintol
4096 3181 : isave(3) = itfile
4097 3181 : isave(4) = iback
4098 3181 : isave(5) = nskip
4099 3181 : isave(6) = head
4100 3181 : isave(7) = col
4101 3181 : isave(8) = itail
4102 3181 : isave(9) = iter
4103 3181 : isave(10) = iupdat
4104 3181 : isave(12) = nseg
4105 3181 : isave(13) = nfgv
4106 3181 : isave(14) = info
4107 3181 : isave(15) = ifun
4108 3181 : isave(16) = iword
4109 3181 : isave(17) = nfree
4110 3181 : isave(18) = nact
4111 3181 : isave(19) = ileave
4112 3181 : isave(20) = nenter
4113 :
4114 3181 : dsave(1) = theta
4115 3181 : dsave(2) = fold
4116 3181 : dsave(3) = tol
4117 3181 : dsave(4) = dnorm
4118 3181 : dsave(5) = epsmch
4119 3181 : dsave(6) = cpu1
4120 3181 : dsave(7) = cachyt
4121 3181 : dsave(8) = sbtime
4122 3181 : dsave(9) = lnscht
4123 3181 : dsave(10) = time1
4124 3181 : dsave(11) = gd
4125 3181 : dsave(12) = step_max
4126 3181 : dsave(13) = g_inf_norm
4127 3181 : dsave(14) = stp
4128 3181 : dsave(15) = gdold
4129 3181 : dsave(16) = dtd
4130 :
4131 3181 : END SUBROUTINE save_local
4132 :
4133 : END MODULE cp_lbfgs
|
