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