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