Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : !******************************************************************************
9 : MODULE powell
10 : USE kinds, ONLY: dp
11 : USE mathconstants, ONLY: twopi
12 : #include "../base/base_uses.f90"
13 :
14 : IMPLICIT NONE
15 :
16 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'powell'
17 :
18 : TYPE opt_state_type
19 : INTEGER :: state = -1
20 : INTEGER :: nvar = -1
21 : INTEGER :: iprint = -1
22 : INTEGER :: unit = -1
23 : INTEGER :: maxfun = -1
24 : REAL(dp) :: rhobeg = 0.0_dp, rhoend = 0.0_dp
25 : REAL(dp), DIMENSION(:), POINTER :: w => NULL()
26 : REAL(dp), DIMENSION(:), POINTER :: xopt => NULL()
27 : ! local variables
28 : INTEGER :: np = -1, nh = -1, nptm = -1, nftest = -1, idz = -1, itest = -1, nf = -1, nfm = -1, nfmm = -1, &
29 : nfsav = -1, knew = -1, kopt = -1, ksave = -1, ktemp = -1
30 : REAL(dp) :: rhosq = 0.0_dp, recip = 0.0_dp, reciq = 0.0_dp, fbeg = 0.0_dp, &
31 : fopt = 0.0_dp, diffa = 0.0_dp, xoptsq = 0.0_dp, &
32 : rho = 0.0_dp, delta = 0.0_dp, dsq = 0.0_dp, dnorm = 0.0_dp, &
33 : ratio = 0.0_dp, temp = 0.0_dp, tempq = 0.0_dp, beta = 0.0_dp, &
34 : dx = 0.0_dp, vquad = 0.0_dp, diff = 0.0_dp, diffc = 0.0_dp, &
35 : diffb = 0.0_dp, fsave = 0.0_dp, detrat = 0.0_dp, hdiag = 0.0_dp, &
36 : distsq = 0.0_dp, gisq = 0.0_dp, gqsq = 0.0_dp, f = 0.0_dp, &
37 : bstep = 0.0_dp, alpha = 0.0_dp, dstep = 0.0_dp
38 : END TYPE opt_state_type
39 :
40 : PRIVATE
41 : PUBLIC :: powell_optimize, opt_state_type
42 :
43 : CONTAINS
44 :
45 : ! **************************************************************************************************
46 : !> \brief ...
47 : !> \param n ...
48 : !> \param x ...
49 : !> \param optstate ...
50 : ! **************************************************************************************************
51 60386695 : SUBROUTINE powell_optimize(n, x, optstate)
52 : INTEGER, INTENT(IN) :: n
53 : REAL(dp), DIMENSION(*), INTENT(INOUT) :: x
54 : TYPE(opt_state_type), INTENT(INOUT) :: optstate
55 :
56 : CHARACTER(len=*), PARAMETER :: routineN = 'powell_optimize'
57 :
58 : INTEGER :: handle, npt
59 :
60 60386695 : CALL timeset(routineN, handle)
61 :
62 61151520 : SELECT CASE (optstate%state)
63 : CASE (0)
64 764825 : npt = 2*n + 1
65 2294475 : ALLOCATE (optstate%w((npt + 13)*(npt + n) + 3*n*(n + 3)/2))
66 2294475 : ALLOCATE (optstate%xopt(n))
67 : ! Initialize w
68 108742927 : optstate%w = 0.0_dp
69 764825 : optstate%state = 1
70 764825 : CALL newuoa(n, x, optstate)
71 : CASE (1, 2)
72 58092223 : CALL newuoa(n, x, optstate)
73 : CASE (3)
74 9 : IF (optstate%unit > 0) THEN
75 6 : WRITE (optstate%unit, *) "POWELL| Exceeding maximum number of steps"
76 : END IF
77 9 : optstate%state = -1
78 : CASE (4)
79 4 : IF (optstate%unit > 0) THEN
80 4 : WRITE (optstate%unit, *) "POWELL| Error in trust region"
81 : END IF
82 4 : optstate%state = -1
83 : CASE (5)
84 0 : IF (optstate%unit > 0) THEN
85 0 : WRITE (optstate%unit, *) "POWELL| N out of range"
86 : END IF
87 0 : optstate%state = -1
88 : CASE (6, 7)
89 764809 : optstate%state = -1
90 : CASE (8)
91 2294797 : x(1:n) = optstate%xopt(1:n)
92 764825 : DEALLOCATE (optstate%w)
93 764825 : DEALLOCATE (optstate%xopt)
94 764825 : optstate%state = -1
95 : CASE DEFAULT
96 60386695 : CPABORT("Unknown optimization state")
97 : END SELECT
98 :
99 60386695 : CALL timestop(handle)
100 :
101 60386695 : END SUBROUTINE powell_optimize
102 : ! **************************************************************************************************
103 : !> \brief ...
104 : !> \param n ...
105 : !> \param x ...
106 : !> \param optstate ...
107 : ! **************************************************************************************************
108 58857048 : SUBROUTINE newuoa(n, x, optstate)
109 :
110 : INTEGER, INTENT(IN) :: n
111 : REAL(dp), DIMENSION(*), INTENT(INOUT) :: x
112 : TYPE(opt_state_type), INTENT(INOUT) :: optstate
113 :
114 : INTEGER :: ibmat, id, ifv, igq, ihq, ipq, ivl, iw, &
115 : ixb, ixn, ixo, ixp, izmat, maxfun, &
116 : ndim, np, npt, nptm
117 : REAL(dp) :: rhobeg, rhoend
118 :
119 58857048 : maxfun = optstate%maxfun
120 58857048 : rhobeg = optstate%rhobeg
121 58857048 : rhoend = optstate%rhoend
122 :
123 : !
124 : ! This subroutine seeks the least value of a function of many variab
125 : ! by a trust region method that forms quadratic models by interpolat
126 : ! There can be some freedom in the interpolation conditions, which i
127 : ! taken up by minimizing the Frobenius norm of the change to the sec
128 : ! derivative of the quadratic model, beginning with a zero matrix. T
129 : ! arguments of the subroutine are as follows.
130 : !
131 : ! N must be set to the number of variables and must be at least two.
132 : ! NPT is the number of interpolation conditions. Its value must be i
133 : ! interval [N+2,(N+1)(N+2)/2].
134 : ! Initial values of the variables must be set in X(1),X(2),...,X(N).
135 : ! will be changed to the values that give the least calculated F.
136 : ! RHOBEG and RHOEND must be set to the initial and final values of a
137 : ! region radius, so both must be positive with RHOEND<=RHOBEG. Typ
138 : ! RHOBEG should be about one tenth of the greatest expected change
139 : ! variable, and RHOEND should indicate the accuracy that is requir
140 : ! the final values of the variables.
141 : ! The value of IPRINT should be set to 0, 1, 2 or 3, which controls
142 : ! amount of printing. Specifically, there is no output if IPRINT=0
143 : ! there is output only at the return if IPRINT=1. Otherwise, each
144 : ! value of RHO is printed, with the best vector of variables so fa
145 : ! the corresponding value of the objective function. Further, each
146 : ! value of F with its variables are output if IPRINT=3.
147 : ! MAXFUN must be set to an upper bound on the number of calls of CAL
148 : ! The array W will be used for working space. Its length must be at
149 : ! (NPT+13)*(NPT+N)+3*N*(N+3)/2.
150 : !
151 : ! SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must se
152 : ! the value of the objective function for the variables X(1),X(2),..
153 : !
154 : ! Partition the working space array, so that different parts of it c
155 : ! treated separately by the subroutine that performs the main calcul
156 : !
157 58857048 : np = n + 1
158 58857048 : npt = 2*n + 1
159 58857048 : nptm = npt - np
160 58857048 : IF (npt < n + 2 .OR. npt > ((n + 2)*np)/2) THEN
161 0 : optstate%state = 5
162 0 : RETURN
163 : END IF
164 58857048 : ndim = npt + n
165 58857048 : ixb = 1
166 58857048 : ixo = ixb + n
167 58857048 : ixn = ixo + n
168 58857048 : ixp = ixn + n
169 58857048 : ifv = ixp + n*npt
170 58857048 : igq = ifv + npt
171 58857048 : ihq = igq + n
172 58857048 : ipq = ihq + (n*np)/2
173 58857048 : ibmat = ipq + npt
174 58857048 : izmat = ibmat + ndim*n
175 58857048 : id = izmat + npt*nptm
176 58857048 : ivl = id + n
177 58857048 : iw = ivl + ndim
178 : !
179 : ! The above settings provide a partition of W for subroutine NEWUOB.
180 : ! The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements
181 : ! W plus the space that is needed by the last array of NEWUOB.
182 : !
183 : CALL newuob(n, npt, x, rhobeg, rhoend, maxfun, optstate%w(ixb:), optstate%w(ixo:), &
184 : optstate%w(ixn:), optstate%w(ixp:), optstate%w(ifv:), optstate%w(igq:), optstate%w(ihq:), &
185 : optstate%w(ipq:), optstate%w(ibmat:), optstate%w(izmat:), ndim, optstate%w(id:), &
186 58857048 : optstate%w(ivl:), optstate%w(iw:), optstate)
187 :
188 176582780 : optstate%xopt(1:n) = optstate%w(ixb:ixb + n - 1) + optstate%w(ixo:ixo + n - 1)
189 :
190 : END SUBROUTINE newuoa
191 :
192 : ! **************************************************************************************************
193 : !> \brief ...
194 : !> \param n ...
195 : !> \param npt ...
196 : !> \param x ...
197 : !> \param rhobeg ...
198 : !> \param rhoend ...
199 : !> \param maxfun ...
200 : !> \param xbase ...
201 : !> \param xopt ...
202 : !> \param xnew ...
203 : !> \param xpt ...
204 : !> \param fval ...
205 : !> \param gq ...
206 : !> \param hq ...
207 : !> \param pq ...
208 : !> \param bmat ...
209 : !> \param zmat ...
210 : !> \param ndim ...
211 : !> \param d ...
212 : !> \param vlag ...
213 : !> \param w ...
214 : !> \param opt ...
215 : ! **************************************************************************************************
216 58857048 : SUBROUTINE newuob(n, npt, x, rhobeg, rhoend, maxfun, xbase, &
217 58857048 : xopt, xnew, xpt, fval, gq, hq, pq, bmat, zmat, ndim, d, vlag, w, opt)
218 :
219 : INTEGER, INTENT(in) :: n, npt
220 : REAL(dp), DIMENSION(1:n), INTENT(inout) :: x
221 : REAL(dp), INTENT(in) :: rhobeg, rhoend
222 : INTEGER, INTENT(in) :: maxfun
223 : REAL(dp), DIMENSION(*), INTENT(inout) :: xbase, xopt, xnew
224 : REAL(dp), DIMENSION(npt, *), &
225 : INTENT(inout) :: xpt
226 : REAL(dp), DIMENSION(*), INTENT(inout) :: fval, gq, hq, pq
227 : INTEGER, INTENT(in) :: ndim
228 : REAL(dp), DIMENSION(npt, *), &
229 : INTENT(inout) :: zmat
230 : REAL(dp), DIMENSION(ndim, *), &
231 : INTENT(inout) :: bmat
232 : REAL(dp), DIMENSION(*), INTENT(inout) :: d, vlag, w
233 : TYPE(opt_state_type) :: opt
234 :
235 : INTEGER :: i, idz, ih, ip, ipt, itemp, &
236 : itest, j, jp, jpt, k, knew, &
237 : kopt, ksave, ktemp, nf, nfm, &
238 : nfmm, nfsav, nftest, nh, np, &
239 : nptm
240 : LOGICAL :: do_return, skip_check
241 : REAL(dp) :: alpha, beta, bstep, bsum, crvmin, delta, detrat, diff, diffa, &
242 : diffb, diffc, distsq, dnorm, dsq, dstep, dx, f, fbeg, fopt, fsave, &
243 : gisq, gqsq, half, hdiag, one, ratio, recip, reciq, rho, rhosq, sum, &
244 : suma, sumb, sumz, temp, tempq, tenth, vquad, xipt, xjpt, xoptsq, zero
245 :
246 : !
247 : ! The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are ide
248 : ! to the corresponding arguments in SUBROUTINE NEWUOA.
249 : ! XBASE will hold a shift of origin that should reduce the contribut
250 : ! from rounding errors to values of the model and Lagrange functio
251 : ! XOPT will be set to the displacement from XBASE of the vector of
252 : ! variables that provides the least calculated F so far.
253 : ! XNEW will be set to the displacement from XBASE of the vector of
254 : ! variables for the current calculation of F.
255 : ! XPT will contain the interpolation point coordinates relative to X
256 : ! FVAL will hold the values of F at the interpolation points.
257 : ! GQ will hold the gradient of the quadratic model at XBASE.
258 : ! HQ will hold the explicit second derivatives of the quadratic mode
259 : ! PQ will contain the parameters of the implicit second derivatives
260 : ! the quadratic model.
261 : ! BMAT will hold the last N columns of H.
262 : ! ZMAT will hold the factorization of the leading NPT by NPT submatr
263 : ! H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, wh
264 : ! the elements of DZ are plus or minus one, as specified by IDZ.
265 : ! NDIM is the first dimension of BMAT and has the value NPT+N.
266 : ! D is reserved for trial steps from XOPT.
267 : ! VLAG will contain the values of the Lagrange functions at a new po
268 : ! They are part of a product that requires VLAG to be of length ND
269 : ! The array W will be used for working space. Its length must be at
270 : ! 10*NDIM = 10*(NPT+N).
271 :
272 58857048 : skip_check = .FALSE.
273 58857048 : do_return = .FALSE.
274 :
275 58857048 : IF (opt%state == 1) THEN
276 : ! initialize all variable that will be stored
277 : np = 0
278 : nh = 0
279 : nptm = 0
280 : nftest = 0
281 764825 : idz = 0
282 764825 : itest = 0
283 764825 : nf = 0
284 764825 : nfm = 0
285 764825 : nfmm = 0
286 764825 : nfsav = 0
287 764825 : knew = 0
288 764825 : kopt = 0
289 764825 : ksave = 0
290 764825 : ktemp = 0
291 764825 : rhosq = 0._dp
292 764825 : recip = 0._dp
293 764825 : reciq = 0._dp
294 764825 : fbeg = 0._dp
295 764825 : fopt = 0._dp
296 764825 : diffa = 0._dp
297 764825 : xoptsq = 0._dp
298 764825 : rho = 0._dp
299 764825 : delta = 0._dp
300 764825 : dsq = 0._dp
301 764825 : dnorm = 0._dp
302 764825 : ratio = 0._dp
303 764825 : temp = 0._dp
304 764825 : tempq = 0._dp
305 764825 : beta = 0._dp
306 764825 : dx = 0._dp
307 764825 : vquad = 0._dp
308 764825 : diff = 0._dp
309 764825 : diffc = 0._dp
310 764825 : diffb = 0._dp
311 764825 : fsave = 0._dp
312 764825 : detrat = 0._dp
313 764825 : hdiag = 0._dp
314 764825 : distsq = 0._dp
315 764825 : gisq = 0._dp
316 764825 : gqsq = 0._dp
317 764825 : f = 0._dp
318 764825 : bstep = 0._dp
319 764825 : alpha = 0._dp
320 764825 : dstep = 0._dp
321 : !
322 : END IF
323 :
324 58857048 : ipt = 0
325 58857048 : jpt = 0
326 58857048 : xipt = 0._dp
327 58857048 : xjpt = 0._dp
328 :
329 58857048 : half = 0.5_dp
330 58857048 : one = 1.0_dp
331 58857048 : tenth = 0.1_dp
332 58857048 : zero = 0.0_dp
333 58857048 : np = n + 1
334 58857048 : nh = (n*np)/2
335 58857048 : nptm = npt - np
336 58857048 : nftest = MAX(maxfun, 1)
337 :
338 58857048 : IF (opt%state /= 2) THEN
339 : !
340 : ! Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero.
341 : !
342 2294797 : DO j = 1, n
343 1529972 : xbase(j) = x(j)
344 9211120 : DO k = 1, npt
345 9211120 : xpt(k, j) = zero
346 : END DO
347 13051533 : DO i = 1, ndim
348 12286708 : bmat(i, j) = zero
349 : END DO
350 : END DO
351 3067605 : DO ih = 1, nh
352 3067605 : hq(ih) = zero
353 : END DO
354 4589594 : DO k = 1, npt
355 3824769 : pq(k) = zero
356 12270742 : DO j = 1, nptm
357 11505917 : zmat(k, j) = zero
358 : END DO
359 : END DO
360 : !
361 : ! Begin the initialization procedure. NF becomes one more than the n
362 : ! of function values so far. The coordinates of the displacement of
363 : ! next initial interpolation point from XBASE are set in XPT(NF,.).
364 : !
365 764825 : rhosq = rhobeg*rhobeg
366 764825 : recip = one/rhosq
367 764825 : reciq = SQRT(half)/rhosq
368 764825 : nf = 0
369 764825 : CALL begin_initialization_of_interpolation
370 24194113 : RETURN
371 : END IF
372 :
373 58092223 : CALL set_state
374 :
375 58092223 : preparation: IF (nf <= npt) THEN
376 3824292 : fval(nf) = f
377 3824292 : IF (nf == 1) THEN
378 764825 : fbeg = f
379 764825 : fopt = f
380 764825 : kopt = 1
381 3059467 : ELSE IF (f < fopt) THEN
382 853903 : fopt = f
383 853903 : kopt = nf
384 : END IF
385 : !
386 : ! Set the nonzero initial elements of BMAT and the quadratic model i
387 : ! the cases when NF is at most 2*N+1.
388 : !
389 3824292 : IF (NFM <= 2*N) THEN
390 3824292 : IF (nfm >= 1 .AND. nfm <= n) THEN
391 1529766 : gq(nfm) = (f - fbeg)/rhobeg
392 1529766 : IF (npt < nf + n) THEN
393 0 : bmat(1, nfm) = -one/rhobeg
394 0 : bmat(nf, nfm) = one/rhobeg
395 0 : bmat(npt + nfm, nfm) = -half*rhosq
396 : END IF
397 2294526 : ELSE IF (nfm > n) THEN
398 1529701 : bmat(nf - n, nfmm) = half/rhobeg
399 1529701 : bmat(nf, nfmm) = -half/rhobeg
400 1529701 : zmat(1, nfmm) = -reciq - reciq
401 1529701 : zmat(nf - n, nfmm) = reciq
402 1529701 : zmat(nf, nfmm) = reciq
403 1529701 : ih = (nfmm*(nfmm + 1))/2
404 1529701 : temp = (fbeg - f)/rhobeg
405 1529701 : hq(ih) = (gq(nfmm) - temp)/rhobeg
406 1529701 : gq(nfmm) = half*(gq(nfmm) + temp)
407 : END IF
408 : !
409 : ! Set the off-diagonal second derivatives of the Lagrange functions
410 : ! the initial quadratic model.
411 : !
412 : ELSE
413 0 : ih = (ipt*(ipt - 1))/2 + jpt
414 : IF (xipt < zero) ipt = ipt + n
415 : IF (xjpt < zero) jpt = jpt + n
416 0 : zmat(1, nfmm) = recip
417 0 : zmat(nf, nfmm) = recip
418 0 : zmat(ipt + 1, nfmm) = -recip
419 : zmat(jpt + 1, nfmm) = -recip
420 0 : hq(ih) = (fbeg - fval(ipt + 1) - fval(jpt + 1) + f)/(xipt*xjpt)
421 : END IF
422 3824292 : IF (nf < npt) THEN
423 3059476 : CALL begin_initialization_of_interpolation
424 3059476 : RETURN
425 : END IF
426 : !
427 : ! Begin the iterative procedure, because the initial model is comple
428 : !
429 764816 : rho = rhobeg
430 764816 : delta = rho
431 764816 : idz = 1
432 764816 : diffa = zero
433 764816 : diffb = zero
434 764816 : itest = 0
435 764816 : xoptsq = zero
436 2294517 : DO i = 1, n
437 1529701 : xopt(i) = xpt(kopt, i)
438 2294517 : xoptsq = xoptsq + xopt(i)**2
439 : END DO
440 764816 : nfsav = nf
441 : !
442 : ! Generate the next trust region step and test its length. Set KNEW
443 : ! to -1 if the purpose of the next F will be to improve the model.
444 : !
445 764816 : skip_check = .TRUE.
446 : ELSE preparation
447 54267931 : IF (knew == -1) THEN
448 638430 : opt%state = 6
449 638430 : CALL get_state
450 638430 : CALL finalize_optimization_after_failure
451 638430 : RETURN
452 : END IF
453 : !
454 : ! Use the quadratic model to predict the change in F due to the step
455 : ! and set DIFF to the error of this prediction.
456 : !
457 53629501 : vquad = zero
458 53629501 : ih = 0
459 160894356 : DO j = 1, n
460 107264855 : vquad = vquad + d(j)*gq(j)
461 321807118 : DO i = 1, j
462 160912762 : ih = ih + 1
463 160912762 : temp = d(i)*xnew(j) + d(j)*xopt(i)
464 160912762 : IF (i == j) temp = half*temp
465 268177617 : vquad = vquad + temp*hq(ih)
466 : END DO
467 : END DO
468 321788712 : DO k = 1, npt
469 321788712 : vquad = vquad + pq(k)*w(k)
470 : END DO
471 53629501 : diff = f - fopt - vquad
472 53629501 : diffc = diffb
473 53629501 : diffb = diffa
474 53629501 : diffa = ABS(diff)
475 53629501 : IF (dnorm > rho) nfsav = nf
476 : !
477 : ! Update FOPT and XOPT if the new F is the least value of the object
478 : ! function so far. The branch when KNEW is positive occurs if D is n
479 : ! a trust region step.
480 : !
481 53629501 : fsave = fopt
482 53629501 : IF (f < fopt) THEN
483 27399973 : fopt = f
484 27399973 : xoptsq = zero
485 82201822 : DO i = 1, n
486 54801849 : xopt(i) = xnew(i)
487 82201822 : xoptsq = xoptsq + xopt(i)**2
488 : END DO
489 : END IF
490 53629501 : ksave = knew
491 53629501 : IF (knew <= 0) THEN
492 : !
493 : ! Pick the next value of DELTA after a trust region step.
494 : !
495 34662934 : IF (vquad >= zero) THEN
496 : ! Return because a trust region step has failed to reduce Q
497 4 : opt%state = 4
498 4 : CALL get_state
499 4 : CALL finalize_optimization_after_failure
500 4 : RETURN
501 : END IF
502 34662930 : ratio = (f - fsave)/vquad
503 34662930 : IF (ratio <= tenth) THEN
504 12837415 : delta = half*dnorm
505 21825515 : ELSE IF (ratio <= 0.7_dp) THEN
506 3775267 : delta = MAX(half*delta, dnorm)
507 : ELSE
508 18050248 : delta = MAX(half*delta, dnorm + dnorm)
509 : END IF
510 34662930 : IF (delta <= 1.5_dp*rho) delta = rho
511 : !
512 : ! Set KNEW to the index of the next interpolation point to be delete
513 : !
514 34662930 : rhosq = MAX(tenth*delta, rho)**2
515 34662930 : ktemp = 0
516 34662930 : detrat = zero
517 34662930 : IF (f >= fsave) THEN
518 11021104 : ktemp = kopt
519 11021104 : detrat = one
520 : END IF
521 207984594 : DO k = 1, npt
522 173321664 : hdiag = zero
523 520005361 : DO j = 1, nptm
524 346683697 : temp = one
525 346683697 : IF (j < idz) temp = -one
526 520005361 : hdiag = hdiag + temp*zmat(k, j)**2
527 : END DO
528 173321664 : temp = ABS(beta*hdiag + vlag(k)**2)
529 173321664 : distsq = zero
530 520005361 : DO j = 1, n
531 520005361 : distsq = distsq + (xpt(k, j) - xopt(j))**2
532 : END DO
533 173321664 : IF (distsq > rhosq) temp = temp*(distsq/rhosq)**3
534 207984594 : IF (temp > detrat .AND. k /= ktemp) THEN
535 74634776 : detrat = temp
536 74634776 : knew = k
537 : END IF
538 : END DO
539 34662930 : IF (knew == 0) EXIT preparation
540 : !
541 : ! Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point
542 : ! can be moved. Begin the updating of the quadratic model, starting
543 : ! with the explicit second derivative term.
544 : !
545 : END IF
546 53176237 : CALL update(n, npt, bmat, zmat, idz, ndim, vlag, beta, knew, w)
547 53176237 : fval(knew) = f
548 53176237 : ih = 0
549 159534510 : DO i = 1, n
550 106358273 : temp = pq(knew)*xpt(knew, i)
551 319087252 : DO j = 1, i
552 159552742 : ih = ih + 1
553 265911015 : hq(ih) = hq(ih) + temp*xpt(knew, j)
554 : END DO
555 : END DO
556 53176237 : pq(knew) = zero
557 : !
558 : ! Update the other second derivative parameters, and then the gradie
559 : ! vector of the model. Also include the new interpolation point.
560 : !
561 159534510 : DO j = 1, nptm
562 106358273 : temp = diff*zmat(knew, j)
563 106358273 : IF (j < idz) temp = -temp
564 691387205 : DO k = 1, npt
565 638210968 : pq(k) = pq(k) + temp*zmat(k, j)
566 : END DO
567 : END DO
568 53176237 : gqsq = zero
569 159534510 : DO i = 1, n
570 106358273 : gq(i) = gq(i) + diff*bmat(knew, i)
571 106358273 : gqsq = gqsq + gq(i)**2
572 159534510 : xpt(knew, i) = xnew(i)
573 : END DO
574 : !
575 : ! If a trust region step makes a small change to the objective funct
576 : ! then calculate the gradient of the least Frobenius norm interpolan
577 : ! XBASE, and store it in W, using VLAG for a vector of right hand si
578 : !
579 53176237 : IF (ksave == 0 .AND. delta == rho) THEN
580 6876833 : IF (ABS(ratio) > 1.0e-2_dp) THEN
581 4623844 : itest = 0
582 : ELSE
583 13517976 : DO k = 1, npt
584 13517976 : vlag(k) = fval(k) - fval(kopt)
585 : END DO
586 2252989 : gisq = zero
587 6758988 : DO i = 1, n
588 4505999 : sum = zero
589 27036208 : DO k = 1, npt
590 27036208 : sum = sum + bmat(k, i)*vlag(k)
591 : END DO
592 4505999 : gisq = gisq + sum*sum
593 6758988 : w(i) = sum
594 : END DO
595 : !
596 : ! Test whether to replace the new quadratic model by the least Frobe
597 : ! norm interpolant, making the replacement if the test is satisfied.
598 : !
599 2252989 : itest = itest + 1
600 2252989 : IF (gqsq < 1.0e2_dp*gisq) itest = 0
601 2252989 : IF (itest >= 3) THEN
602 376974 : DO i = 1, n
603 376974 : gq(i) = w(i)
604 : END DO
605 502632 : DO ih = 1, nh
606 502632 : hq(ih) = zero
607 : END DO
608 376974 : DO j = 1, nptm
609 251316 : w(j) = zero
610 1507896 : DO k = 1, npt
611 1507896 : w(j) = w(j) + vlag(k)*zmat(k, j)
612 : END DO
613 376974 : IF (j < idz) w(j) = -w(j)
614 : END DO
615 753948 : DO k = 1, npt
616 628290 : pq(k) = zero
617 2010528 : DO j = 1, nptm
618 1884870 : pq(k) = pq(k) + zmat(k, j)*w(j)
619 : END DO
620 : END DO
621 125658 : itest = 0
622 : END IF
623 : END IF
624 : END IF
625 53176237 : IF (f < fsave) kopt = knew
626 : !
627 : ! If a trust region step has provided a sufficient decrease in F, th
628 : ! branch for another trust region calculation. The case KSAVE>0 occu
629 : ! when the new function value was calculated by a model step.
630 : !
631 53176237 : IF (f <= fsave + tenth*vquad .OR. ksave > 0) THEN
632 : skip_check = .TRUE.
633 : ELSE
634 12384149 : knew = 0
635 : END IF
636 : END IF preparation
637 : !
638 : ! Alternatively, find out if the interpolation points are close enough
639 : ! to the best point so far.
640 : !
641 : outer: DO
642 63960943 : IF (.NOT. skip_check) THEN
643 22404039 : skip_check = .FALSE.
644 22404039 : distsq = 4.0_dp*delta*delta
645 134429542 : DO k = 1, npt
646 112025503 : sum = zero
647 336107763 : DO j = 1, n
648 336107763 : sum = sum + (xpt(k, j) - xopt(j))**2
649 : END DO
650 134429542 : IF (sum > distsq) THEN
651 30625552 : knew = k
652 30625552 : distsq = sum
653 : END IF
654 : END DO
655 : !
656 : ! If KNEW is positive, then set DSTEP, and branch back for the next
657 : ! iteration, which will generate a "model step".
658 : !
659 22404039 : IF (knew > 0) THEN
660 18966569 : dstep = MAX(MIN(tenth*SQRT(distsq), half*delta), rho)
661 18966569 : dsq = dstep*dstep
662 18966569 : CALL generate_next_model_step
663 18966569 : RETURN
664 : END IF
665 3437470 : IF (ratio <= zero .AND. MAX(delta, dnorm) <= rho) THEN
666 : !
667 : ! The calculations with the current value of RHO are complete. Pick
668 : ! next values of RHO and DELTA.
669 : !
670 2995176 : CALL update_rho()
671 2995176 : IF (do_return) RETURN
672 : END IF
673 : END IF
674 :
675 45823474 : skip_check = .FALSE.
676 : inner: DO
677 45823474 : knew = 0
678 45823474 : CALL trsapp(n, npt, xopt, xpt, gq, hq, pq, delta, d, w, w(np), w(np + n), w(np + 2*n), crvmin)
679 45823474 : dsq = zero
680 137474763 : DO i = 1, n
681 137474763 : dsq = dsq + d(i)**2
682 : END DO
683 45823474 : dnorm = MIN(delta, SQRT(dsq))
684 45823474 : IF (dnorm < half*rho) THEN
685 11160539 : knew = -1
686 11160539 : delta = tenth*delta
687 11160539 : ratio = -1.0_dp
688 11160539 : IF (delta <= 1.5_dp*rho) delta = rho
689 11160539 : IF (nf <= nfsav + 2) CYCLE outer
690 3091409 : temp = 0.125_dp*crvmin*rho*rho
691 3091409 : IF (temp <= MAX(diffa, diffb, diffc)) CYCLE outer
692 1593909 : CALL update_rho()
693 1593909 : IF (do_return) RETURN
694 : CYCLE inner
695 : END IF
696 : EXIT inner
697 : END DO inner
698 9566630 : EXIT outer
699 : END DO outer
700 : !
701 : ! Shift XBASE if XOPT may be too far from XBASE. First make the chan
702 : ! to BMAT that do not depend on ZMAT.
703 34662935 : CALL generate_next_model_step
704 :
705 : CONTAINS
706 : ! **************************************************************************************************
707 : !> \brief ...
708 : ! **************************************************************************************************
709 4589085 : SUBROUTINE update_rho()
710 4589085 : IF (rho > rhoend) THEN
711 3824276 : delta = half*rho
712 3824276 : ratio = rho/rhoend
713 3824276 : IF (ratio <= 16.0_dp) THEN
714 764803 : rho = rhoend
715 3059473 : ELSE IF (ratio <= 250.0_dp) THEN
716 764803 : rho = SQRT(ratio)*rhoend
717 : ELSE
718 2294670 : rho = tenth*rho
719 : END IF
720 3824276 : delta = MAX(delta, rho)
721 3824276 : nfsav = nf
722 : ELSE
723 : !
724 : ! Return from the calculation, after another Newton-Raphson step, if
725 : ! it is too short to have been tried before.
726 : !
727 764809 : IF (knew == -1) THEN
728 638430 : CALL calc_next_value_of_objective_func
729 : ELSE
730 126379 : opt%state = 7
731 126379 : CALL get_state
732 :
733 126379 : CALL finalize_optimization_after_failure
734 : END IF
735 764809 : do_return = .TRUE.
736 : END IF
737 4589085 : END SUBROUTINE update_rho
738 : ! **************************************************************************************************
739 : !> \brief ...
740 : ! **************************************************************************************************
741 3824301 : SUBROUTINE begin_initialization_of_interpolation()
742 3824301 : nfm = nf
743 3824301 : nfmm = nf - n
744 3824301 : nf = nf + 1
745 3824301 : IF (nfm <= 2*n) THEN
746 3824301 : IF (nfm >= 1 .AND. nfm <= N) THEN
747 1529775 : xpt(nf, nfm) = rhobeg
748 2294526 : ELSE IF (nfm > n) THEN
749 1529701 : xpt(nf, nfmm) = -rhobeg
750 : END IF
751 : ELSE
752 0 : itemp = (nfmm - 1)/n
753 0 : jpt = nfm - itemp*n - n
754 0 : ipt = jpt + itemp
755 0 : IF (ipt > n) THEN
756 0 : itemp = jpt
757 0 : jpt = ipt - n
758 0 : ipt = itemp
759 : END IF
760 0 : xipt = rhobeg
761 0 : IF (fval(ipt + np) < fval(ipt + 1)) xipt = -xipt
762 0 : XJPT = RHOBEG
763 0 : IF (fval(jpt + np) < fval(jpt + 1)) xjpt = -xjpt
764 0 : xpt(nf, ipt) = xipt
765 0 : xpt(nf, jpt) = xjpt
766 : END IF
767 : !
768 : ! Calculate the next value of F, label 70 being reached immediately
769 : ! after this calculation. The least function value so far and its in
770 : ! are required.
771 : !
772 11478579 : DO j = 1, n
773 11478579 : x(j) = xpt(nf, j) + xbase(j)
774 : END DO
775 3824301 : CALL check_number_of_steps
776 3824301 : END SUBROUTINE begin_initialization_of_interpolation
777 : ! **************************************************************************************************
778 : !> \brief ...
779 : ! **************************************************************************************************
780 53629504 : SUBROUTINE generate_next_model_step()
781 53629504 : IF (dsq <= 1.0e-3_dp*xoptsq) THEN
782 3285953 : tempq = 0.25_dp*xoptsq
783 19716006 : DO k = 1, npt
784 16430053 : sum = zero
785 49291775 : DO i = 1, n
786 49291775 : sum = sum + xpt(k, i)*xopt(i)
787 : END DO
788 16430053 : temp = pq(k)*sum
789 16430053 : sum = sum - half*xoptsq
790 16430053 : w(npt + k) = sum
791 52577728 : DO i = 1, n
792 32861722 : gq(i) = gq(i) + temp*xpt(k, i)
793 32861722 : xpt(k, i) = xpt(k, i) - half*xopt(i)
794 32861722 : vlag(i) = bmat(k, i)
795 32861722 : w(i) = sum*xpt(k, i) + tempq*xopt(i)
796 32861722 : ip = npt + i
797 98588558 : DO j = 1, i
798 82158505 : bmat(ip, j) = bmat(ip, j) + vlag(i)*w(j) + w(i)*vlag(j)
799 : END DO
800 : END DO
801 : END DO
802 : !
803 : ! Then the revisions of BMAT that depend on ZMAT are calculated.
804 : !
805 9858003 : DO k = 1, nptm
806 6572050 : sumz = zero
807 39433772 : DO i = 1, npt
808 32861722 : sumz = sumz + zmat(i, k)
809 39433772 : w(i) = w(npt + i)*zmat(i, k)
810 : END DO
811 19716886 : DO j = 1, n
812 13144836 : sum = tempq*sumz*xopt(j)
813 78876680 : DO i = 1, npt
814 65731844 : sum = sum + w(i)*xpt(i, j)
815 65731844 : vlag(j) = sum
816 78876680 : IF (k < idz) sum = -sum
817 : END DO
818 85448730 : DO i = 1, npt
819 78876680 : bmat(i, j) = bmat(i, j) + sum*zmat(i, k)
820 : END DO
821 : END DO
822 23002839 : DO i = 1, n
823 13144836 : ip = i + npt
824 13144836 : temp = vlag(i)
825 13144836 : IF (k < idz) temp = -temp
826 39436056 : DO j = 1, i
827 32864006 : bmat(ip, j) = bmat(ip, j) + temp*vlag(j)
828 : END DO
829 : END DO
830 : END DO
831 : !
832 : ! The following instructions complete the shift of XBASE, including
833 : ! the changes to the parameters of the quadratic model.
834 : !
835 3285953 : ih = 0
836 9858003 : DO j = 1, n
837 6572050 : w(j) = zero
838 39433772 : DO k = 1, npt
839 32861722 : w(j) = w(j) + pq(k)*xpt(k, j)
840 39433772 : xpt(k, j) = xpt(k, j) - half*xopt(j)
841 : END DO
842 19716446 : DO i = 1, j
843 9858443 : ih = ih + 1
844 9858443 : IF (i < j) gq(j) = gq(j) + hq(ih)*xopt(i)
845 9858443 : gq(i) = gq(i) + hq(ih)*xopt(j)
846 9858443 : hq(ih) = hq(ih) + w(i)*xopt(j) + xopt(i)*w(j)
847 16430493 : bmat(npt + i, j) = bmat(npt + j, i)
848 : END DO
849 : END DO
850 9858003 : DO j = 1, n
851 6572050 : xbase(j) = xbase(j) + xopt(j)
852 9858003 : xopt(j) = zero
853 : END DO
854 3285953 : xoptsq = zero
855 : END IF
856 : !
857 : ! Pick the model step if KNEW is positive. A different choice of D
858 : ! may be made later, if the choice of D by BIGLAG causes substantial
859 : ! cancellation in DENOM.
860 : !
861 53629504 : IF (knew > 0) THEN
862 : CALL biglag(n, npt, xopt, xpt, bmat, zmat, idz, ndim, knew, dstep, &
863 18966569 : d, alpha, vlag, vlag(npt + 1), w, w(np), w(np + n))
864 : END IF
865 : !
866 : ! Calculate VLAG and BETA for the current choice of D. The first NPT
867 : ! components of W_check will be held in W.
868 : !
869 321788740 : DO k = 1, npt
870 268159236 : suma = zero
871 268159236 : sumb = zero
872 268159236 : sum = zero
873 804545554 : DO j = 1, n
874 536386318 : suma = suma + xpt(k, j)*d(j)
875 536386318 : sumb = sumb + xpt(k, j)*xopt(j)
876 804545554 : sum = sum + bmat(k, j)*d(j)
877 : END DO
878 268159236 : w(k) = suma*(half*suma + sumb)
879 321788740 : vlag(k) = sum
880 : END DO
881 53629504 : beta = zero
882 160894370 : DO k = 1, nptm
883 107264866 : sum = zero
884 643651184 : DO i = 1, npt
885 643651184 : sum = sum + zmat(i, k)*w(i)
886 : END DO
887 107264866 : IF (k < idz) THEN
888 0 : beta = beta + sum*sum
889 0 : sum = -sum
890 : ELSE
891 107264866 : beta = beta - sum*sum
892 : END IF
893 697280688 : DO i = 1, npt
894 643651184 : vlag(i) = vlag(i) + sum*zmat(i, k)
895 : END DO
896 : END DO
897 53629504 : bsum = zero
898 53629504 : dx = zero
899 160894370 : DO j = 1, n
900 107264866 : sum = zero
901 643651184 : DO i = 1, npt
902 643651184 : sum = sum + w(i)*bmat(i, j)
903 : END DO
904 107264866 : bsum = bsum + sum*d(j)
905 107264866 : jp = npt + j
906 321825592 : DO k = 1, n
907 321825592 : sum = sum + bmat(jp, k)*d(k)
908 : END DO
909 107264866 : vlag(jp) = sum
910 107264866 : bsum = bsum + sum*d(j)
911 160894370 : dx = dx + d(j)*xopt(j)
912 : END DO
913 53629504 : beta = dx*dx + dsq*(xoptsq + dx + dx + half*dsq) + beta - bsum
914 53629504 : vlag(kopt) = vlag(kopt) + one
915 : !
916 : ! If KNEW is positive and if the cancellation in DENOM is unacceptab
917 : ! then BIGDEN calculates an alternative model step, XNEW being used
918 : ! working space.
919 : !
920 53629504 : IF (knew > 0) THEN
921 18966569 : temp = one + alpha*beta/vlag(knew)**2
922 18966569 : IF (ABS(temp) <= 0.8_dp) THEN
923 : CALL bigden(n, npt, xopt, xpt, bmat, zmat, idz, ndim, kopt, &
924 110 : knew, d, w, vlag, beta, xnew, w(ndim + 1), w(6*ndim + 1))
925 : END IF
926 : END IF
927 : !
928 : ! Calculate the next value of the objective function.
929 : !
930 53629504 : CALL calc_next_value_of_objective_func
931 53629504 : END SUBROUTINE generate_next_model_step
932 : ! **************************************************************************************************
933 : !> \brief ...
934 : ! **************************************************************************************************
935 54267934 : SUBROUTINE calc_next_value_of_objective_func()
936 162809698 : DO i = 1, n
937 108541764 : xnew(i) = xopt(i) + d(i)
938 162809698 : x(i) = xbase(i) + xnew(i)
939 : END DO
940 54267934 : nf = nf + 1
941 54267934 : CALL check_number_of_steps
942 54267934 : END SUBROUTINE calc_next_value_of_objective_func
943 : ! **************************************************************************************************
944 : !> \brief ...
945 : ! **************************************************************************************************
946 58092235 : SUBROUTINE check_number_of_steps()
947 58092235 : IF (nf > nftest) THEN
948 : ! return to many steps
949 11 : nf = nf - 1
950 11 : opt%state = 3
951 11 : CALL get_state
952 11 : CALL finalize_optimization_after_failure
953 : ELSE
954 :
955 58092224 : CALL get_state
956 :
957 58092224 : opt%state = 2
958 : END IF
959 58092235 : END SUBROUTINE check_number_of_steps
960 : ! **************************************************************************************************
961 : !> \brief ...
962 : ! **************************************************************************************************
963 764824 : SUBROUTINE finalize_optimization_after_failure()
964 764824 : IF (fopt <= f) THEN
965 714084 : DO i = 1, n
966 714084 : x(i) = xbase(i) + xopt(i)
967 : END DO
968 237921 : f = fopt
969 : END IF
970 :
971 764824 : CALL get_state
972 764824 : END SUBROUTINE finalize_optimization_after_failure
973 : ! **************************************************************************************************
974 : !> \brief ...
975 : ! **************************************************************************************************
976 59621872 : SUBROUTINE get_state()
977 59621872 : opt%np = np
978 59621872 : opt%nh = nh
979 59621872 : opt%nptm = nptm
980 59621872 : opt%nftest = nftest
981 59621872 : opt%idz = idz
982 59621872 : opt%itest = itest
983 59621872 : opt%nf = nf
984 59621872 : opt%nfm = nfm
985 59621872 : opt%nfmm = nfmm
986 59621872 : opt%nfsav = nfsav
987 59621872 : opt%knew = knew
988 59621872 : opt%kopt = kopt
989 59621872 : opt%ksave = ksave
990 59621872 : opt%ktemp = ktemp
991 59621872 : opt%rhosq = rhosq
992 59621872 : opt%recip = recip
993 59621872 : opt%reciq = reciq
994 59621872 : opt%fbeg = fbeg
995 59621872 : opt%fopt = fopt
996 59621872 : opt%diffa = diffa
997 59621872 : opt%xoptsq = xoptsq
998 59621872 : opt%rho = rho
999 59621872 : opt%delta = delta
1000 59621872 : opt%dsq = dsq
1001 59621872 : opt%dnorm = dnorm
1002 59621872 : opt%ratio = ratio
1003 59621872 : opt%temp = temp
1004 59621872 : opt%tempq = tempq
1005 59621872 : opt%beta = beta
1006 59621872 : opt%dx = dx
1007 59621872 : opt%vquad = vquad
1008 59621872 : opt%diff = diff
1009 59621872 : opt%diffc = diffc
1010 59621872 : opt%diffb = diffb
1011 59621872 : opt%fsave = fsave
1012 59621872 : opt%detrat = detrat
1013 59621872 : opt%hdiag = hdiag
1014 59621872 : opt%distsq = distsq
1015 59621872 : opt%gisq = gisq
1016 59621872 : opt%gqsq = gqsq
1017 59621872 : opt%f = f
1018 59621872 : opt%bstep = bstep
1019 59621872 : opt%alpha = alpha
1020 59621872 : opt%dstep = dstep
1021 59621872 : END SUBROUTINE get_state
1022 : ! **************************************************************************************************
1023 : !> \brief ...
1024 : ! **************************************************************************************************
1025 58092223 : SUBROUTINE set_state()
1026 58092223 : np = opt%np
1027 58092223 : nh = opt%nh
1028 58092223 : nptm = opt%nptm
1029 58092223 : nftest = opt%nftest
1030 58092223 : idz = opt%idz
1031 58092223 : itest = opt%itest
1032 58092223 : nf = opt%nf
1033 58092223 : nfm = opt%nfm
1034 58092223 : nfmm = opt%nfmm
1035 58092223 : nfsav = opt%nfsav
1036 58092223 : knew = opt%knew
1037 58092223 : kopt = opt%kopt
1038 58092223 : ksave = opt%ksave
1039 58092223 : ktemp = opt%ktemp
1040 58092223 : rhosq = opt%rhosq
1041 58092223 : recip = opt%recip
1042 58092223 : reciq = opt%reciq
1043 58092223 : fbeg = opt%fbeg
1044 58092223 : fopt = opt%fopt
1045 58092223 : diffa = opt%diffa
1046 58092223 : xoptsq = opt%xoptsq
1047 58092223 : rho = opt%rho
1048 58092223 : delta = opt%delta
1049 58092223 : dsq = opt%dsq
1050 58092223 : dnorm = opt%dnorm
1051 58092223 : ratio = opt%ratio
1052 58092223 : temp = opt%temp
1053 58092223 : tempq = opt%tempq
1054 58092223 : beta = opt%beta
1055 58092223 : dx = opt%dx
1056 58092223 : vquad = opt%vquad
1057 58092223 : diff = opt%diff
1058 58092223 : diffc = opt%diffc
1059 58092223 : diffb = opt%diffb
1060 58092223 : fsave = opt%fsave
1061 58092223 : detrat = opt%detrat
1062 58092223 : hdiag = opt%hdiag
1063 58092223 : distsq = opt%distsq
1064 58092223 : gisq = opt%gisq
1065 58092223 : gqsq = opt%gqsq
1066 58092223 : f = opt%f
1067 58092223 : bstep = opt%bstep
1068 58092223 : alpha = opt%alpha
1069 58092223 : dstep = opt%dstep
1070 58092223 : END SUBROUTINE set_state
1071 :
1072 : END SUBROUTINE newuob
1073 :
1074 : ! **************************************************************************************************
1075 : !> \brief ...
1076 : !> \param n ...
1077 : !> \param npt ...
1078 : !> \param xopt ...
1079 : !> \param xpt ...
1080 : !> \param bmat ...
1081 : !> \param zmat ...
1082 : !> \param idz ...
1083 : !> \param ndim ...
1084 : !> \param kopt ...
1085 : !> \param knew ...
1086 : !> \param d ...
1087 : !> \param w ...
1088 : !> \param vlag ...
1089 : !> \param beta ...
1090 : !> \param s ...
1091 : !> \param wvec ...
1092 : !> \param prod ...
1093 : ! **************************************************************************************************
1094 110 : SUBROUTINE bigden(n, npt, xopt, xpt, bmat, zmat, idz, ndim, kopt, &
1095 110 : knew, d, w, vlag, beta, s, wvec, prod)
1096 :
1097 : INTEGER, INTENT(in) :: n, npt
1098 : REAL(dp), DIMENSION(*), INTENT(in) :: xopt
1099 : REAL(dp), DIMENSION(npt, *), INTENT(in) :: xpt
1100 : INTEGER, INTENT(in) :: ndim, idz
1101 : REAL(dp), DIMENSION(npt, *), INTENT(inout) :: zmat
1102 : REAL(dp), DIMENSION(ndim, *), INTENT(inout) :: bmat
1103 : INTEGER, INTENT(inout) :: kopt, knew
1104 : REAL(dp), DIMENSION(*), INTENT(inout) :: d, w, vlag
1105 : REAL(dp), INTENT(inout) :: beta
1106 : REAL(dp), DIMENSION(*), INTENT(inout) :: s
1107 : REAL(dp), DIMENSION(ndim, *), INTENT(inout) :: wvec, prod
1108 :
1109 : REAL(dp), PARAMETER :: half = 0.5_dp, one = 1._dp, &
1110 : quart = 0.25_dp, two = 2._dp, &
1111 : zero = 0._dp
1112 :
1113 : INTEGER :: i, ip, isave, iterc, iu, j, jc, k, ksav, &
1114 : nptm, nw
1115 : REAL(dp) :: alpha, angle, dd, denmax, denold, densav, diff, ds, dstemp, dtest, ss, ssden, &
1116 : sstemp, step, sum, sumold, tau, temp, tempa, tempb, tempc, xoptd, xopts, xoptsq
1117 : REAL(dp), DIMENSION(9) :: den, denex, par
1118 :
1119 : !
1120 : ! N is the number of variables.
1121 : ! NPT is the number of interpolation equations.
1122 : ! XOPT is the best interpolation point so far.
1123 : ! XPT contains the coordinates of the current interpolation points.
1124 : ! BMAT provides the last N columns of H.
1125 : ! ZMAT and IDZ give a factorization of the first NPT by NPT submatri
1126 : ! NDIM is the first dimension of BMAT and has the value NPT+N.
1127 : ! KOPT is the index of the optimal interpolation point.
1128 : ! KNEW is the index of the interpolation point that is going to be m
1129 : ! D will be set to the step from XOPT to the new point, and on entry
1130 : ! should be the D that was calculated by the last call of BIGLAG.
1131 : ! length of the initial D provides a trust region bound on the fin
1132 : ! W will be set to Wcheck for the final choice of D.
1133 : ! VLAG will be set to Theta*Wcheck+e_b for the final choice of D.
1134 : ! BETA will be set to the value that will occur in the updating form
1135 : ! when the KNEW-th interpolation point is moved to its new positio
1136 : ! S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be us
1137 : ! for working space.
1138 : !
1139 : ! D is calculated in a way that should provide a denominator with a
1140 : ! modulus in the updating formula when the KNEW-th interpolation poi
1141 : ! shifted to the new position XOPT+D.
1142 : !
1143 :
1144 110 : nptm = npt - n - 1
1145 : !
1146 : ! Store the first NPT elements of the KNEW-th column of H in W(N+1)
1147 : ! to W(N+NPT).
1148 : !
1149 660 : DO k = 1, npt
1150 660 : w(n + k) = zero
1151 : END DO
1152 330 : DO j = 1, nptm
1153 220 : temp = zmat(knew, j)
1154 220 : IF (j < idz) temp = -temp
1155 1430 : DO k = 1, npt
1156 1320 : w(n + k) = w(n + k) + temp*zmat(k, j)
1157 : END DO
1158 : END DO
1159 110 : alpha = w(n + knew)
1160 : !
1161 : ! The initial search direction D is taken from the last call of BIGL
1162 : ! and the initial S is set below, usually to the direction from X_OP
1163 : ! to X_KNEW, but a different direction to an interpolation point may
1164 : ! be chosen, in order to prevent S from being nearly parallel to D.
1165 : !
1166 110 : dd = zero
1167 110 : ds = zero
1168 110 : ss = zero
1169 110 : xoptsq = zero
1170 330 : DO i = 1, n
1171 220 : dd = dd + d(i)**2
1172 220 : s(i) = xpt(knew, i) - xopt(i)
1173 220 : ds = ds + d(i)*s(i)
1174 220 : ss = ss + s(i)**2
1175 330 : xoptsq = xoptsq + xopt(i)**2
1176 : END DO
1177 110 : IF (ds*ds > 0.99_dp*dd*ss) THEN
1178 0 : ksav = knew
1179 0 : dtest = ds*ds/ss
1180 0 : DO k = 1, npt
1181 0 : IF (k /= kopt) THEN
1182 : dstemp = zero
1183 : sstemp = zero
1184 0 : DO i = 1, n
1185 0 : diff = xpt(k, i) - xopt(i)
1186 0 : dstemp = dstemp + d(i)*diff
1187 0 : sstemp = sstemp + diff*diff
1188 : END DO
1189 0 : IF (dstemp*dstemp/sstemp < dtest) THEN
1190 0 : ksav = k
1191 0 : dtest = dstemp*dstemp/sstemp
1192 0 : ds = dstemp
1193 0 : ss = sstemp
1194 : END IF
1195 : END IF
1196 : END DO
1197 0 : DO i = 1, n
1198 0 : s(i) = xpt(ksav, i) - xopt(i)
1199 : END DO
1200 : END IF
1201 110 : ssden = dd*ss - ds*ds
1202 110 : iterc = 0
1203 110 : densav = zero
1204 : !
1205 : ! Begin the iteration by overwriting S with a vector that has the
1206 : ! required length and direction.
1207 : !
1208 : mainloop: DO
1209 176 : iterc = iterc + 1
1210 176 : temp = one/SQRT(ssden)
1211 176 : xoptd = zero
1212 176 : xopts = zero
1213 528 : DO i = 1, n
1214 352 : s(i) = temp*(dd*s(i) - ds*d(i))
1215 352 : xoptd = xoptd + xopt(i)*d(i)
1216 528 : xopts = xopts + xopt(i)*s(i)
1217 : END DO
1218 : !
1219 : ! Set the coefficients of the first two terms of BETA.
1220 : !
1221 176 : tempa = half*xoptd*xoptd
1222 176 : tempb = half*xopts*xopts
1223 176 : den(1) = dd*(xoptsq + half*dd) + tempa + tempb
1224 176 : den(2) = two*xoptd*dd
1225 176 : den(3) = two*xopts*dd
1226 176 : den(4) = tempa - tempb
1227 176 : den(5) = xoptd*xopts
1228 880 : DO i = 6, 9
1229 880 : den(i) = zero
1230 : END DO
1231 : !
1232 : ! Put the coefficients of Wcheck in WVEC.
1233 : !
1234 1056 : DO k = 1, npt
1235 : tempa = zero
1236 : tempb = zero
1237 : tempc = zero
1238 2640 : DO i = 1, n
1239 1760 : tempa = tempa + xpt(k, i)*d(i)
1240 1760 : tempb = tempb + xpt(k, i)*s(i)
1241 2640 : tempc = tempc + xpt(k, i)*xopt(i)
1242 : END DO
1243 880 : wvec(k, 1) = quart*(tempa*tempa + tempb*tempb)
1244 880 : wvec(k, 2) = tempa*tempc
1245 880 : wvec(k, 3) = tempb*tempc
1246 880 : wvec(k, 4) = quart*(tempa*tempa - tempb*tempb)
1247 1056 : wvec(k, 5) = half*tempa*tempb
1248 : END DO
1249 528 : DO i = 1, n
1250 352 : ip = i + npt
1251 352 : wvec(ip, 1) = zero
1252 352 : wvec(ip, 2) = d(i)
1253 352 : wvec(ip, 3) = s(i)
1254 352 : wvec(ip, 4) = zero
1255 528 : wvec(ip, 5) = zero
1256 : END DO
1257 : !
1258 : ! Put the coefficients of THETA*Wcheck in PROD.
1259 : !
1260 1056 : DO jc = 1, 5
1261 880 : nw = npt
1262 880 : IF (jc == 2 .OR. jc == 3) nw = ndim
1263 5280 : DO k = 1, npt
1264 5280 : prod(k, jc) = zero
1265 : END DO
1266 2640 : DO j = 1, nptm
1267 : sum = zero
1268 10560 : DO k = 1, npt
1269 10560 : sum = sum + zmat(k, j)*wvec(k, jc)
1270 : END DO
1271 1760 : IF (j < idz) sum = -sum
1272 11440 : DO k = 1, npt
1273 10560 : prod(k, jc) = prod(k, jc) + sum*zmat(k, j)
1274 : END DO
1275 : END DO
1276 880 : IF (nw == ndim) THEN
1277 2112 : DO k = 1, npt
1278 : sum = zero
1279 5280 : DO j = 1, n
1280 5280 : sum = sum + bmat(k, j)*wvec(npt + j, jc)
1281 : END DO
1282 2112 : prod(k, jc) = prod(k, jc) + sum
1283 : END DO
1284 : END IF
1285 2816 : DO j = 1, n
1286 : sum = zero
1287 11968 : DO i = 1, nw
1288 11968 : sum = sum + bmat(i, j)*wvec(i, jc)
1289 : END DO
1290 2640 : prod(npt + j, jc) = sum
1291 : END DO
1292 : END DO
1293 : !
1294 : ! Include in DEN the part of BETA that depends on THETA.
1295 : !
1296 1408 : DO k = 1, ndim
1297 : sum = zero
1298 7392 : DO I = 1, 5
1299 6160 : par(i) = half*prod(k, i)*wvec(k, i)
1300 7392 : sum = sum + par(i)
1301 : END DO
1302 1232 : den(1) = den(1) - par(1) - sum
1303 1232 : tempa = prod(k, 1)*wvec(k, 2) + prod(k, 2)*wvec(k, 1)
1304 1232 : tempb = prod(k, 2)*wvec(k, 4) + prod(k, 4)*wvec(k, 2)
1305 1232 : tempc = prod(k, 3)*wvec(k, 5) + prod(k, 5)*wvec(k, 3)
1306 1232 : den(2) = den(2) - tempa - half*(tempb + tempc)
1307 1232 : den(6) = den(6) - half*(tempb - tempc)
1308 1232 : tempa = prod(k, 1)*wvec(k, 3) + prod(k, 3)*wvec(k, 1)
1309 1232 : tempb = prod(k, 2)*wvec(k, 5) + prod(k, 5)*wvec(k, 2)
1310 1232 : tempc = prod(k, 3)*wvec(k, 4) + prod(k, 4)*wvec(k, 3)
1311 1232 : den(3) = den(3) - tempa - half*(tempb - tempc)
1312 1232 : den(7) = den(7) - half*(tempb + tempc)
1313 1232 : tempa = prod(k, 1)*wvec(k, 4) + prod(k, 4)*wvec(k, 1)
1314 1232 : den(4) = den(4) - tempa - par(2) + par(3)
1315 1232 : tempa = prod(k, 1)*wvec(k, 5) + prod(k, 5)*wvec(k, 1)
1316 1232 : tempb = prod(k, 2)*wvec(k, 3) + prod(k, 3)*wvec(k, 2)
1317 1232 : den(5) = den(5) - tempa - half*tempb
1318 1232 : den(8) = den(8) - par(4) + par(5)
1319 1232 : tempa = prod(k, 4)*wvec(k, 5) + prod(k, 5)*wvec(k, 4)
1320 1408 : den(9) = den(9) - half*tempa
1321 : END DO
1322 : !
1323 : ! Extend DEN so that it holds all the coefficients of DENOM.
1324 : !
1325 : sum = zero
1326 1056 : DO i = 1, 5
1327 880 : par(i) = half*prod(knew, i)**2
1328 1056 : sum = sum + par(i)
1329 : END DO
1330 176 : denex(1) = alpha*den(1) + par(1) + sum
1331 176 : tempa = two*prod(knew, 1)*prod(knew, 2)
1332 176 : tempb = prod(knew, 2)*prod(knew, 4)
1333 176 : tempc = prod(knew, 3)*prod(knew, 5)
1334 176 : denex(2) = alpha*den(2) + tempa + tempb + tempc
1335 176 : denex(6) = alpha*den(6) + tempb - tempc
1336 176 : tempa = two*prod(knew, 1)*prod(knew, 3)
1337 176 : tempb = prod(knew, 2)*prod(knew, 5)
1338 176 : tempc = prod(knew, 3)*prod(knew, 4)
1339 176 : denex(3) = alpha*den(3) + tempa + tempb - tempc
1340 176 : denex(7) = alpha*den(7) + tempb + tempc
1341 176 : tempa = two*prod(knew, 1)*prod(knew, 4)
1342 176 : denex(4) = alpha*den(4) + tempa + par(2) - par(3)
1343 176 : tempa = two*prod(knew, 1)*prod(knew, 5)
1344 176 : denex(5) = alpha*den(5) + tempa + prod(knew, 2)*prod(knew, 3)
1345 176 : denex(8) = alpha*den(8) + par(4) - par(5)
1346 176 : denex(9) = alpha*den(9) + prod(knew, 4)*prod(knew, 5)
1347 : !
1348 : ! Seek the value of the angle that maximizes the modulus of DENOM.
1349 : !
1350 176 : sum = denex(1) + denex(2) + denex(4) + denex(6) + denex(8)
1351 176 : denold = sum
1352 176 : denmax = sum
1353 176 : isave = 0
1354 176 : iu = 49
1355 176 : temp = twopi/REAL(IU + 1, dp)
1356 176 : par(1) = one
1357 8800 : DO i = 1, iu
1358 8624 : angle = REAL(i, dp)*temp
1359 8624 : par(2) = COS(angle)
1360 8624 : par(3) = SIN(angle)
1361 8624 : DO j = 4, 8, 2
1362 25872 : par(j) = par(2)*par(j - 2) - par(3)*par(j - 1)
1363 25872 : par(j + 1) = par(2)*par(j - 1) + par(3)*par(j - 2)
1364 : END DO
1365 : sumold = sum
1366 : sum = zero
1367 86240 : DO j = 1, 9
1368 86240 : sum = sum + denex(j)*par(j)
1369 : END DO
1370 8800 : IF (ABS(sum) > ABS(denmax)) THEN
1371 : denmax = sum
1372 : isave = i
1373 : tempa = sumold
1374 8020 : ELSE IF (i == isave + 1) THEN
1375 326 : tempb = sum
1376 : END IF
1377 : END DO
1378 176 : IF (isave == 0) tempa = sum
1379 86 : IF (isave == iu) tempb = denold
1380 176 : step = zero
1381 176 : IF (tempa /= tempb) THEN
1382 176 : tempa = tempa - denmax
1383 176 : tempb = tempb - denmax
1384 176 : step = half*(tempa - tempb)/(tempa + tempb)
1385 : END IF
1386 176 : angle = temp*(REAL(isave, dp) + step)
1387 : !
1388 : ! Calculate the new parameters of the denominator, the new VLAG vect
1389 : ! and the new D. Then test for convergence.
1390 : !
1391 176 : par(2) = COS(angle)
1392 176 : par(3) = SIN(angle)
1393 176 : DO j = 4, 8, 2
1394 528 : par(j) = par(2)*par(j - 2) - par(3)*par(j - 1)
1395 528 : par(j + 1) = par(2)*par(j - 1) + par(3)*par(j - 2)
1396 : END DO
1397 176 : beta = zero
1398 176 : denmax = zero
1399 1760 : DO j = 1, 9
1400 1584 : beta = beta + den(j)*par(j)
1401 1760 : denmax = denmax + denex(j)*par(j)
1402 : END DO
1403 1408 : DO k = 1, ndim
1404 1232 : vlag(k) = zero
1405 7568 : DO j = 1, 5
1406 7392 : vlag(k) = vlag(k) + prod(k, j)*par(j)
1407 : END DO
1408 : END DO
1409 176 : tau = vlag(knew)
1410 176 : dd = zero
1411 176 : tempa = zero
1412 176 : tempb = zero
1413 528 : DO i = 1, n
1414 352 : d(i) = par(2)*d(i) + par(3)*s(i)
1415 352 : w(i) = xopt(i) + d(i)
1416 352 : dd = dd + d(i)**2
1417 352 : tempa = tempa + d(i)*w(i)
1418 528 : tempb = tempb + w(i)*w(i)
1419 : END DO
1420 176 : IF (iterc >= n) EXIT mainloop
1421 110 : IF (iterc >= 1) densav = MAX(densav, denold)
1422 110 : IF (ABS(denmax) <= 1.1_dp*ABS(densav)) EXIT mainloop
1423 198 : densav = denmax
1424 : !
1425 : ! Set S to half the gradient of the denominator with respect to D.
1426 : ! Then branch for the next iteration.
1427 : !
1428 198 : DO i = 1, n
1429 132 : temp = tempa*xopt(i) + tempb*d(i) - vlag(npt + i)
1430 198 : s(i) = tau*bmat(knew, i) + alpha*temp
1431 : END DO
1432 396 : DO k = 1, npt
1433 : sum = zero
1434 990 : DO j = 1, n
1435 990 : sum = sum + xpt(k, j)*w(j)
1436 : END DO
1437 330 : temp = (tau*w(n + k) - alpha*vlag(k))*sum
1438 1056 : DO i = 1, n
1439 990 : s(i) = s(i) + temp*xpt(k, i)
1440 : END DO
1441 : END DO
1442 : ss = zero
1443 : ds = zero
1444 198 : DO i = 1, n
1445 132 : ss = ss + s(i)**2
1446 198 : ds = ds + d(i)*s(i)
1447 : END DO
1448 66 : ssden = dd*ss - ds*ds
1449 176 : IF (ssden < 1.0e-8_dp*dd*ss) EXIT mainloop
1450 : END DO mainloop
1451 : !
1452 : ! Set the vector W before the RETURN from the subroutine.
1453 : !
1454 880 : DO k = 1, ndim
1455 770 : w(k) = zero
1456 4730 : DO j = 1, 5
1457 4620 : w(k) = w(k) + wvec(k, j)*par(j)
1458 : END DO
1459 : END DO
1460 110 : vlag(kopt) = vlag(kopt) + one
1461 :
1462 110 : END SUBROUTINE bigden
1463 :
1464 : ! **************************************************************************************************
1465 : !> \brief ...
1466 : !> \param n ...
1467 : !> \param npt ...
1468 : !> \param xopt ...
1469 : !> \param xpt ...
1470 : !> \param bmat ...
1471 : !> \param zmat ...
1472 : !> \param idz ...
1473 : !> \param ndim ...
1474 : !> \param knew ...
1475 : !> \param delta ...
1476 : !> \param d ...
1477 : !> \param alpha ...
1478 : !> \param hcol ...
1479 : !> \param gc ...
1480 : !> \param gd ...
1481 : !> \param s ...
1482 : !> \param w ...
1483 : ! **************************************************************************************************
1484 18966569 : SUBROUTINE biglag(n, npt, xopt, xpt, bmat, zmat, idz, ndim, knew, &
1485 : delta, d, alpha, hcol, gc, gd, s, w)
1486 : INTEGER, INTENT(in) :: n, npt
1487 : REAL(dp), DIMENSION(*), INTENT(in) :: xopt
1488 : REAL(dp), DIMENSION(npt, *), INTENT(in) :: xpt
1489 : INTEGER, INTENT(in) :: ndim, idz
1490 : REAL(dp), DIMENSION(npt, *), INTENT(inout) :: zmat
1491 : REAL(dp), DIMENSION(ndim, *), INTENT(inout) :: bmat
1492 : INTEGER, INTENT(inout) :: knew
1493 : REAL(dp), INTENT(inout) :: delta
1494 : REAL(dp), DIMENSION(*), INTENT(inout) :: d
1495 : REAL(dp), INTENT(inout) :: alpha
1496 : REAL(dp), DIMENSION(*), INTENT(inout) :: hcol, gc, gd, s, w
1497 :
1498 : REAL(dp), PARAMETER :: half = 0.5_dp, one = 1._dp, zero = 0._dp
1499 :
1500 : INTEGER :: i, isave, iterc, iu, j, k, nptm
1501 : REAL(dp) :: angle, cf1, cf2, cf3, cf4, cf5, cth, dd, &
1502 : delsq, denom, dhd, gg, scale, sp, ss, &
1503 : step, sth, sum, tau, taubeg, taumax, &
1504 : tauold, temp, tempa, tempb
1505 :
1506 : !
1507 : !
1508 : ! N is the number of variables.
1509 : ! NPT is the number of interpolation equations.
1510 : ! XOPT is the best interpolation point so far.
1511 : ! XPT contains the coordinates of the current interpolation points.
1512 : ! BMAT provides the last N columns of H.
1513 : ! ZMAT and IDZ give a factorization of the first NPT by NPT submatrix
1514 : ! NDIM is the first dimension of BMAT and has the value NPT+N.
1515 : ! KNEW is the index of the interpolation point that is going to be m
1516 : ! DELTA is the current trust region bound.
1517 : ! D will be set to the step from XOPT to the new point.
1518 : ! ALPHA will be set to the KNEW-th diagonal element of the H matrix.
1519 : ! HCOL, GC, GD, S and W will be used for working space.
1520 : !
1521 : ! The step D is calculated in a way that attempts to maximize the mo
1522 : ! of LFUNC(XOPT+D), subject to the bound ||D|| <= DELTA, where LFU
1523 : ! the KNEW-th Lagrange function.
1524 : !
1525 :
1526 18966569 : delsq = delta*delta
1527 18966569 : nptm = npt - n - 1
1528 : !
1529 : ! Set the first NPT components of HCOL to the leading elements of th
1530 : ! KNEW-th column of H.
1531 : !
1532 18966569 : iterc = 0
1533 113804106 : DO k = 1, npt
1534 113804106 : hcol(k) = zero
1535 : END DO
1536 56902053 : DO j = 1, nptm
1537 37935484 : temp = zmat(knew, j)
1538 37935484 : IF (j < idz) temp = -temp
1539 246604529 : DO k = 1, npt
1540 227637960 : hcol(k) = hcol(k) + temp*zmat(k, j)
1541 : END DO
1542 : END DO
1543 18966569 : alpha = hcol(knew)
1544 : !
1545 : ! Set the unscaled initial direction D. Form the gradient of LFUNC a
1546 : ! XOPT, and multiply D by the second derivative matrix of LFUNC.
1547 : !
1548 18966569 : dd = zero
1549 56902053 : DO i = 1, n
1550 37935484 : d(i) = xpt(knew, i) - xopt(i)
1551 37935484 : gc(i) = bmat(knew, i)
1552 37935484 : gd(i) = zero
1553 56902053 : dd = dd + d(i)**2
1554 : END DO
1555 113804106 : DO k = 1, npt
1556 : temp = zero
1557 : sum = zero
1558 284540013 : DO j = 1, n
1559 189702476 : temp = temp + xpt(k, j)*xopt(j)
1560 284540013 : sum = sum + xpt(k, j)*d(j)
1561 : END DO
1562 94837537 : temp = hcol(k)*temp
1563 94837537 : sum = hcol(k)*sum
1564 303506582 : DO i = 1, n
1565 189702476 : gc(i) = gc(i) + temp*xpt(k, i)
1566 284540013 : gd(i) = gd(i) + sum*xpt(k, i)
1567 : END DO
1568 : END DO
1569 : !
1570 : ! Scale D and GD, with a sign change if required. Set S to another
1571 : ! vector in the initial two dimensional subspace.
1572 : !
1573 : gg = zero
1574 : sp = zero
1575 : dhd = zero
1576 56902053 : DO i = 1, n
1577 37935484 : gg = gg + gc(i)**2
1578 37935484 : sp = sp + d(i)*gc(i)
1579 56902053 : dhd = dhd + d(i)*gd(i)
1580 : END DO
1581 18966569 : scale = delta/SQRT(dd)
1582 18966569 : IF (sp*dhd < zero) scale = -scale
1583 18966569 : temp = zero
1584 18966569 : IF (sp*sp > 0.99_dp*dd*gg) temp = one
1585 18966569 : tau = scale*(ABS(sp) + half*scale*ABS(dhd))
1586 18966569 : IF (gg*delsq < 0.01_dp*tau*tau) temp = one
1587 56902053 : DO i = 1, n
1588 37935484 : d(i) = scale*d(i)
1589 37935484 : gd(i) = scale*gd(i)
1590 56902053 : s(i) = gc(i) + temp*gd(i)
1591 : END DO
1592 : !
1593 : ! Begin the iteration by overwriting S with a vector that has the
1594 : ! required length and direction, except that termination occurs if
1595 : ! the given D and S are nearly parallel.
1596 : !
1597 : mainloop: DO
1598 26854239 : iterc = iterc + 1
1599 26854239 : dd = zero
1600 26854239 : sp = zero
1601 26854239 : ss = zero
1602 80566399 : DO i = 1, n
1603 53712160 : dd = dd + d(i)**2
1604 53712160 : sp = sp + d(i)*s(i)
1605 80566399 : ss = ss + s(i)**2
1606 : END DO
1607 26854239 : temp = dd*ss - sp*sp
1608 26854239 : IF (temp <= 1.0e-8_dp*dd*ss) EXIT mainloop
1609 25686236 : denom = SQRT(temp)
1610 77062116 : DO i = 1, n
1611 51375880 : s(i) = (dd*s(i) - sp*d(i))/denom
1612 77062116 : w(i) = zero
1613 : END DO
1614 : !
1615 : ! Calculate the coefficients of the objective function on the circle
1616 : ! beginning with the multiplication of S by the second derivative ma
1617 : !
1618 154124232 : DO k = 1, npt
1619 : sum = zero
1620 385353776 : DO j = 1, n
1621 385353776 : sum = sum + xpt(k, j)*s(j)
1622 : END DO
1623 128437996 : sum = hcol(k)*sum
1624 411040012 : DO i = 1, n
1625 385353776 : w(i) = w(i) + sum*xpt(k, i)
1626 : END DO
1627 : END DO
1628 : cf1 = zero
1629 : cf2 = zero
1630 : cf3 = zero
1631 : cf4 = zero
1632 : cf5 = zero
1633 77062116 : DO i = 1, n
1634 51375880 : cf1 = cf1 + s(i)*w(i)
1635 51375880 : cf2 = cf2 + d(i)*gc(i)
1636 51375880 : cf3 = cf3 + s(i)*gc(i)
1637 51375880 : cf4 = cf4 + d(i)*gd(i)
1638 77062116 : cf5 = cf5 + s(i)*gd(i)
1639 : END DO
1640 25686236 : cf1 = half*cf1
1641 25686236 : cf4 = half*cf4 - cf1
1642 : !
1643 : ! Seek the value of the angle that maximizes the modulus of TAU.
1644 : !
1645 25686236 : taubeg = cf1 + cf2 + cf4
1646 25686236 : taumax = taubeg
1647 25686236 : tauold = taubeg
1648 25686236 : isave = 0
1649 25686236 : iu = 49
1650 25686236 : temp = twopi/REAL(iu + 1, DP)
1651 1284311800 : DO i = 1, iu
1652 1258625564 : angle = REAL(i, dp)*temp
1653 1258625564 : cth = COS(angle)
1654 1258625564 : sth = SIN(angle)
1655 1258625564 : tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1656 1258625564 : IF (ABS(tau) > ABS(taumax)) THEN
1657 : taumax = tau
1658 : isave = i
1659 : tempa = tauold
1660 1192879501 : ELSE IF (i == isave + 1) THEN
1661 27926669 : tempb = taU
1662 : END IF
1663 1284311800 : tauold = tau
1664 : END DO
1665 25686236 : IF (isave == 0) tempa = tau
1666 15085913 : IF (isave == iu) tempb = taubeg
1667 25686236 : step = zero
1668 25686236 : IF (tempa /= tempb) THEN
1669 25686236 : tempa = tempa - taumax
1670 25686236 : tempb = tempb - taumax
1671 25686236 : step = half*(tempa - tempb)/(tempa + tempb)
1672 : END IF
1673 25686236 : angle = temp*(REAL(isave, DP) + step)
1674 : !
1675 : ! Calculate the new D and GD. Then test for convergence.
1676 : !
1677 25686236 : cth = COS(angle)
1678 25686236 : sth = SIN(angle)
1679 25686236 : tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1680 77062116 : DO i = 1, n
1681 51375880 : d(i) = cth*d(i) + sth*s(i)
1682 51375880 : gd(i) = cth*gd(i) + sth*w(i)
1683 77062116 : s(i) = gc(i) + gd(i)
1684 : END DO
1685 25686236 : IF (ABS(tau) <= 1.1_dp*ABS(taubeg)) EXIT mainloop
1686 26854239 : IF (iterc >= n) EXIT mainloop
1687 : END DO mainloop
1688 :
1689 18966569 : END SUBROUTINE biglag
1690 :
1691 : ! **************************************************************************************************
1692 : !> \brief ...
1693 : !> \param n ...
1694 : !> \param npt ...
1695 : !> \param xopt ...
1696 : !> \param xpt ...
1697 : !> \param gq ...
1698 : !> \param hq ...
1699 : !> \param pq ...
1700 : !> \param delta ...
1701 : !> \param step ...
1702 : !> \param d ...
1703 : !> \param g ...
1704 : !> \param hd ...
1705 : !> \param hs ...
1706 : !> \param crvmin ...
1707 : ! **************************************************************************************************
1708 45823474 : SUBROUTINE trsapp(n, npt, xopt, xpt, gq, hq, pq, delta, step, d, g, hd, hs, crvmin)
1709 :
1710 : INTEGER, INTENT(IN) :: n, npt
1711 : REAL(dp), DIMENSION(*), INTENT(IN) :: xopt
1712 : REAL(dp), DIMENSION(npt, *), &
1713 : INTENT(IN) :: xpt
1714 : REAL(dp), DIMENSION(*), INTENT(INOUT) :: gq, hq, pq
1715 : REAL(dp), INTENT(IN) :: delta
1716 : REAL(dp), DIMENSION(*), INTENT(INOUT) :: step, d, g, hd, hs
1717 : REAL(dp), INTENT(INOUT) :: crvmin
1718 :
1719 : REAL(dp), PARAMETER :: half = 0.5_dp, zero = 0.0_dp
1720 :
1721 : INTEGER :: i, isave, iterc, itermax, &
1722 : itersw, iu, j
1723 : LOGICAL :: jump1, jump2
1724 : REAL(dp) :: alpha, angle, angtest, bstep, cf, cth, dd, delsq, dg, dhd, &
1725 : dhs, ds, gg, ggbeg, ggsav, qadd, qbeg, qmin, qnew, qred, qsav, ratio, &
1726 : reduc, sg, sgk, shs, ss, sth, temp, tempa, tempb
1727 :
1728 : !
1729 : ! N is the number of variables of a quadratic objective function, Q
1730 : ! The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meani
1731 : ! in order to define the current quadratic model Q.
1732 : ! DELTA is the trust region radius, and has to be positive.
1733 : ! STEP will be set to the calculated trial step.
1734 : ! The arrays D, G, HD and HS will be used for working space.
1735 : ! CRVMIN will be set to the least curvature of H along the conjugate
1736 : ! directions that occur, except that it is set to zero if STEP goe
1737 : ! all the way to the trust region boundary.
1738 : !
1739 : ! The calculation of STEP begins with the truncated conjugate gradient
1740 : ! method. If the boundary of the trust region is reached, then further
1741 : ! changes to STEP may be made, each one being in the 2D space spanned
1742 : ! by the current STEP and the corresponding gradient of Q. Thus STEP
1743 : ! should provide a substantial reduction to Q within the trust region
1744 : !
1745 : ! Initialization, which includes setting HD to H times XOPT.
1746 : !
1747 :
1748 45823474 : delsq = delta*delta
1749 45823474 : iterc = 0
1750 45823474 : itermax = n
1751 45823474 : itersw = itermax
1752 137474763 : DO i = 1, n
1753 137474763 : d(i) = xopt(i)
1754 : END DO
1755 45823474 : CALL updatehd
1756 : !
1757 : ! Prepare for the first line search.
1758 : !
1759 45823474 : qred = zero
1760 45823474 : dd = zero
1761 137474763 : DO i = 1, n
1762 91651289 : step(i) = zero
1763 91651289 : hs(i) = zero
1764 91651289 : g(i) = gq(i) + hd(i)
1765 91651289 : d(i) = -g(i)
1766 137474763 : dd = dd + d(i)**2
1767 : END DO
1768 45823474 : crvmin = zero
1769 45823474 : IF (dd == zero) RETURN
1770 : ds = zero
1771 : ss = zero
1772 : gg = dd
1773 : ggbeg = gg
1774 : !
1775 : ! Calculate the step to the trust region boundary and the product HD
1776 : !
1777 : jump1 = .FALSE.
1778 : jump2 = .FALSE.
1779 : mainloop: DO
1780 77350711 : IF (.NOT. jump2) THEN
1781 77349944 : IF (.NOT. jump1) THEN
1782 77349944 : iterc = iterc + 1
1783 77349944 : temp = delsq - ss
1784 77349944 : bstep = temp/(ds + SQRT(ds*ds + dd*temp))
1785 77349944 : CALL updatehd
1786 : END IF
1787 77349944 : jump1 = .FALSE.
1788 77349944 : IF (iterc <= itersw) THEN
1789 77349944 : dhd = zero
1790 232063028 : DO j = 1, n
1791 232063028 : dhd = dhd + d(j)*hd(j)
1792 : END DO
1793 : !
1794 : ! Update CRVMIN and set the step-length ALPHA.
1795 : !
1796 77349944 : alpha = bstep
1797 77349944 : IF (dhd > zero) THEN
1798 68358138 : temp = dhd/dd
1799 68358138 : IF (iterc == 1) crvmin = temp
1800 68358138 : crvmin = MIN(crvmin, temp)
1801 68358138 : alpha = MIN(alpha, gg/dhd)
1802 : END IF
1803 77349944 : qadd = alpha*(gg - half*alpha*dhd)
1804 77349944 : qred = qred + qadd
1805 : !
1806 : ! Update STEP and HS.
1807 : !
1808 77349944 : ggsav = gg
1809 77349944 : gg = zero
1810 232063028 : DO i = 1, n
1811 154713084 : step(i) = step(i) + alpha*d(i)
1812 154713084 : hs(i) = hs(i) + alpha*hd(i)
1813 232063028 : gg = gg + (g(i) + hs(i))**2
1814 : END DO
1815 : !
1816 : ! Begin another conjugate direction iteration if required.
1817 : !
1818 77349944 : IF (alpha < bstep) THEN
1819 50847910 : IF (qadd <= 0.01_dp*qred) EXIT mainloop
1820 50076577 : IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1821 31526675 : IF (iterc == itermax) EXIT mainloop
1822 31526623 : temp = gg/ggsav
1823 31526623 : dd = zero
1824 31526623 : ds = zero
1825 31526623 : ss = zero
1826 94589055 : DO i = 1, n
1827 63062432 : d(i) = temp*d(i) - g(i) - hs(i)
1828 63062432 : dd = dd + d(i)**2
1829 63062432 : ds = ds + d(i)*step(I)
1830 94589055 : ss = ss + step(i)**2
1831 : END DO
1832 31526623 : IF (ds <= zero) EXIT mainloop
1833 31526623 : IF (ss < delsq) CYCLE mainloop
1834 : END IF
1835 26502034 : crvmin = zero
1836 26502034 : itersw = iterc
1837 26502034 : jump2 = .TRUE.
1838 26502034 : IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1839 : ELSE
1840 : jump2 = .FALSE.
1841 : END IF
1842 : END IF
1843 : !
1844 : ! Test whether an alternative iteration is required.
1845 : !
1846 : !!!! IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1847 : IF (jump2) THEN
1848 26106955 : sg = zero
1849 26106955 : shs = zero
1850 78326318 : DO i = 1, n
1851 52219363 : sg = sg + step(i)*g(i)
1852 78326318 : shs = shs + step(i)*hs(i)
1853 : END DO
1854 26106955 : sgk = sg + shs
1855 26106955 : angtest = sgk/SQRT(gg*delsq)
1856 26106955 : IF (angtest <= -0.99_dp) EXIT mainloop
1857 : !
1858 : ! Begin the alternative iteration by calculating D and HD and some
1859 : ! scalar products.
1860 : !
1861 22794760 : iterc = iterc + 1
1862 22794760 : temp = SQRT(delsq*gg - sgk*sgk)
1863 22794760 : tempa = delsq/temp
1864 22794760 : tempb = sgk/temp
1865 68389297 : DO i = 1, n
1866 68389297 : d(i) = tempa*(g(i) + hs(i)) - tempb*step(i)
1867 : END DO
1868 22794760 : CALL updatehd
1869 22794760 : IF (iterc <= itersw) THEN
1870 : jump1 = .TRUE.
1871 : CYCLE mainloop
1872 : END IF
1873 : END IF
1874 22794760 : dg = zero
1875 22794760 : dhd = zero
1876 22794760 : dhs = zero
1877 68389297 : DO i = 1, n
1878 45594537 : dg = dg + d(i)*g(i)
1879 45594537 : dhd = dhd + hd(i)*d(i)
1880 68389297 : dhs = dhs + hd(i)*step(i)
1881 : END DO
1882 : !
1883 : ! Seek the value of the angle that minimizes Q.
1884 : !
1885 22794760 : cf = half*(shs - dhd)
1886 22794760 : qbeg = sg + cf
1887 22794760 : qsav = qbeg
1888 22794760 : qmin = qbeg
1889 22794760 : isave = 0
1890 22794760 : iu = 49
1891 : temp = twopi/REAL(iu + 1, dp)
1892 1139738000 : DO i = 1, iu
1893 1116943240 : angle = REAL(i, dp)*temp
1894 1116943240 : cth = COS(angle)
1895 1116943240 : sth = SIN(angle)
1896 1116943240 : qnew = (sg + cf*cth)*cth + (dg + dhs*cth)*sth
1897 1116943240 : IF (qnew < qmin) THEN
1898 : qmin = qnew
1899 : isave = i
1900 : tempa = qsav
1901 1093200935 : ELSE IF (i == isave + 1) THEN
1902 28439255 : tempb = qnew
1903 : END IF
1904 1139738000 : qsav = qnew
1905 : END DO
1906 22794760 : IF (isave == zero) tempa = qnew
1907 9879702 : IF (isave == iu) tempb = qbeg
1908 22794760 : angle = zero
1909 22794760 : IF (tempa /= tempb) THEN
1910 22794760 : tempa = tempa - qmin
1911 22794760 : tempb = tempb - qmin
1912 22794760 : angle = half*(tempa - tempb)/(tempa + tempb)
1913 : END IF
1914 22794760 : angle = temp*(REAL(isave, DP) + angle)
1915 : !
1916 : ! Calculate the new STEP and HS. Then test for convergence.
1917 : !
1918 22794760 : cth = COS(angle)
1919 22794760 : sth = SIN(angle)
1920 22794760 : reduc = qbeg - (sg + cf*cth)*cth - (dg + dhs*cth)*sth
1921 22794760 : gg = zero
1922 68389297 : DO i = 1, n
1923 45594537 : step(i) = cth*step(i) + sth*d(i)
1924 45594537 : hs(i) = cth*hs(i) + sth*hd(i)
1925 68389297 : gg = gg + (g(i) + hs(i))**2
1926 : END DO
1927 22794760 : qred = qred + reduc
1928 22794760 : ratio = reduc/qred
1929 22794760 : IF (iterc < itermax .AND. ratio > 0.01_dp) THEN
1930 767 : jump2 = .TRUE.
1931 : ELSE
1932 : EXIT mainloop
1933 : END IF
1934 :
1935 45824088 : IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1936 :
1937 : END DO mainloop
1938 :
1939 : CONTAINS
1940 : ! **************************************************************************************************
1941 : !> \brief ...
1942 : ! **************************************************************************************************
1943 145968178 : SUBROUTINE updatehd
1944 : INTEGER :: i, ih, j, k
1945 :
1946 437927088 : DO i = 1, n
1947 437927088 : hd(i) = zero
1948 : END DO
1949 875854176 : DO k = 1, npt
1950 729885998 : temp = zero
1951 2189923376 : DO j = 1, n
1952 2189923376 : temp = temp + xpt(k, j)*d(j)
1953 : END DO
1954 729885998 : temp = temp*pq(k)
1955 2335891554 : DO i = 1, n
1956 2189923376 : hd(i) = hd(i) + temp*xpt(k, i)
1957 : END DO
1958 : END DO
1959 145968178 : ih = 0
1960 437927088 : DO j = 1, n
1961 875926160 : DO i = 1, j
1962 437999072 : ih = ih + 1
1963 437999072 : IF (i < j) hd(j) = hd(j) + hq(ih)*d(i)
1964 729957982 : hd(i) = hd(i) + hq(ih)*d(j)
1965 : END DO
1966 : END DO
1967 145968178 : END SUBROUTINE updatehd
1968 :
1969 : END SUBROUTINE trsapp
1970 :
1971 : ! **************************************************************************************************
1972 : !> \brief ...
1973 : !> \param n ...
1974 : !> \param npt ...
1975 : !> \param bmat ...
1976 : !> \param zmat ...
1977 : !> \param idz ...
1978 : !> \param ndim ...
1979 : !> \param vlag ...
1980 : !> \param beta ...
1981 : !> \param knew ...
1982 : !> \param w ...
1983 : ! **************************************************************************************************
1984 53176237 : SUBROUTINE update(n, npt, bmat, zmat, idz, ndim, vlag, beta, knew, w)
1985 :
1986 : INTEGER, INTENT(IN) :: n, npt, ndim
1987 : INTEGER, INTENT(INOUT) :: idz
1988 : REAL(dp), DIMENSION(npt, *), INTENT(INOUT) :: zmat
1989 : REAL(dp), DIMENSION(ndim, *), INTENT(INOUT) :: bmat
1990 : REAL(dp), DIMENSION(*), INTENT(INOUT) :: vlag
1991 : REAL(dp), INTENT(INOUT) :: beta
1992 : INTEGER, INTENT(INOUT) :: knew
1993 : REAL(dp), DIMENSION(*), INTENT(INOUT) :: w
1994 :
1995 : REAL(dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
1996 :
1997 : INTEGER :: i, iflag, j, ja, jb, jl, jp, nptm
1998 : REAL(dp) :: alpha, denom, scala, scalb, tau, tausq, &
1999 : temp, tempa, tempb
2000 :
2001 : ! The arrays BMAT and ZMAT with IDZ are updated, in order to shift the
2002 : ! interpolation point that has index KNEW. On entry, VLAG contains the
2003 : ! components of the vector Theta*Wcheck+e_b of the updating formula
2004 : ! (6.11), and BETA holds the value of the parameter that has this na
2005 : ! The vector W is used for working space.
2006 : !
2007 :
2008 53176237 : nptm = npt - n - 1
2009 : !
2010 : ! Apply the rotations that put zeros in the KNEW-th row of ZMAT.
2011 : !
2012 53176237 : jl = 1
2013 106358273 : DO j = 2, nptm
2014 106358273 : IF (j == idz) THEN
2015 : jl = idz
2016 53182036 : ELSE IF (zmat(knew, j) /= zero) THEN
2017 51737992 : temp = SQRT(zmat(knew, jl)**2 + zmat(knew, j)**2)
2018 51737992 : tempa = zmat(knew, jl)/temp
2019 51737992 : tempb = zmat(knew, j)/temp
2020 310474944 : DO I = 1, NPT
2021 258736952 : temp = tempa*zmat(i, jl) + tempb*zmat(i, j)
2022 258736952 : zmat(i, j) = tempa*zmat(i, j) - tempb*zmat(i, jl)
2023 310474944 : zmat(i, jl) = temp
2024 : END DO
2025 51737992 : zmat(knew, j) = zero
2026 : END IF
2027 : END DO
2028 : !
2029 : ! Put the first NPT components of the KNEW-th column of HLAG into W,
2030 : ! and calculate the parameters of the updating formula.
2031 : !
2032 53176237 : tempa = zmat(knew, 1)
2033 53176237 : IF (idz >= 2) tempa = -tempa
2034 53176237 : IF (jl > 1) tempb = zmat(knew, jl)
2035 319069020 : DO i = 1, npt
2036 265892783 : w(i) = tempa*zmat(i, 1)
2037 319069020 : IF (jl > 1) w(i) = w(i) + tempb*zmat(i, jl)
2038 : END DO
2039 53176237 : alpha = w(knew)
2040 53176237 : tau = vlag(knew)
2041 53176237 : tausq = tau*tau
2042 53176237 : denom = alpha*beta + tausq
2043 53176237 : vlag(knew) = vlag(knew) - one
2044 : !
2045 : ! Complete the updating of ZMAT when there is only one nonzero eleme
2046 : ! in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to
2047 : ! then the first column of ZMAT will be exchanged with another one l
2048 : !
2049 53176237 : iflag = 0
2050 53176237 : IF (JL == 1) THEN
2051 53176237 : temp = SQRT(ABS(denom))
2052 53176237 : tempb = tempa/temp
2053 53176237 : tempa = tau/temp
2054 319069020 : DO i = 1, npt
2055 319069020 : zmat(i, 1) = tempa*zmat(i, 1) - tempb*vlag(i)
2056 : END DO
2057 53176237 : IF (idz == 1 .AND. temp < zero) idz = 2
2058 53176237 : IF (idz >= 2 .AND. temp >= zero) iflag = 1
2059 : ELSE
2060 : !
2061 : ! Complete the updating of ZMAT in the alternative case.
2062 : !
2063 0 : ja = 1
2064 0 : IF (beta >= zero) ja = jl
2065 0 : jb = jl + 1 - ja
2066 0 : temp = zmat(knew, jb)/denom
2067 0 : tempa = temp*beta
2068 0 : tempb = temp*tau
2069 0 : temp = zmat(knew, ja)
2070 0 : scala = one/SQRT(ABS(beta)*temp*temp + tausq)
2071 0 : scalb = scala*SQRT(ABS(denom))
2072 0 : DO i = 1, npt
2073 0 : zmat(i, ja) = scala*(tau*zmat(i, ja) - temp*vlag(i))
2074 0 : zmat(i, jb) = scalb*(zmat(i, jb) - tempa*w(i) - tempb*vlag(i))
2075 : END DO
2076 0 : IF (denom <= zero) THEN
2077 0 : IF (beta < zero) idz = idz + 1
2078 0 : IF (beta >= zero) iflag = 1
2079 : END IF
2080 : END IF
2081 : !
2082 : ! IDZ is reduced in the following case, and usually the first column
2083 : ! of ZMAT is exchanged with a later one.
2084 : !
2085 : IF (iflag == 1) THEN
2086 0 : idz = idz - 1
2087 0 : DO i = 1, npt
2088 0 : temp = zmat(i, 1)
2089 0 : zmat(i, 1) = zmat(i, idz)
2090 0 : zmat(i, idz) = temp
2091 : END DO
2092 : END IF
2093 : !
2094 : ! Finally, update the matrix BMAT.
2095 : !
2096 159534510 : DO j = 1, n
2097 106358273 : jp = npt + j
2098 106358273 : w(jp) = bmat(knew, j)
2099 106358273 : tempa = (alpha*vlag(jp) - tau*w(jp))/denom
2100 106358273 : tempb = (-beta*w(jp) - tau*vlag(jp))/denom
2101 850939947 : DO i = 1, jp
2102 691405437 : bmat(i, j) = bmat(i, j) + tempa*vlag(i) + tempb*w(i)
2103 797763710 : IF (i > npt) bmat(jp, i - npt) = bmat(i, j)
2104 : END DO
2105 : END DO
2106 :
2107 53176237 : END SUBROUTINE update
2108 :
2109 0 : END MODULE powell
|
