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 64030364 : 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 64030364 : CALL timeset(routineN, handle)
61 :
62 64845697 : SELECT CASE (optstate%state)
63 : CASE (0)
64 815333 : npt = 2*n + 1
65 2445999 : ALLOCATE (optstate%w((npt + 13)*(npt + n) + 3*n*(n + 3)/2))
66 2445999 : ALLOCATE (optstate%xopt(n))
67 : ! Initialize w
68 115921415 : optstate%w = 0.0_dp
69 815333 : optstate%state = 1
70 815333 : CALL newuoa(n, x, optstate)
71 : CASE (1, 2)
72 61584370 : CALL newuoa(n, x, optstate)
73 : CASE (3)
74 11 : IF (optstate%unit > 0) THEN
75 7 : WRITE (optstate%unit, *) "POWELL| Exceeding maximum number of steps"
76 : END IF
77 11 : optstate%state = -1
78 : CASE (4)
79 0 : IF (optstate%unit > 0) THEN
80 0 : WRITE (optstate%unit, *) "POWELL| Error in trust region"
81 : END IF
82 0 : 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 815317 : optstate%state = -1
90 : CASE (8)
91 2446353 : x(1:n) = optstate%xopt(1:n)
92 815333 : DEALLOCATE (optstate%w)
93 815333 : DEALLOCATE (optstate%xopt)
94 815333 : optstate%state = -1
95 : CASE DEFAULT
96 64030364 : CPABORT("Unknown optimization state")
97 : END SELECT
98 :
99 64030364 : CALL timestop(handle)
100 :
101 64030364 : END SUBROUTINE powell_optimize
102 : ! **************************************************************************************************
103 : !> \brief ...
104 : !> \param n ...
105 : !> \param x ...
106 : !> \param optstate ...
107 : ! **************************************************************************************************
108 62399703 : 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 62399703 : maxfun = optstate%maxfun
120 62399703 : rhobeg = optstate%rhobeg
121 62399703 : 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 62399703 : np = n + 1
158 62399703 : npt = 2*n + 1
159 62399703 : nptm = npt - np
160 62399703 : IF (npt < n + 2 .OR. npt > ((n + 2)*np)/2) THEN
161 0 : optstate%state = 5
162 0 : RETURN
163 : END IF
164 62399703 : ndim = npt + n
165 62399703 : ixb = 1
166 62399703 : ixo = ixb + n
167 62399703 : ixn = ixo + n
168 62399703 : ixp = ixn + n
169 62399703 : ifv = ixp + n*npt
170 62399703 : igq = ifv + npt
171 62399703 : ihq = igq + n
172 62399703 : ipq = ihq + (n*np)/2
173 62399703 : ibmat = ipq + npt
174 62399703 : izmat = ibmat + ndim*n
175 62399703 : id = izmat + npt*nptm
176 62399703 : ivl = id + n
177 62399703 : 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 62399703 : optstate%w(ivl:), optstate%w(iw:), optstate)
187 :
188 187214345 : 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 62399703 : SUBROUTINE newuob(n, npt, x, rhobeg, rhoend, maxfun, xbase, &
217 62399703 : 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 62399703 : skip_check = .FALSE.
273 62399703 : do_return = .FALSE.
274 :
275 62399703 : 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 815333 : idz = 0
282 815333 : itest = 0
283 815333 : nf = 0
284 815333 : nfm = 0
285 815333 : nfmm = 0
286 815333 : nfsav = 0
287 815333 : knew = 0
288 815333 : kopt = 0
289 815333 : ksave = 0
290 815333 : ktemp = 0
291 815333 : rhosq = 0._dp
292 815333 : recip = 0._dp
293 815333 : reciq = 0._dp
294 815333 : fbeg = 0._dp
295 815333 : fopt = 0._dp
296 815333 : diffa = 0._dp
297 815333 : xoptsq = 0._dp
298 815333 : rho = 0._dp
299 815333 : delta = 0._dp
300 815333 : dsq = 0._dp
301 815333 : dnorm = 0._dp
302 815333 : ratio = 0._dp
303 815333 : temp = 0._dp
304 815333 : tempq = 0._dp
305 815333 : beta = 0._dp
306 815333 : dx = 0._dp
307 815333 : vquad = 0._dp
308 815333 : diff = 0._dp
309 815333 : diffc = 0._dp
310 815333 : diffb = 0._dp
311 815333 : fsave = 0._dp
312 815333 : detrat = 0._dp
313 815333 : hdiag = 0._dp
314 815333 : distsq = 0._dp
315 815333 : gisq = 0._dp
316 815333 : gqsq = 0._dp
317 815333 : f = 0._dp
318 815333 : bstep = 0._dp
319 815333 : alpha = 0._dp
320 815333 : dstep = 0._dp
321 : !
322 : END IF
323 :
324 62399703 : ipt = 0
325 62399703 : jpt = 0
326 62399703 : xipt = 0._dp
327 62399703 : xjpt = 0._dp
328 :
329 62399703 : half = 0.5_dp
330 62399703 : one = 1.0_dp
331 62399703 : tenth = 0.1_dp
332 62399703 : zero = 0.0_dp
333 62399703 : np = n + 1
334 62399703 : nh = (n*np)/2
335 62399703 : nptm = npt - np
336 62399703 : nftest = MAX(maxfun, 1)
337 :
338 62399703 : IF (opt%state /= 2) THEN
339 : !
340 : ! Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero.
341 : !
342 2446353 : DO j = 1, n
343 1631020 : xbase(j) = x(j)
344 9818560 : DO k = 1, npt
345 9818560 : xpt(k, j) = zero
346 : END DO
347 13912153 : DO i = 1, ndim
348 13096820 : bmat(i, j) = zero
349 : END DO
350 : END DO
351 3269973 : DO ih = 1, nh
352 3269973 : hq(ih) = zero
353 : END DO
354 4892706 : DO k = 1, npt
355 4077373 : pq(k) = zero
356 13080246 : DO j = 1, nptm
357 12264913 : 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 815333 : rhosq = rhobeg*rhobeg
366 815333 : recip = one/rhosq
367 815333 : reciq = SQRT(half)/rhosq
368 815333 : nf = 0
369 815333 : CALL begin_initialization_of_interpolation
370 25774785 : RETURN
371 : END IF
372 :
373 61584370 : CALL set_state
374 :
375 61584370 : preparation: IF (nf <= npt) THEN
376 4076896 : fval(nf) = f
377 4076896 : IF (nf == 1) THEN
378 815333 : fbeg = f
379 815333 : fopt = f
380 815333 : kopt = 1
381 3261563 : ELSE IF (f < fopt) THEN
382 910366 : fopt = f
383 910366 : 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 4076896 : IF (NFM <= 2*N) THEN
390 4076896 : IF (nfm >= 1 .AND. nfm <= n) THEN
391 1630814 : gq(nfm) = (f - fbeg)/rhobeg
392 1630814 : 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 2446082 : ELSE IF (nfm > n) THEN
398 1630749 : bmat(nf - n, nfmm) = half/rhobeg
399 1630749 : bmat(nf, nfmm) = -half/rhobeg
400 1630749 : zmat(1, nfmm) = -reciq - reciq
401 1630749 : zmat(nf - n, nfmm) = reciq
402 1630749 : zmat(nf, nfmm) = reciq
403 1630749 : ih = (nfmm*(nfmm + 1))/2
404 1630749 : temp = (fbeg - f)/rhobeg
405 1630749 : hq(ih) = (gq(nfmm) - temp)/rhobeg
406 1630749 : 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 4076896 : IF (nf < npt) THEN
423 3261572 : CALL begin_initialization_of_interpolation
424 3261572 : RETURN
425 : END IF
426 : !
427 : ! Begin the iterative procedure, because the initial model is comple
428 : !
429 815324 : rho = rhobeg
430 815324 : delta = rho
431 815324 : idz = 1
432 815324 : diffa = zero
433 815324 : diffb = zero
434 815324 : itest = 0
435 815324 : xoptsq = zero
436 2446073 : DO i = 1, n
437 1630749 : xopt(i) = xpt(kopt, i)
438 2446073 : xoptsq = xoptsq + xopt(i)**2
439 : END DO
440 815324 : 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 815324 : skip_check = .TRUE.
446 : ELSE preparation
447 57507474 : IF (knew == -1) THEN
448 680613 : opt%state = 6
449 680613 : CALL get_state
450 680613 : CALL finalize_optimization_after_failure
451 680613 : 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 56826861 : vquad = zero
458 56826861 : ih = 0
459 170488818 : DO j = 1, n
460 113661957 : vquad = vquad + d(j)*gq(j)
461 341016365 : DO i = 1, j
462 170527547 : ih = ih + 1
463 170527547 : temp = d(i)*xnew(j) + d(j)*xopt(i)
464 170527547 : IF (i == j) temp = half*temp
465 284189504 : vquad = vquad + temp*hq(ih)
466 : END DO
467 : END DO
468 340977636 : DO k = 1, npt
469 340977636 : vquad = vquad + pq(k)*w(k)
470 : END DO
471 56826861 : diff = f - fopt - vquad
472 56826861 : diffc = diffb
473 56826861 : diffb = diffa
474 56826861 : diffa = ABS(diff)
475 56826861 : 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 56826861 : fsave = fopt
482 56826861 : IF (f < fopt) THEN
483 28878807 : fopt = f
484 28878807 : xoptsq = zero
485 86639758 : DO i = 1, n
486 57760951 : xopt(i) = xnew(i)
487 86639758 : xoptsq = xoptsq + xopt(i)**2
488 : END DO
489 : END IF
490 56826861 : ksave = knew
491 56826861 : IF (knew <= 0) THEN
492 : !
493 : ! Pick the next value of DELTA after a trust region step.
494 : !
495 36624914 : IF (vquad >= zero) THEN
496 : ! Return because a trust region step has failed to reduce Q
497 0 : opt%state = 4
498 0 : CALL get_state
499 0 : CALL finalize_optimization_after_failure
500 0 : RETURN
501 : END IF
502 36624914 : ratio = (f - fsave)/vquad
503 36624914 : IF (ratio <= tenth) THEN
504 13674812 : delta = half*dnorm
505 22950102 : ELSE IF (ratio <= 0.7_dp) THEN
506 4002728 : delta = MAX(half*delta, dnorm)
507 : ELSE
508 18947374 : delta = MAX(half*delta, dnorm + dnorm)
509 : END IF
510 36624914 : 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 36624914 : rhosq = MAX(tenth*delta, rho)**2
515 36624914 : ktemp = 0
516 36624914 : detrat = zero
517 36624914 : IF (f >= fsave) THEN
518 11738879 : ktemp = kopt
519 11738879 : detrat = one
520 : END IF
521 219759782 : DO k = 1, npt
522 183134868 : hdiag = zero
523 549500147 : DO j = 1, nptm
524 366365279 : temp = one
525 366365279 : IF (j < idz) temp = -one
526 549500147 : hdiag = hdiag + temp*zmat(k, j)**2
527 : END DO
528 183134868 : temp = ABS(beta*hdiag + vlag(k)**2)
529 183134868 : distsq = zero
530 549500147 : DO j = 1, n
531 549500147 : distsq = distsq + (xpt(k, j) - xopt(j))**2
532 : END DO
533 183134868 : IF (distsq > rhosq) temp = temp*(distsq/rhosq)**3
534 219759782 : IF (temp > detrat .AND. k /= ktemp) THEN
535 78735579 : detrat = temp
536 78735579 : knew = k
537 : END IF
538 : END DO
539 36624914 : 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 56332932 : CALL update(n, npt, bmat, zmat, idz, ndim, vlag, beta, knew, w)
547 56332932 : fval(knew) = f
548 56332932 : ih = 0
549 169006995 : DO i = 1, n
550 112674063 : temp = pq(knew)*xpt(knew, i)
551 338052604 : DO j = 1, i
552 169045609 : ih = ih + 1
553 281719672 : hq(ih) = hq(ih) + temp*xpt(knew, j)
554 : END DO
555 : END DO
556 56332932 : 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 169006995 : DO j = 1, nptm
562 112674063 : temp = diff*zmat(knew, j)
563 112674063 : IF (j < idz) temp = -temp
564 732515368 : DO k = 1, npt
565 676182436 : pq(k) = pq(k) + temp*zmat(k, j)
566 : END DO
567 : END DO
568 56332932 : gqsq = zero
569 169006995 : DO i = 1, n
570 112674063 : gq(i) = gq(i) + diff*bmat(knew, i)
571 112674063 : gqsq = gqsq + gq(i)**2
572 169006995 : 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 56332932 : IF (ksave == 0 .AND. delta == rho) THEN
580 7331185 : IF (ABS(ratio) > 1.0e-2_dp) THEN
581 4921771 : itest = 0
582 : ELSE
583 14456578 : DO k = 1, npt
584 14456578 : vlag(k) = fval(k) - fval(kopt)
585 : END DO
586 2409414 : gisq = zero
587 7228289 : DO i = 1, n
588 4818875 : sum = zero
589 28914176 : DO k = 1, npt
590 28914176 : sum = sum + bmat(k, i)*vlag(k)
591 : END DO
592 4818875 : gisq = gisq + sum*sum
593 7228289 : 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 2409414 : itest = itest + 1
600 2409414 : IF (gqsq < 1.0e2_dp*gisq) itest = 0
601 2409414 : IF (itest >= 3) THEN
602 394434 : DO i = 1, n
603 394434 : gq(i) = w(i)
604 : END DO
605 525912 : DO ih = 1, nh
606 525912 : hq(ih) = zero
607 : END DO
608 394434 : DO j = 1, nptm
609 262956 : w(j) = zero
610 1577736 : DO k = 1, npt
611 1577736 : w(j) = w(j) + vlag(k)*zmat(k, j)
612 : END DO
613 394434 : IF (j < idz) w(j) = -w(j)
614 : END DO
615 788868 : DO k = 1, npt
616 657390 : pq(k) = zero
617 2103648 : DO j = 1, nptm
618 1972170 : pq(k) = pq(k) + zmat(k, j)*w(j)
619 : END DO
620 : END DO
621 131478 : itest = 0
622 : END IF
623 : END IF
624 : END IF
625 56332932 : 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 56332932 : IF (f <= fsave + tenth*vquad .OR. ksave > 0) THEN
632 : skip_check = .TRUE.
633 : ELSE
634 13180873 : 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 67839957 : IF (.NOT. skip_check) THEN
643 23872574 : skip_check = .FALSE.
644 23872574 : distsq = 4.0_dp*delta*delta
645 143242320 : DO k = 1, npt
646 119369746 : sum = zero
647 358165544 : DO j = 1, n
648 358165544 : sum = sum + (xpt(k, j) - xopt(j))**2
649 : END DO
650 143242320 : IF (sum > distsq) THEN
651 32595256 : knew = k
652 32595256 : 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 23872574 : IF (knew > 0) THEN
660 20201950 : dstep = MAX(MIN(tenth*SQRT(distsq), half*delta), rho)
661 20201950 : dsq = dstep*dstep
662 20201950 : CALL generate_next_model_step
663 20201950 : RETURN
664 : END IF
665 3670624 : 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 3202034 : CALL update_rho()
671 3202034 : IF (do_return) RETURN
672 : END IF
673 : END IF
674 :
675 48512844 : skip_check = .FALSE.
676 : inner: DO
677 48512844 : knew = 0
678 48512844 : CALL trsapp(n, npt, xopt, xpt, gq, hq, pq, delta, d, w, w(np), w(np + n), w(np + 2*n), crvmin)
679 48512844 : dsq = zero
680 145544523 : DO i = 1, n
681 145544523 : dsq = dsq + d(i)**2
682 : END DO
683 48512844 : dnorm = MIN(delta, SQRT(dsq))
684 48512844 : IF (dnorm < half*rho) THEN
685 11887926 : knew = -1
686 11887926 : delta = tenth*delta
687 11887926 : ratio = -1.0_dp
688 11887926 : IF (delta <= 1.5_dp*rho) delta = rho
689 11887926 : IF (nf <= nfsav + 2) CYCLE outer
690 3310881 : temp = 0.125_dp*crvmin*rho*rho
691 3310881 : IF (temp <= MAX(diffa, diffb, diffc)) CYCLE outer
692 1690154 : CALL update_rho()
693 1690154 : IF (do_return) RETURN
694 : CYCLE inner
695 : END IF
696 : EXIT inner
697 : END DO inner
698 10197772 : 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 36624918 : CALL generate_next_model_step
704 :
705 : CONTAINS
706 : ! **************************************************************************************************
707 : !> \brief ...
708 : ! **************************************************************************************************
709 4892188 : SUBROUTINE update_rho()
710 4892188 : IF (rho > rhoend) THEN
711 4076871 : delta = half*rho
712 4076871 : ratio = rho/rhoend
713 4076871 : IF (ratio <= 16.0_dp) THEN
714 815307 : rho = rhoend
715 3261564 : ELSE IF (ratio <= 250.0_dp) THEN
716 815307 : rho = SQRT(ratio)*rhoend
717 : ELSE
718 2446257 : rho = tenth*rho
719 : END IF
720 4076871 : delta = MAX(delta, rho)
721 4076871 : 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 815317 : IF (knew == -1) THEN
728 680613 : CALL calc_next_value_of_objective_func
729 : ELSE
730 134704 : opt%state = 7
731 134704 : CALL get_state
732 :
733 134704 : CALL finalize_optimization_after_failure
734 : END IF
735 815317 : do_return = .TRUE.
736 : END IF
737 4892188 : END SUBROUTINE update_rho
738 : ! **************************************************************************************************
739 : !> \brief ...
740 : ! **************************************************************************************************
741 4076905 : SUBROUTINE begin_initialization_of_interpolation()
742 4076905 : nfm = nf
743 4076905 : nfmm = nf - n
744 4076905 : nf = nf + 1
745 4076905 : IF (nfm <= 2*n) THEN
746 4076905 : IF (nfm >= 1 .AND. nfm <= N) THEN
747 1630823 : xpt(nf, nfm) = rhobeg
748 2446082 : ELSE IF (nfm > n) THEN
749 1630749 : 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 12237575 : DO j = 1, n
773 12237575 : x(j) = xpt(nf, j) + xbase(j)
774 : END DO
775 4076905 : CALL check_number_of_steps
776 4076905 : END SUBROUTINE begin_initialization_of_interpolation
777 : ! **************************************************************************************************
778 : !> \brief ...
779 : ! **************************************************************************************************
780 56826868 : SUBROUTINE generate_next_model_step()
781 56826868 : IF (dsq <= 1.0e-3_dp*xoptsq) THEN
782 3464347 : tempq = 0.25_dp*xoptsq
783 20786394 : DO k = 1, npt
784 17322047 : sum = zero
785 51968289 : DO i = 1, n
786 51968289 : sum = sum + xpt(k, i)*xopt(i)
787 : END DO
788 17322047 : temp = pq(k)*sum
789 17322047 : sum = sum - half*xoptsq
790 17322047 : w(npt + k) = sum
791 55432636 : DO i = 1, n
792 34646242 : gq(i) = gq(i) + temp*xpt(k, i)
793 34646242 : xpt(k, i) = xpt(k, i) - half*xopt(i)
794 34646242 : vlag(i) = bmat(k, i)
795 34646242 : w(i) = sum*xpt(k, i) + tempq*xopt(i)
796 34646242 : ip = npt + i
797 103946982 : DO j = 1, i
798 86624935 : 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 10393197 : DO k = 1, nptm
806 6928850 : sumz = zero
807 41575092 : DO i = 1, npt
808 34646242 : sumz = sumz + zmat(i, k)
809 41575092 : w(i) = w(npt + i)*zmat(i, k)
810 : END DO
811 20787546 : DO j = 1, n
812 13858696 : sum = tempq*sumz*xopt(j)
813 83169840 : DO i = 1, npt
814 69311144 : sum = sum + w(i)*xpt(i, j)
815 69311144 : vlag(j) = sum
816 83169840 : IF (k < idz) sum = -sum
817 : END DO
818 90098690 : DO i = 1, npt
819 83169840 : bmat(i, j) = bmat(i, j) + sum*zmat(i, k)
820 : END DO
821 : END DO
822 24251893 : DO i = 1, n
823 13858696 : ip = i + npt
824 13858696 : temp = vlag(i)
825 13858696 : IF (k < idz) temp = -temp
826 41580006 : DO j = 1, i
827 34651156 : 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 3464347 : ih = 0
836 10393197 : DO j = 1, n
837 6928850 : w(j) = zero
838 41575092 : DO k = 1, npt
839 34646242 : w(j) = w(j) + pq(k)*xpt(k, j)
840 41575092 : xpt(k, j) = xpt(k, j) - half*xopt(j)
841 : END DO
842 20786970 : DO i = 1, j
843 10393773 : ih = ih + 1
844 10393773 : IF (i < j) gq(j) = gq(j) + hq(ih)*xopt(i)
845 10393773 : gq(i) = gq(i) + hq(ih)*xopt(j)
846 10393773 : hq(ih) = hq(ih) + w(i)*xopt(j) + xopt(i)*w(j)
847 17322623 : bmat(npt + i, j) = bmat(npt + j, i)
848 : END DO
849 : END DO
850 10393197 : DO j = 1, n
851 6928850 : xbase(j) = xbase(j) + xopt(j)
852 10393197 : xopt(j) = zero
853 : END DO
854 3464347 : 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 56826868 : IF (knew > 0) THEN
862 : CALL biglag(n, npt, xopt, xpt, bmat, zmat, idz, ndim, knew, dstep, &
863 20201950 : 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 340977752 : DO k = 1, npt
870 284150884 : suma = zero
871 284150884 : sumb = zero
872 284150884 : sum = zero
873 852600592 : DO j = 1, n
874 568449708 : suma = suma + xpt(k, j)*d(j)
875 568449708 : sumb = sumb + xpt(k, j)*xopt(j)
876 852600592 : sum = sum + bmat(k, j)*d(j)
877 : END DO
878 284150884 : w(k) = suma*(half*suma + sumb)
879 340977752 : vlag(k) = sum
880 : END DO
881 56826868 : beta = zero
882 170488876 : DO k = 1, nptm
883 113662008 : sum = zero
884 682111716 : DO i = 1, npt
885 682111716 : sum = sum + zmat(i, k)*w(i)
886 : END DO
887 113662008 : IF (k < idz) THEN
888 0 : beta = beta + sum*sum
889 0 : sum = -sum
890 : ELSE
891 113662008 : beta = beta - sum*sum
892 : END IF
893 738938584 : DO i = 1, npt
894 682111716 : vlag(i) = vlag(i) + sum*zmat(i, k)
895 : END DO
896 : END DO
897 56826868 : bsum = zero
898 56826868 : dx = zero
899 170488876 : DO j = 1, n
900 113662008 : sum = zero
901 682111716 : DO i = 1, npt
902 682111716 : sum = sum + w(i)*bmat(i, j)
903 : END DO
904 113662008 : bsum = bsum + sum*d(j)
905 113662008 : jp = npt + j
906 341055858 : DO k = 1, n
907 341055858 : sum = sum + bmat(jp, k)*d(k)
908 : END DO
909 113662008 : vlag(jp) = sum
910 113662008 : bsum = bsum + sum*d(j)
911 170488876 : dx = dx + d(j)*xopt(j)
912 : END DO
913 56826868 : beta = dx*dx + dsq*(xoptsq + dx + dx + half*dsq) + beta - bsum
914 56826868 : 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 56826868 : IF (knew > 0) THEN
921 20201950 : temp = one + alpha*beta/vlag(knew)**2
922 20201950 : IF (ABS(temp) <= 0.8_dp) THEN
923 : CALL bigden(n, npt, xopt, xpt, bmat, zmat, idz, ndim, kopt, &
924 72 : 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 56826868 : CALL calc_next_value_of_objective_func
931 56826868 : END SUBROUTINE generate_next_model_step
932 : ! **************************************************************************************************
933 : !> \brief ...
934 : ! **************************************************************************************************
935 57507481 : SUBROUTINE calc_next_value_of_objective_func()
936 172530755 : DO i = 1, n
937 115023274 : xnew(i) = xopt(i) + d(i)
938 172530755 : x(i) = xbase(i) + xnew(i)
939 : END DO
940 57507481 : nf = nf + 1
941 57507481 : CALL check_number_of_steps
942 57507481 : END SUBROUTINE calc_next_value_of_objective_func
943 : ! **************************************************************************************************
944 : !> \brief ...
945 : ! **************************************************************************************************
946 61584386 : SUBROUTINE check_number_of_steps()
947 61584386 : IF (nf > nftest) THEN
948 : ! return to many steps
949 15 : nf = nf - 1
950 15 : opt%state = 3
951 15 : CALL get_state
952 15 : CALL finalize_optimization_after_failure
953 : ELSE
954 :
955 61584371 : CALL get_state
956 :
957 61584371 : opt%state = 2
958 : END IF
959 61584386 : END SUBROUTINE check_number_of_steps
960 : ! **************************************************************************************************
961 : !> \brief ...
962 : ! **************************************************************************************************
963 815332 : SUBROUTINE finalize_optimization_after_failure()
964 815332 : IF (fopt <= f) THEN
965 757557 : DO i = 1, n
966 757557 : x(i) = xbase(i) + xopt(i)
967 : END DO
968 252402 : f = fopt
969 : END IF
970 :
971 815332 : CALL get_state
972 815332 : END SUBROUTINE finalize_optimization_after_failure
973 : ! **************************************************************************************************
974 : !> \brief ...
975 : ! **************************************************************************************************
976 63215035 : SUBROUTINE get_state()
977 63215035 : opt%np = np
978 63215035 : opt%nh = nh
979 63215035 : opt%nptm = nptm
980 63215035 : opt%nftest = nftest
981 63215035 : opt%idz = idz
982 63215035 : opt%itest = itest
983 63215035 : opt%nf = nf
984 63215035 : opt%nfm = nfm
985 63215035 : opt%nfmm = nfmm
986 63215035 : opt%nfsav = nfsav
987 63215035 : opt%knew = knew
988 63215035 : opt%kopt = kopt
989 63215035 : opt%ksave = ksave
990 63215035 : opt%ktemp = ktemp
991 63215035 : opt%rhosq = rhosq
992 63215035 : opt%recip = recip
993 63215035 : opt%reciq = reciq
994 63215035 : opt%fbeg = fbeg
995 63215035 : opt%fopt = fopt
996 63215035 : opt%diffa = diffa
997 63215035 : opt%xoptsq = xoptsq
998 63215035 : opt%rho = rho
999 63215035 : opt%delta = delta
1000 63215035 : opt%dsq = dsq
1001 63215035 : opt%dnorm = dnorm
1002 63215035 : opt%ratio = ratio
1003 63215035 : opt%temp = temp
1004 63215035 : opt%tempq = tempq
1005 63215035 : opt%beta = beta
1006 63215035 : opt%dx = dx
1007 63215035 : opt%vquad = vquad
1008 63215035 : opt%diff = diff
1009 63215035 : opt%diffc = diffc
1010 63215035 : opt%diffb = diffb
1011 63215035 : opt%fsave = fsave
1012 63215035 : opt%detrat = detrat
1013 63215035 : opt%hdiag = hdiag
1014 63215035 : opt%distsq = distsq
1015 63215035 : opt%gisq = gisq
1016 63215035 : opt%gqsq = gqsq
1017 63215035 : opt%f = f
1018 63215035 : opt%bstep = bstep
1019 63215035 : opt%alpha = alpha
1020 63215035 : opt%dstep = dstep
1021 63215035 : END SUBROUTINE get_state
1022 : ! **************************************************************************************************
1023 : !> \brief ...
1024 : ! **************************************************************************************************
1025 61584370 : SUBROUTINE set_state()
1026 61584370 : np = opt%np
1027 61584370 : nh = opt%nh
1028 61584370 : nptm = opt%nptm
1029 61584370 : nftest = opt%nftest
1030 61584370 : idz = opt%idz
1031 61584370 : itest = opt%itest
1032 61584370 : nf = opt%nf
1033 61584370 : nfm = opt%nfm
1034 61584370 : nfmm = opt%nfmm
1035 61584370 : nfsav = opt%nfsav
1036 61584370 : knew = opt%knew
1037 61584370 : kopt = opt%kopt
1038 61584370 : ksave = opt%ksave
1039 61584370 : ktemp = opt%ktemp
1040 61584370 : rhosq = opt%rhosq
1041 61584370 : recip = opt%recip
1042 61584370 : reciq = opt%reciq
1043 61584370 : fbeg = opt%fbeg
1044 61584370 : fopt = opt%fopt
1045 61584370 : diffa = opt%diffa
1046 61584370 : xoptsq = opt%xoptsq
1047 61584370 : rho = opt%rho
1048 61584370 : delta = opt%delta
1049 61584370 : dsq = opt%dsq
1050 61584370 : dnorm = opt%dnorm
1051 61584370 : ratio = opt%ratio
1052 61584370 : temp = opt%temp
1053 61584370 : tempq = opt%tempq
1054 61584370 : beta = opt%beta
1055 61584370 : dx = opt%dx
1056 61584370 : vquad = opt%vquad
1057 61584370 : diff = opt%diff
1058 61584370 : diffc = opt%diffc
1059 61584370 : diffb = opt%diffb
1060 61584370 : fsave = opt%fsave
1061 61584370 : detrat = opt%detrat
1062 61584370 : hdiag = opt%hdiag
1063 61584370 : distsq = opt%distsq
1064 61584370 : gisq = opt%gisq
1065 61584370 : gqsq = opt%gqsq
1066 61584370 : f = opt%f
1067 61584370 : bstep = opt%bstep
1068 61584370 : alpha = opt%alpha
1069 61584370 : dstep = opt%dstep
1070 61584370 : 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 72 : SUBROUTINE bigden(n, npt, xopt, xpt, bmat, zmat, idz, ndim, kopt, &
1095 72 : 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 72 : 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 432 : DO k = 1, npt
1150 432 : w(n + k) = zero
1151 : END DO
1152 216 : DO j = 1, nptm
1153 144 : temp = zmat(knew, j)
1154 144 : IF (j < idz) temp = -temp
1155 936 : DO k = 1, npt
1156 864 : w(n + k) = w(n + k) + temp*zmat(k, j)
1157 : END DO
1158 : END DO
1159 72 : 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 72 : dd = zero
1167 72 : ds = zero
1168 72 : ss = zero
1169 72 : xoptsq = zero
1170 216 : DO i = 1, n
1171 144 : dd = dd + d(i)**2
1172 144 : s(i) = xpt(knew, i) - xopt(i)
1173 144 : ds = ds + d(i)*s(i)
1174 144 : ss = ss + s(i)**2
1175 216 : xoptsq = xoptsq + xopt(i)**2
1176 : END DO
1177 72 : 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 72 : ssden = dd*ss - ds*ds
1202 72 : iterc = 0
1203 72 : 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 72 : iterc = iterc + 1
1210 72 : temp = one/SQRT(ssden)
1211 72 : xoptd = zero
1212 72 : xopts = zero
1213 216 : DO i = 1, n
1214 144 : s(i) = temp*(dd*s(i) - ds*d(i))
1215 144 : xoptd = xoptd + xopt(i)*d(i)
1216 216 : xopts = xopts + xopt(i)*s(i)
1217 : END DO
1218 : !
1219 : ! Set the coefficients of the first two terms of BETA.
1220 : !
1221 72 : tempa = half*xoptd*xoptd
1222 72 : tempb = half*xopts*xopts
1223 72 : den(1) = dd*(xoptsq + half*dd) + tempa + tempb
1224 72 : den(2) = two*xoptd*dd
1225 72 : den(3) = two*xopts*dd
1226 72 : den(4) = tempa - tempb
1227 72 : den(5) = xoptd*xopts
1228 360 : DO i = 6, 9
1229 360 : den(i) = zero
1230 : END DO
1231 : !
1232 : ! Put the coefficients of Wcheck in WVEC.
1233 : !
1234 432 : DO k = 1, npt
1235 : tempa = zero
1236 : tempb = zero
1237 : tempc = zero
1238 1080 : DO i = 1, n
1239 720 : tempa = tempa + xpt(k, i)*d(i)
1240 720 : tempb = tempb + xpt(k, i)*s(i)
1241 1080 : tempc = tempc + xpt(k, i)*xopt(i)
1242 : END DO
1243 360 : wvec(k, 1) = quart*(tempa*tempa + tempb*tempb)
1244 360 : wvec(k, 2) = tempa*tempc
1245 360 : wvec(k, 3) = tempb*tempc
1246 360 : wvec(k, 4) = quart*(tempa*tempa - tempb*tempb)
1247 432 : wvec(k, 5) = half*tempa*tempb
1248 : END DO
1249 216 : DO i = 1, n
1250 144 : ip = i + npt
1251 144 : wvec(ip, 1) = zero
1252 144 : wvec(ip, 2) = d(i)
1253 144 : wvec(ip, 3) = s(i)
1254 144 : wvec(ip, 4) = zero
1255 216 : wvec(ip, 5) = zero
1256 : END DO
1257 : !
1258 : ! Put the coefficients of THETA*Wcheck in PROD.
1259 : !
1260 432 : DO jc = 1, 5
1261 360 : nw = npt
1262 360 : IF (jc == 2 .OR. jc == 3) nw = ndim
1263 2160 : DO k = 1, npt
1264 2160 : prod(k, jc) = zero
1265 : END DO
1266 1080 : DO j = 1, nptm
1267 : sum = zero
1268 4320 : DO k = 1, npt
1269 4320 : sum = sum + zmat(k, j)*wvec(k, jc)
1270 : END DO
1271 720 : IF (j < idz) sum = -sum
1272 4680 : DO k = 1, npt
1273 4320 : prod(k, jc) = prod(k, jc) + sum*zmat(k, j)
1274 : END DO
1275 : END DO
1276 360 : IF (nw == ndim) THEN
1277 864 : DO k = 1, npt
1278 : sum = zero
1279 2160 : DO j = 1, n
1280 2160 : sum = sum + bmat(k, j)*wvec(npt + j, jc)
1281 : END DO
1282 864 : prod(k, jc) = prod(k, jc) + sum
1283 : END DO
1284 : END IF
1285 1152 : DO j = 1, n
1286 : sum = zero
1287 4896 : DO i = 1, nw
1288 4896 : sum = sum + bmat(i, j)*wvec(i, jc)
1289 : END DO
1290 1080 : 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 576 : DO k = 1, ndim
1297 : sum = zero
1298 3024 : DO I = 1, 5
1299 2520 : par(i) = half*prod(k, i)*wvec(k, i)
1300 3024 : sum = sum + par(i)
1301 : END DO
1302 504 : den(1) = den(1) - par(1) - sum
1303 504 : tempa = prod(k, 1)*wvec(k, 2) + prod(k, 2)*wvec(k, 1)
1304 504 : tempb = prod(k, 2)*wvec(k, 4) + prod(k, 4)*wvec(k, 2)
1305 504 : tempc = prod(k, 3)*wvec(k, 5) + prod(k, 5)*wvec(k, 3)
1306 504 : den(2) = den(2) - tempa - half*(tempb + tempc)
1307 504 : den(6) = den(6) - half*(tempb - tempc)
1308 504 : tempa = prod(k, 1)*wvec(k, 3) + prod(k, 3)*wvec(k, 1)
1309 504 : tempb = prod(k, 2)*wvec(k, 5) + prod(k, 5)*wvec(k, 2)
1310 504 : tempc = prod(k, 3)*wvec(k, 4) + prod(k, 4)*wvec(k, 3)
1311 504 : den(3) = den(3) - tempa - half*(tempb - tempc)
1312 504 : den(7) = den(7) - half*(tempb + tempc)
1313 504 : tempa = prod(k, 1)*wvec(k, 4) + prod(k, 4)*wvec(k, 1)
1314 504 : den(4) = den(4) - tempa - par(2) + par(3)
1315 504 : tempa = prod(k, 1)*wvec(k, 5) + prod(k, 5)*wvec(k, 1)
1316 504 : tempb = prod(k, 2)*wvec(k, 3) + prod(k, 3)*wvec(k, 2)
1317 504 : den(5) = den(5) - tempa - half*tempb
1318 504 : den(8) = den(8) - par(4) + par(5)
1319 504 : tempa = prod(k, 4)*wvec(k, 5) + prod(k, 5)*wvec(k, 4)
1320 576 : 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 432 : DO i = 1, 5
1327 360 : par(i) = half*prod(knew, i)**2
1328 432 : sum = sum + par(i)
1329 : END DO
1330 72 : denex(1) = alpha*den(1) + par(1) + sum
1331 72 : tempa = two*prod(knew, 1)*prod(knew, 2)
1332 72 : tempb = prod(knew, 2)*prod(knew, 4)
1333 72 : tempc = prod(knew, 3)*prod(knew, 5)
1334 72 : denex(2) = alpha*den(2) + tempa + tempb + tempc
1335 72 : denex(6) = alpha*den(6) + tempb - tempc
1336 72 : tempa = two*prod(knew, 1)*prod(knew, 3)
1337 72 : tempb = prod(knew, 2)*prod(knew, 5)
1338 72 : tempc = prod(knew, 3)*prod(knew, 4)
1339 72 : denex(3) = alpha*den(3) + tempa + tempb - tempc
1340 72 : denex(7) = alpha*den(7) + tempb + tempc
1341 72 : tempa = two*prod(knew, 1)*prod(knew, 4)
1342 72 : denex(4) = alpha*den(4) + tempa + par(2) - par(3)
1343 72 : tempa = two*prod(knew, 1)*prod(knew, 5)
1344 72 : denex(5) = alpha*den(5) + tempa + prod(knew, 2)*prod(knew, 3)
1345 72 : denex(8) = alpha*den(8) + par(4) - par(5)
1346 72 : 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 72 : sum = denex(1) + denex(2) + denex(4) + denex(6) + denex(8)
1351 72 : denold = sum
1352 72 : denmax = sum
1353 72 : isave = 0
1354 72 : iu = 49
1355 72 : temp = twopi/REAL(IU + 1, dp)
1356 72 : par(1) = one
1357 3600 : DO i = 1, iu
1358 3528 : angle = REAL(i, dp)*temp
1359 3528 : par(2) = COS(angle)
1360 3528 : par(3) = SIN(angle)
1361 3528 : DO j = 4, 8, 2
1362 10584 : par(j) = par(2)*par(j - 2) - par(3)*par(j - 1)
1363 10584 : par(j + 1) = par(2)*par(j - 1) + par(3)*par(j - 2)
1364 : END DO
1365 : sumold = sum
1366 : sum = zero
1367 35280 : DO j = 1, 9
1368 35280 : sum = sum + denex(j)*par(j)
1369 : END DO
1370 3600 : IF (ABS(sum) > ABS(denmax)) THEN
1371 : denmax = sum
1372 : isave = i
1373 : tempa = sumold
1374 3512 : ELSE IF (i == isave + 1) THEN
1375 88 : tempb = sum
1376 : END IF
1377 : END DO
1378 72 : IF (isave == 0) tempa = sum
1379 16 : IF (isave == iu) tempb = denold
1380 72 : step = zero
1381 72 : IF (tempa /= tempb) THEN
1382 72 : tempa = tempa - denmax
1383 72 : tempb = tempb - denmax
1384 72 : step = half*(tempa - tempb)/(tempa + tempb)
1385 : END IF
1386 72 : 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 72 : par(2) = COS(angle)
1392 72 : par(3) = SIN(angle)
1393 72 : DO j = 4, 8, 2
1394 216 : par(j) = par(2)*par(j - 2) - par(3)*par(j - 1)
1395 216 : par(j + 1) = par(2)*par(j - 1) + par(3)*par(j - 2)
1396 : END DO
1397 72 : beta = zero
1398 72 : denmax = zero
1399 720 : DO j = 1, 9
1400 648 : beta = beta + den(j)*par(j)
1401 720 : denmax = denmax + denex(j)*par(j)
1402 : END DO
1403 576 : DO k = 1, ndim
1404 504 : vlag(k) = zero
1405 3096 : DO j = 1, 5
1406 3024 : vlag(k) = vlag(k) + prod(k, j)*par(j)
1407 : END DO
1408 : END DO
1409 72 : tau = vlag(knew)
1410 72 : dd = zero
1411 72 : tempa = zero
1412 72 : tempb = zero
1413 216 : DO i = 1, n
1414 144 : d(i) = par(2)*d(i) + par(3)*s(i)
1415 144 : w(i) = xopt(i) + d(i)
1416 144 : dd = dd + d(i)**2
1417 144 : tempa = tempa + d(i)*w(i)
1418 216 : tempb = tempb + w(i)*w(i)
1419 : END DO
1420 72 : IF (iterc >= n) EXIT mainloop
1421 72 : IF (iterc >= 1) densav = MAX(densav, denold)
1422 72 : IF (ABS(denmax) <= 1.1_dp*ABS(densav)) EXIT mainloop
1423 0 : 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 0 : DO i = 1, n
1429 0 : temp = tempa*xopt(i) + tempb*d(i) - vlag(npt + i)
1430 0 : s(i) = tau*bmat(knew, i) + alpha*temp
1431 : END DO
1432 0 : DO k = 1, npt
1433 : sum = zero
1434 0 : DO j = 1, n
1435 0 : sum = sum + xpt(k, j)*w(j)
1436 : END DO
1437 0 : temp = (tau*w(n + k) - alpha*vlag(k))*sum
1438 0 : DO i = 1, n
1439 0 : s(i) = s(i) + temp*xpt(k, i)
1440 : END DO
1441 : END DO
1442 : ss = zero
1443 : ds = zero
1444 0 : DO i = 1, n
1445 0 : ss = ss + s(i)**2
1446 0 : ds = ds + d(i)*s(i)
1447 : END DO
1448 0 : ssden = dd*ss - ds*ds
1449 72 : 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 576 : DO k = 1, ndim
1455 504 : w(k) = zero
1456 3096 : DO j = 1, 5
1457 3024 : w(k) = w(k) + wvec(k, j)*par(j)
1458 : END DO
1459 : END DO
1460 72 : vlag(kopt) = vlag(kopt) + one
1461 :
1462 72 : 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 20201950 : 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 20201950 : delsq = delta*delta
1527 20201950 : 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 20201950 : iterc = 0
1533 121217904 : DO k = 1, npt
1534 121217904 : hcol(k) = zero
1535 : END DO
1536 60608952 : DO j = 1, nptm
1537 40407002 : temp = zmat(knew, j)
1538 40407002 : IF (j < idz) temp = -temp
1539 262692590 : DO k = 1, npt
1540 242490640 : hcol(k) = hcol(k) + temp*zmat(k, j)
1541 : END DO
1542 : END DO
1543 20201950 : 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 20201950 : dd = zero
1549 60608952 : DO i = 1, n
1550 40407002 : d(i) = xpt(knew, i) - xopt(i)
1551 40407002 : gc(i) = bmat(knew, i)
1552 40407002 : gd(i) = zero
1553 60608952 : dd = dd + d(i)**2
1554 : END DO
1555 121217904 : DO k = 1, npt
1556 : temp = zero
1557 : sum = zero
1558 303099592 : DO j = 1, n
1559 202083638 : temp = temp + xpt(k, j)*xopt(j)
1560 303099592 : sum = sum + xpt(k, j)*d(j)
1561 : END DO
1562 101015954 : temp = hcol(k)*temp
1563 101015954 : sum = hcol(k)*sum
1564 323301542 : DO i = 1, n
1565 202083638 : gc(i) = gc(i) + temp*xpt(k, i)
1566 303099592 : 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 60608952 : DO i = 1, n
1577 40407002 : gg = gg + gc(i)**2
1578 40407002 : sp = sp + d(i)*gc(i)
1579 60608952 : dhd = dhd + d(i)*gd(i)
1580 : END DO
1581 20201950 : scale = delta/SQRT(dd)
1582 20201950 : IF (sp*dhd < zero) scale = -scale
1583 20201950 : temp = zero
1584 20201950 : IF (sp*sp > 0.99_dp*dd*gg) temp = one
1585 20201950 : tau = scale*(ABS(sp) + half*scale*ABS(dhd))
1586 20201950 : IF (gg*delsq < 0.01_dp*tau*tau) temp = one
1587 60608952 : DO i = 1, n
1588 40407002 : d(i) = scale*d(i)
1589 40407002 : gd(i) = scale*gd(i)
1590 60608952 : 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 28583397 : iterc = iterc + 1
1599 28583397 : dd = zero
1600 28583397 : sp = zero
1601 28583397 : ss = zero
1602 85755515 : DO i = 1, n
1603 57172118 : dd = dd + d(i)**2
1604 57172118 : sp = sp + d(i)*s(i)
1605 85755515 : ss = ss + s(i)**2
1606 : END DO
1607 28583397 : temp = dd*ss - sp*sp
1608 28583397 : IF (temp <= 1.0e-8_dp*dd*ss) EXIT mainloop
1609 27358353 : denom = SQRT(temp)
1610 82080109 : DO i = 1, n
1611 54721756 : s(i) = (dd*s(i) - sp*d(i))/denom
1612 82080109 : 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 164160218 : DO k = 1, npt
1619 : sum = zero
1620 410498545 : DO j = 1, n
1621 410498545 : sum = sum + xpt(k, j)*s(j)
1622 : END DO
1623 136801865 : sum = hcol(k)*sum
1624 437856898 : DO i = 1, n
1625 410498545 : 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 82080109 : DO i = 1, n
1634 54721756 : cf1 = cf1 + s(i)*w(i)
1635 54721756 : cf2 = cf2 + d(i)*gc(i)
1636 54721756 : cf3 = cf3 + s(i)*gc(i)
1637 54721756 : cf4 = cf4 + d(i)*gd(i)
1638 82080109 : cf5 = cf5 + s(i)*gd(i)
1639 : END DO
1640 27358353 : cf1 = half*cf1
1641 27358353 : cf4 = half*cf4 - cf1
1642 : !
1643 : ! Seek the value of the angle that maximizes the modulus of TAU.
1644 : !
1645 27358353 : taubeg = cf1 + cf2 + cf4
1646 27358353 : taumax = taubeg
1647 27358353 : tauold = taubeg
1648 27358353 : isave = 0
1649 27358353 : iu = 49
1650 27358353 : temp = twopi/REAL(iu + 1, DP)
1651 1367917650 : DO i = 1, iu
1652 1340559297 : angle = REAL(i, dp)*temp
1653 1340559297 : cth = COS(angle)
1654 1340559297 : sth = SIN(angle)
1655 1340559297 : tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1656 1340559297 : IF (ABS(tau) > ABS(taumax)) THEN
1657 : taumax = tau
1658 : isave = i
1659 : tempa = tauold
1660 1270830995 : ELSE IF (i == isave + 1) THEN
1661 29724085 : tempb = taU
1662 : END IF
1663 1367917650 : tauold = tau
1664 : END DO
1665 27358353 : IF (isave == 0) tempa = tau
1666 16050116 : IF (isave == iu) tempb = taubeg
1667 27358353 : step = zero
1668 27358353 : IF (tempa /= tempb) THEN
1669 27358353 : tempa = tempa - taumax
1670 27358353 : tempb = tempb - taumax
1671 27358353 : step = half*(tempa - tempb)/(tempa + tempb)
1672 : END IF
1673 27358353 : angle = temp*(REAL(isave, DP) + step)
1674 : !
1675 : ! Calculate the new D and GD. Then test for convergence.
1676 : !
1677 27358353 : cth = COS(angle)
1678 27358353 : sth = SIN(angle)
1679 27358353 : tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1680 82080109 : DO i = 1, n
1681 54721756 : d(i) = cth*d(i) + sth*s(i)
1682 54721756 : gd(i) = cth*gd(i) + sth*w(i)
1683 82080109 : s(i) = gc(i) + gd(i)
1684 : END DO
1685 27358353 : IF (ABS(tau) <= 1.1_dp*ABS(taubeg)) EXIT mainloop
1686 28583397 : IF (iterc >= n) EXIT mainloop
1687 : END DO mainloop
1688 :
1689 20201950 : 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 48512844 : 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 48512844 : delsq = delta*delta
1749 48512844 : iterc = 0
1750 48512844 : itermax = n
1751 48512844 : itersw = itermax
1752 145544523 : DO i = 1, n
1753 145544523 : d(i) = xopt(i)
1754 : END DO
1755 48512844 : CALL updatehd
1756 : !
1757 : ! Prepare for the first line search.
1758 : !
1759 48512844 : qred = zero
1760 48512844 : dd = zero
1761 145544523 : DO i = 1, n
1762 97031679 : step(i) = zero
1763 97031679 : hs(i) = zero
1764 97031679 : g(i) = gq(i) + hd(i)
1765 97031679 : d(i) = -g(i)
1766 145544523 : dd = dd + d(i)**2
1767 : END DO
1768 48512844 : crvmin = zero
1769 48512844 : 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 81997759 : IF (.NOT. jump2) THEN
1781 81996565 : IF (.NOT. jump1) THEN
1782 81996565 : iterc = iterc + 1
1783 81996565 : temp = delsq - ss
1784 81996565 : bstep = temp/(ds + SQRT(ds*ds + dd*temp))
1785 81996565 : CALL updatehd
1786 : END IF
1787 81996565 : jump1 = .FALSE.
1788 81996565 : IF (iterc <= itersw) THEN
1789 81996565 : dhd = zero
1790 246005935 : DO j = 1, n
1791 246005935 : dhd = dhd + d(j)*hd(j)
1792 : END DO
1793 : !
1794 : ! Update CRVMIN and set the step-length ALPHA.
1795 : !
1796 81996565 : alpha = bstep
1797 81996565 : IF (dhd > zero) THEN
1798 72438003 : temp = dhd/dd
1799 72438003 : IF (iterc == 1) crvmin = temp
1800 72438003 : crvmin = MIN(crvmin, temp)
1801 72438003 : alpha = MIN(alpha, gg/dhd)
1802 : END IF
1803 81996565 : qadd = alpha*(gg - half*alpha*dhd)
1804 81996565 : qred = qred + qadd
1805 : !
1806 : ! Update STEP and HS.
1807 : !
1808 81996565 : ggsav = gg
1809 81996565 : gg = zero
1810 246005935 : DO i = 1, n
1811 164009370 : step(i) = step(i) + alpha*d(i)
1812 164009370 : hs(i) = hs(i) + alpha*hd(i)
1813 246005935 : gg = gg + (g(i) + hs(i))**2
1814 : END DO
1815 : !
1816 : ! Begin another conjugate direction iteration if required.
1817 : !
1818 81996565 : IF (alpha < bstep) THEN
1819 53992835 : IF (qadd <= 0.01_dp*qred) EXIT mainloop
1820 53179355 : IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1821 33484074 : IF (iterc == itermax) EXIT mainloop
1822 33483874 : temp = gg/ggsav
1823 33483874 : dd = zero
1824 33483874 : ds = zero
1825 33483874 : ss = zero
1826 100462202 : DO i = 1, n
1827 66978328 : d(i) = temp*d(i) - g(i) - hs(i)
1828 66978328 : dd = dd + d(i)**2
1829 66978328 : ds = ds + d(i)*step(I)
1830 100462202 : ss = ss + step(i)**2
1831 : END DO
1832 33483874 : IF (ds <= zero) EXIT mainloop
1833 33483874 : IF (ss < delsq) CYCLE mainloop
1834 : END IF
1835 28003730 : crvmin = zero
1836 28003730 : itersw = iterc
1837 28003730 : jump2 = .TRUE.
1838 28003730 : 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 27584300 : sg = zero
1849 27584300 : shs = zero
1850 82764729 : DO i = 1, n
1851 55180429 : sg = sg + step(i)*g(i)
1852 82764729 : shs = shs + step(i)*hs(i)
1853 : END DO
1854 27584300 : sgk = sg + shs
1855 27584300 : angtest = sgk/SQRT(gg*delsq)
1856 27584300 : 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 24083468 : iterc = iterc + 1
1862 24083468 : temp = SQRT(delsq*gg - sgk*sgk)
1863 24083468 : tempa = delsq/temp
1864 24083468 : tempb = sgk/temp
1865 72261505 : DO i = 1, n
1866 72261505 : d(i) = tempa*(g(i) + hs(i)) - tempb*step(i)
1867 : END DO
1868 24083468 : CALL updatehd
1869 24083468 : IF (iterc <= itersw) THEN
1870 : jump1 = .TRUE.
1871 : CYCLE mainloop
1872 : END IF
1873 : END IF
1874 24083468 : dg = zero
1875 24083468 : dhd = zero
1876 24083468 : dhs = zero
1877 72261505 : DO i = 1, n
1878 48178037 : dg = dg + d(i)*g(i)
1879 48178037 : dhd = dhd + hd(i)*d(i)
1880 72261505 : dhs = dhs + hd(i)*step(i)
1881 : END DO
1882 : !
1883 : ! Seek the value of the angle that minimizes Q.
1884 : !
1885 24083468 : cf = half*(shs - dhd)
1886 24083468 : qbeg = sg + cf
1887 24083468 : qsav = qbeg
1888 24083468 : qmin = qbeg
1889 24083468 : isave = 0
1890 24083468 : iu = 49
1891 : temp = twopi/REAL(iu + 1, dp)
1892 1204173400 : DO i = 1, iu
1893 1180089932 : angle = REAL(i, dp)*temp
1894 1180089932 : cth = COS(angle)
1895 1180089932 : sth = SIN(angle)
1896 1180089932 : qnew = (sg + cf*cth)*cth + (dg + dhs*cth)*sth
1897 1180089932 : IF (qnew < qmin) THEN
1898 : qmin = qnew
1899 : isave = i
1900 : tempa = qsav
1901 1154724168 : ELSE IF (i == isave + 1) THEN
1902 30123085 : tempb = qnew
1903 : END IF
1904 1204173400 : qsav = qnew
1905 : END DO
1906 24083468 : IF (isave == zero) tempa = qnew
1907 10546968 : IF (isave == iu) tempb = qbeg
1908 24083468 : angle = zero
1909 24083468 : IF (tempa /= tempb) THEN
1910 24083468 : tempa = tempa - qmin
1911 24083468 : tempb = tempb - qmin
1912 24083468 : angle = half*(tempa - tempb)/(tempa + tempb)
1913 : END IF
1914 24083468 : angle = temp*(REAL(isave, DP) + angle)
1915 : !
1916 : ! Calculate the new STEP and HS. Then test for convergence.
1917 : !
1918 24083468 : cth = COS(angle)
1919 24083468 : sth = SIN(angle)
1920 24083468 : reduc = qbeg - (sg + cf*cth)*cth - (dg + dhs*cth)*sth
1921 24083468 : gg = zero
1922 72261505 : DO i = 1, n
1923 48178037 : step(i) = cth*step(i) + sth*d(i)
1924 48178037 : hs(i) = cth*hs(i) + sth*hd(i)
1925 72261505 : gg = gg + (g(i) + hs(i))**2
1926 : END DO
1927 24083468 : qred = qred + reduc
1928 24083468 : ratio = reduc/qred
1929 24083468 : IF (iterc < itermax .AND. ratio > 0.01_dp) THEN
1930 1194 : jump2 = .TRUE.
1931 : ELSE
1932 : EXIT mainloop
1933 : END IF
1934 :
1935 48513885 : IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1936 :
1937 : END DO mainloop
1938 :
1939 : CONTAINS
1940 : ! **************************************************************************************************
1941 : !> \brief ...
1942 : ! **************************************************************************************************
1943 154592877 : SUBROUTINE updatehd
1944 : INTEGER :: i, ih, j, k
1945 :
1946 463811963 : DO i = 1, n
1947 463811963 : hd(i) = zero
1948 : END DO
1949 927623926 : DO k = 1, npt
1950 773031049 : temp = zero
1951 2319713107 : DO j = 1, n
1952 2319713107 : temp = temp + xpt(k, j)*d(j)
1953 : END DO
1954 773031049 : temp = temp*pq(k)
1955 2474305984 : DO i = 1, n
1956 2319713107 : hd(i) = hd(i) + temp*xpt(k, i)
1957 : END DO
1958 : END DO
1959 154592877 : ih = 0
1960 463811963 : DO j = 1, n
1961 927787249 : DO i = 1, j
1962 463975286 : ih = ih + 1
1963 463975286 : IF (i < j) hd(j) = hd(j) + hq(ih)*d(i)
1964 773194372 : hd(i) = hd(i) + hq(ih)*d(j)
1965 : END DO
1966 : END DO
1967 154592877 : 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 56332932 : 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 56332932 : nptm = npt - n - 1
2009 : !
2010 : ! Apply the rotations that put zeros in the KNEW-th row of ZMAT.
2011 : !
2012 56332932 : jl = 1
2013 112674063 : DO j = 2, nptm
2014 112674063 : IF (j == idz) THEN
2015 : jl = idz
2016 56341131 : ELSE IF (zmat(knew, j) /= zero) THEN
2017 54804156 : temp = SQRT(zmat(knew, jl)**2 + zmat(knew, j)**2)
2018 54804156 : tempa = zmat(knew, jl)/temp
2019 54804156 : tempb = zmat(knew, j)/temp
2020 328935492 : DO I = 1, NPT
2021 274131336 : temp = tempa*zmat(i, jl) + tempb*zmat(i, j)
2022 274131336 : zmat(i, j) = tempa*zmat(i, j) - tempb*zmat(i, jl)
2023 328935492 : zmat(i, jl) = temp
2024 : END DO
2025 54804156 : 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 56332932 : tempa = zmat(knew, 1)
2033 56332932 : IF (idz >= 2) tempa = -tempa
2034 56332932 : IF (jl > 1) tempb = zmat(knew, jl)
2035 338013990 : DO i = 1, npt
2036 281681058 : w(i) = tempa*zmat(i, 1)
2037 338013990 : IF (jl > 1) w(i) = w(i) + tempb*zmat(i, jl)
2038 : END DO
2039 56332932 : alpha = w(knew)
2040 56332932 : tau = vlag(knew)
2041 56332932 : tausq = tau*tau
2042 56332932 : denom = alpha*beta + tausq
2043 56332932 : 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 56332932 : iflag = 0
2050 56332932 : IF (JL == 1) THEN
2051 56332932 : temp = SQRT(ABS(denom))
2052 56332932 : tempb = tempa/temp
2053 56332932 : tempa = tau/temp
2054 338013990 : DO i = 1, npt
2055 338013990 : zmat(i, 1) = tempa*zmat(i, 1) - tempb*vlag(i)
2056 : END DO
2057 56332932 : IF (idz == 1 .AND. temp < zero) idz = 2
2058 56332932 : 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 169006995 : DO j = 1, n
2097 112674063 : jp = npt + j
2098 112674063 : w(jp) = bmat(knew, j)
2099 112674063 : tempa = (alpha*vlag(jp) - tau*w(jp))/denom
2100 112674063 : tempb = (-beta*w(jp) - tau*vlag(jp))/denom
2101 901560977 : DO i = 1, jp
2102 732553982 : bmat(i, j) = bmat(i, j) + tempa*vlag(i) + tempb*w(i)
2103 845228045 : IF (i > npt) bmat(jp, i - npt) = bmat(i, j)
2104 : END DO
2105 : END DO
2106 :
2107 56332932 : END SUBROUTINE update
2108 :
2109 0 : END MODULE powell
|
