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 61126177 : 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 61126177 : CALL timeset(routineN, handle)
61 :
62 61900036 : SELECT CASE (optstate%state)
63 : CASE (0)
64 773859 : npt = 2*n + 1
65 2321577 : ALLOCATE (optstate%w((npt + 13)*(npt + n) + 3*n*(n + 3)/2))
66 2321577 : ALLOCATE (optstate%xopt(n))
67 : ! Initialize w
68 110032107 : optstate%w = 0.0_dp
69 773859 : optstate%state = 1
70 773859 : CALL newuoa(n, x, optstate)
71 : CASE (1, 2)
72 58804605 : 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 773841 : optstate%state = -1
90 : CASE (8)
91 2321931 : x(1:n) = optstate%xopt(1:n)
92 773859 : DEALLOCATE (optstate%w)
93 773859 : DEALLOCATE (optstate%xopt)
94 773859 : optstate%state = -1
95 : CASE DEFAULT
96 61126177 : CPABORT("Unknown optimization state")
97 : END SELECT
98 :
99 61126177 : CALL timestop(handle)
100 :
101 61126177 : END SUBROUTINE powell_optimize
102 : ! **************************************************************************************************
103 : !> \brief ...
104 : !> \param n ...
105 : !> \param x ...
106 : !> \param optstate ...
107 : ! **************************************************************************************************
108 59578464 : 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 59578464 : maxfun = optstate%maxfun
120 59578464 : rhobeg = optstate%rhobeg
121 59578464 : 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 59578464 : np = n + 1
158 59578464 : npt = 2*n + 1
159 59578464 : nptm = npt - np
160 59578464 : IF (npt < n + 2 .OR. npt > ((n + 2)*np)/2) THEN
161 0 : optstate%state = 5
162 0 : RETURN
163 : END IF
164 59578464 : ndim = npt + n
165 59578464 : ixb = 1
166 59578464 : ixo = ixb + n
167 59578464 : ixn = ixo + n
168 59578464 : ixp = ixn + n
169 59578464 : ifv = ixp + n*npt
170 59578464 : igq = ifv + npt
171 59578464 : ihq = igq + n
172 59578464 : ipq = ihq + (n*np)/2
173 59578464 : ibmat = ipq + npt
174 59578464 : izmat = ibmat + ndim*n
175 59578464 : id = izmat + npt*nptm
176 59578464 : ivl = id + n
177 59578464 : 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 59578464 : optstate%w(ivl:), optstate%w(iw:), optstate)
187 :
188 178750260 : 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 59578464 : SUBROUTINE newuob(n, npt, x, rhobeg, rhoend, maxfun, xbase, &
217 59578464 : 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 59578464 : skip_check = .FALSE.
273 59578464 : do_return = .FALSE.
274 :
275 59578464 : 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 773859 : idz = 0
282 773859 : itest = 0
283 773859 : nf = 0
284 773859 : nfm = 0
285 773859 : nfmm = 0
286 773859 : nfsav = 0
287 773859 : knew = 0
288 773859 : kopt = 0
289 773859 : ksave = 0
290 773859 : ktemp = 0
291 773859 : rhosq = 0._dp
292 773859 : recip = 0._dp
293 773859 : reciq = 0._dp
294 773859 : fbeg = 0._dp
295 773859 : fopt = 0._dp
296 773859 : diffa = 0._dp
297 773859 : xoptsq = 0._dp
298 773859 : rho = 0._dp
299 773859 : delta = 0._dp
300 773859 : dsq = 0._dp
301 773859 : dnorm = 0._dp
302 773859 : ratio = 0._dp
303 773859 : temp = 0._dp
304 773859 : tempq = 0._dp
305 773859 : beta = 0._dp
306 773859 : dx = 0._dp
307 773859 : vquad = 0._dp
308 773859 : diff = 0._dp
309 773859 : diffc = 0._dp
310 773859 : diffb = 0._dp
311 773859 : fsave = 0._dp
312 773859 : detrat = 0._dp
313 773859 : hdiag = 0._dp
314 773859 : distsq = 0._dp
315 773859 : gisq = 0._dp
316 773859 : gqsq = 0._dp
317 773859 : f = 0._dp
318 773859 : bstep = 0._dp
319 773859 : alpha = 0._dp
320 773859 : dstep = 0._dp
321 : !
322 : END IF
323 :
324 59578464 : ipt = 0
325 59578464 : jpt = 0
326 59578464 : xipt = 0._dp
327 59578464 : xjpt = 0._dp
328 :
329 59578464 : half = 0.5_dp
330 59578464 : one = 1.0_dp
331 59578464 : tenth = 0.1_dp
332 59578464 : zero = 0.0_dp
333 59578464 : np = n + 1
334 59578464 : nh = (n*np)/2
335 59578464 : nptm = npt - np
336 59578464 : nftest = MAX(maxfun, 1)
337 :
338 59578464 : IF (opt%state /= 2) THEN
339 : !
340 : ! Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero.
341 : !
342 2321931 : DO j = 1, n
343 1548072 : xbase(j) = x(j)
344 9320872 : DO k = 1, npt
345 9320872 : xpt(k, j) = zero
346 : END DO
347 13207095 : DO i = 1, ndim
348 12433236 : bmat(i, j) = zero
349 : END DO
350 : END DO
351 3104077 : DO ih = 1, nh
352 3104077 : hq(ih) = zero
353 : END DO
354 4643862 : DO k = 1, npt
355 3870003 : pq(k) = zero
356 12416662 : DO j = 1, nptm
357 11642803 : 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 773859 : rhosq = rhobeg*rhobeg
366 773859 : recip = one/rhosq
367 773859 : reciq = SQRT(half)/rhosq
368 773859 : nf = 0
369 773859 : CALL begin_initialization_of_interpolation
370 24484913 : RETURN
371 : END IF
372 :
373 58804605 : CALL set_state
374 :
375 58804605 : preparation: IF (nf <= npt) THEN
376 3869526 : fval(nf) = f
377 3869526 : IF (nf == 1) THEN
378 773859 : fbeg = f
379 773859 : fopt = f
380 773859 : kopt = 1
381 3095667 : ELSE IF (f < fopt) THEN
382 864236 : fopt = f
383 864236 : 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 3869526 : IF (NFM <= 2*N) THEN
390 3869526 : IF (nfm >= 1 .AND. nfm <= n) THEN
391 1547866 : gq(nfm) = (f - fbeg)/rhobeg
392 1547866 : 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 2321660 : ELSE IF (nfm > n) THEN
398 1547801 : bmat(nf - n, nfmm) = half/rhobeg
399 1547801 : bmat(nf, nfmm) = -half/rhobeg
400 1547801 : zmat(1, nfmm) = -reciq - reciq
401 1547801 : zmat(nf - n, nfmm) = reciq
402 1547801 : zmat(nf, nfmm) = reciq
403 1547801 : ih = (nfmm*(nfmm + 1))/2
404 1547801 : temp = (fbeg - f)/rhobeg
405 1547801 : hq(ih) = (gq(nfmm) - temp)/rhobeg
406 1547801 : 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 3869526 : IF (nf < npt) THEN
423 3095676 : CALL begin_initialization_of_interpolation
424 3095676 : RETURN
425 : END IF
426 : !
427 : ! Begin the iterative procedure, because the initial model is comple
428 : !
429 773850 : rho = rhobeg
430 773850 : delta = rho
431 773850 : idz = 1
432 773850 : diffa = zero
433 773850 : diffb = zero
434 773850 : itest = 0
435 773850 : xoptsq = zero
436 2321651 : DO i = 1, n
437 1547801 : xopt(i) = xpt(kopt, i)
438 2321651 : xoptsq = xoptsq + xopt(i)**2
439 : END DO
440 773850 : 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 773850 : skip_check = .TRUE.
446 : ELSE preparation
447 54935079 : IF (knew == -1) THEN
448 645970 : opt%state = 6
449 645970 : CALL get_state
450 645970 : CALL finalize_optimization_after_failure
451 645970 : 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 54289109 : vquad = zero
458 54289109 : ih = 0
459 162875196 : DO j = 1, n
460 108586087 : vquad = vquad + d(j)*gq(j)
461 325787950 : DO i = 1, j
462 162912754 : ih = ih + 1
463 162912754 : temp = d(i)*xnew(j) + d(j)*xopt(i)
464 162912754 : IF (i == j) temp = half*temp
465 271498841 : vquad = vquad + temp*hq(ih)
466 : END DO
467 : END DO
468 325750392 : DO k = 1, npt
469 325750392 : vquad = vquad + pq(k)*w(k)
470 : END DO
471 54289109 : diff = f - fopt - vquad
472 54289109 : diffc = diffb
473 54289109 : diffb = diffa
474 54289109 : diffa = ABS(diff)
475 54289109 : 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 54289109 : fsave = fopt
482 54289109 : IF (f < fopt) THEN
483 27742324 : fopt = f
484 27742324 : xoptsq = zero
485 83230235 : DO i = 1, n
486 55487911 : xopt(i) = xnew(i)
487 83230235 : xoptsq = xoptsq + xopt(i)**2
488 : END DO
489 : END IF
490 54289109 : ksave = knew
491 54289109 : IF (knew <= 0) THEN
492 : !
493 : ! Pick the next value of DELTA after a trust region step.
494 : !
495 35093549 : 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 35093545 : ratio = (f - fsave)/vquad
503 35093545 : IF (ratio <= tenth) THEN
504 12992920 : delta = half*dnorm
505 22100625 : ELSE IF (ratio <= 0.7_dp) THEN
506 3820981 : delta = MAX(half*delta, dnorm)
507 : ELSE
508 18279644 : delta = MAX(half*delta, dnorm + dnorm)
509 : END IF
510 35093545 : 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 35093545 : rhosq = MAX(tenth*delta, rho)**2
515 35093545 : ktemp = 0
516 35093545 : detrat = zero
517 35093545 : IF (f >= fsave) THEN
518 11154722 : ktemp = kopt
519 11154722 : detrat = one
520 : END IF
521 210571132 : DO k = 1, npt
522 175477587 : hdiag = zero
523 526525818 : DO j = 1, nptm
524 351048231 : temp = one
525 351048231 : IF (j < idz) temp = -one
526 526525818 : hdiag = hdiag + temp*zmat(k, j)**2
527 : END DO
528 175477587 : temp = ABS(beta*hdiag + vlag(k)**2)
529 175477587 : distsq = zero
530 526525818 : DO j = 1, n
531 526525818 : distsq = distsq + (xpt(k, j) - xopt(j))**2
532 : END DO
533 175477587 : IF (distsq > rhosq) temp = temp*(distsq/rhosq)**3
534 210571132 : IF (temp > detrat .AND. k /= ktemp) THEN
535 75565612 : detrat = temp
536 75565612 : knew = k
537 : END IF
538 : END DO
539 35093545 : 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 53830933 : CALL update(n, npt, bmat, zmat, idz, ndim, vlag, beta, knew, w)
547 53830933 : fval(knew) = f
548 53830933 : ih = 0
549 161500614 : DO i = 1, n
550 107669681 : temp = pq(knew)*xpt(knew, i)
551 323038612 : DO j = 1, i
552 161537998 : ih = ih + 1
553 269207679 : hq(ih) = hq(ih) + temp*xpt(knew, j)
554 : END DO
555 : END DO
556 53830933 : 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 161500614 : DO j = 1, nptm
562 107669681 : temp = diff*zmat(knew, j)
563 107669681 : IF (j < idz) temp = -temp
564 699982925 : DO k = 1, npt
565 646151992 : pq(k) = pq(k) + temp*zmat(k, j)
566 : END DO
567 : END DO
568 53830933 : gqsq = zero
569 161500614 : DO i = 1, n
570 107669681 : gq(i) = gq(i) + diff*bmat(knew, i)
571 107669681 : gqsq = gqsq + gq(i)**2
572 161500614 : 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 53830933 : IF (ksave == 0 .AND. delta == rho) THEN
580 6958117 : IF (ABS(ratio) > 1.0e-2_dp) THEN
581 4679773 : itest = 0
582 : ELSE
583 13670138 : DO k = 1, npt
584 13670138 : vlag(k) = fval(k) - fval(kopt)
585 : END DO
586 2278344 : gisq = zero
587 6835069 : DO i = 1, n
588 4556725 : sum = zero
589 27341140 : DO k = 1, npt
590 27341140 : sum = sum + bmat(k, i)*vlag(k)
591 : END DO
592 4556725 : gisq = gisq + sum*sum
593 6835069 : 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 2278344 : itest = itest + 1
600 2278344 : IF (gqsq < 1.0e2_dp*gisq) itest = 0
601 2278344 : IF (itest >= 3) THEN
602 381420 : DO i = 1, n
603 381420 : gq(i) = w(i)
604 : END DO
605 508560 : DO ih = 1, nh
606 508560 : hq(ih) = zero
607 : END DO
608 381420 : DO j = 1, nptm
609 254280 : w(j) = zero
610 1525680 : DO k = 1, npt
611 1525680 : w(j) = w(j) + vlag(k)*zmat(k, j)
612 : END DO
613 381420 : IF (j < idz) w(j) = -w(j)
614 : END DO
615 762840 : DO k = 1, npt
616 635700 : pq(k) = zero
617 2034240 : DO j = 1, nptm
618 1907100 : pq(k) = pq(k) + zmat(k, j)*w(j)
619 : END DO
620 : END DO
621 127140 : itest = 0
622 : END IF
623 : END IF
624 : END IF
625 53830933 : 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 53830933 : IF (f <= fsave + tenth*vquad .OR. ksave > 0) THEN
632 : skip_check = .TRUE.
633 : ELSE
634 12534742 : 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 64743625 : IF (.NOT. skip_check) THEN
643 22673584 : skip_check = .FALSE.
644 22673584 : distsq = 4.0_dp*delta*delta
645 136048060 : DO k = 1, npt
646 113374476 : sum = zero
647 340177770 : DO j = 1, n
648 340177770 : sum = sum + (xpt(k, j) - xopt(j))**2
649 : END DO
650 136048060 : IF (sum > distsq) THEN
651 30994976 : knew = k
652 30994976 : 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 22673584 : IF (knew > 0) THEN
660 19195563 : dstep = MAX(MIN(tenth*SQRT(distsq), half*delta), rho)
661 19195563 : dsq = dstep*dstep
662 19195563 : CALL generate_next_model_step
663 19195563 : RETURN
664 : END IF
665 3478021 : 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 3030223 : CALL update_rho()
671 3030223 : IF (do_return) RETURN
672 : END IF
673 : END IF
674 :
675 46387276 : skip_check = .FALSE.
676 : inner: DO
677 46387276 : knew = 0
678 46387276 : CALL trsapp(n, npt, xopt, xpt, gq, hq, pq, delta, d, w, w(np), w(np + n), w(np + 2*n), crvmin)
679 46387276 : dsq = zero
680 139167609 : DO i = 1, n
681 139167609 : dsq = dsq + d(i)**2
682 : END DO
683 46387276 : dnorm = MIN(delta, SQRT(dsq))
684 46387276 : IF (dnorm < half*rho) THEN
685 11293725 : knew = -1
686 11293725 : delta = tenth*delta
687 11293725 : ratio = -1.0_dp
688 11293725 : IF (delta <= 1.5_dp*rho) delta = rho
689 11293725 : IF (nf <= nfsav + 2) CYCLE outer
690 3127961 : temp = 0.125_dp*crvmin*rho*rho
691 3127961 : IF (temp <= MAX(diffa, diffb, diffc)) CYCLE outer
692 1613055 : CALL update_rho()
693 1613055 : IF (do_return) RETURN
694 : CYCLE inner
695 : END IF
696 : EXIT inner
697 : END DO inner
698 9680670 : 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 35093551 : CALL generate_next_model_step
704 :
705 : CONTAINS
706 : ! **************************************************************************************************
707 : !> \brief ...
708 : ! **************************************************************************************************
709 4643278 : SUBROUTINE update_rho()
710 4643278 : IF (rho > rhoend) THEN
711 3869437 : delta = half*rho
712 3869437 : ratio = rho/rhoend
713 3869437 : IF (ratio <= 16.0_dp) THEN
714 773835 : rho = rhoend
715 3095602 : ELSE IF (ratio <= 250.0_dp) THEN
716 773835 : rho = SQRT(ratio)*rhoend
717 : ELSE
718 2321767 : rho = tenth*rho
719 : END IF
720 3869437 : delta = MAX(delta, rho)
721 3869437 : 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 773841 : IF (knew == -1) THEN
728 645970 : CALL calc_next_value_of_objective_func
729 : ELSE
730 127871 : opt%state = 7
731 127871 : CALL get_state
732 :
733 127871 : CALL finalize_optimization_after_failure
734 : END IF
735 773841 : do_return = .TRUE.
736 : END IF
737 4643278 : END SUBROUTINE update_rho
738 : ! **************************************************************************************************
739 : !> \brief ...
740 : ! **************************************************************************************************
741 3869535 : SUBROUTINE begin_initialization_of_interpolation()
742 3869535 : nfm = nf
743 3869535 : nfmm = nf - n
744 3869535 : nf = nf + 1
745 3869535 : IF (nfm <= 2*n) THEN
746 3869535 : IF (nfm >= 1 .AND. nfm <= N) THEN
747 1547875 : xpt(nf, nfm) = rhobeg
748 2321660 : ELSE IF (nfm > n) THEN
749 1547801 : 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 11615465 : DO j = 1, n
773 11615465 : x(j) = xpt(nf, j) + xbase(j)
774 : END DO
775 3869535 : CALL check_number_of_steps
776 3869535 : END SUBROUTINE begin_initialization_of_interpolation
777 : ! **************************************************************************************************
778 : !> \brief ...
779 : ! **************************************************************************************************
780 54289114 : SUBROUTINE generate_next_model_step()
781 54289114 : IF (dsq <= 1.0e-3_dp*xoptsq) THEN
782 3326598 : tempq = 0.25_dp*xoptsq
783 19959908 : DO k = 1, npt
784 16633310 : sum = zero
785 49902138 : DO i = 1, n
786 49902138 : sum = sum + xpt(k, i)*xopt(i)
787 : END DO
788 16633310 : temp = pq(k)*sum
789 16633310 : sum = sum - half*xoptsq
790 16633310 : w(npt + k) = sum
791 53228736 : DO i = 1, n
792 33268828 : gq(i) = gq(i) + temp*xpt(k, i)
793 33268828 : xpt(k, i) = xpt(k, i) - half*xopt(i)
794 33268828 : vlag(i) = bmat(k, i)
795 33268828 : w(i) = sum*xpt(k, i) + tempq*xopt(i)
796 33268828 : ip = npt + i
797 99814908 : DO j = 1, i
798 83181598 : 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 9979954 : DO k = 1, nptm
806 6653356 : sumz = zero
807 39922184 : DO i = 1, npt
808 33268828 : sumz = sumz + zmat(i, k)
809 39922184 : w(i) = w(npt + i)*zmat(i, k)
810 : END DO
811 19961092 : DO j = 1, n
812 13307736 : sum = tempq*sumz*xopt(j)
813 79864448 : DO i = 1, npt
814 66556712 : sum = sum + w(i)*xpt(i, j)
815 66556712 : vlag(j) = sum
816 79864448 : IF (k < idz) sum = -sum
817 : END DO
818 86517804 : DO i = 1, npt
819 79864448 : bmat(i, j) = bmat(i, j) + sum*zmat(i, k)
820 : END DO
821 : END DO
822 23287690 : DO i = 1, n
823 13307736 : ip = i + npt
824 13307736 : temp = vlag(i)
825 13307736 : IF (k < idz) temp = -temp
826 39927204 : DO j = 1, i
827 33273848 : 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 3326598 : ih = 0
836 9979954 : DO j = 1, n
837 6653356 : w(j) = zero
838 39922184 : DO k = 1, npt
839 33268828 : w(j) = w(j) + pq(k)*xpt(k, j)
840 39922184 : xpt(k, j) = xpt(k, j) - half*xopt(j)
841 : END DO
842 19960500 : DO i = 1, j
843 9980546 : ih = ih + 1
844 9980546 : IF (i < j) gq(j) = gq(j) + hq(ih)*xopt(i)
845 9980546 : gq(i) = gq(i) + hq(ih)*xopt(j)
846 9980546 : hq(ih) = hq(ih) + w(i)*xopt(j) + xopt(i)*w(j)
847 16633902 : bmat(npt + i, j) = bmat(npt + j, i)
848 : END DO
849 : END DO
850 9979954 : DO j = 1, n
851 6653356 : xbase(j) = xbase(j) + xopt(j)
852 9979954 : xopt(j) = zero
853 : END DO
854 3326598 : 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 54289114 : IF (knew > 0) THEN
862 : CALL biglag(n, npt, xopt, xpt, bmat, zmat, idz, ndim, knew, dstep, &
863 19195563 : 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 325750496 : DO k = 1, npt
870 271461382 : suma = zero
871 271461382 : sumb = zero
872 271461382 : sum = zero
873 814527768 : DO j = 1, n
874 543066386 : suma = suma + xpt(k, j)*d(j)
875 543066386 : sumb = sumb + xpt(k, j)*xopt(j)
876 814527768 : sum = sum + bmat(k, j)*d(j)
877 : END DO
878 271461382 : w(k) = suma*(half*suma + sumb)
879 325750496 : vlag(k) = sum
880 : END DO
881 54289114 : beta = zero
882 162875248 : DO k = 1, nptm
883 108586134 : sum = zero
884 651652520 : DO i = 1, npt
885 651652520 : sum = sum + zmat(i, k)*w(i)
886 : END DO
887 108586134 : IF (k < idz) THEN
888 0 : beta = beta + sum*sum
889 0 : sum = -sum
890 : ELSE
891 108586134 : beta = beta - sum*sum
892 : END IF
893 705941634 : DO i = 1, npt
894 651652520 : vlag(i) = vlag(i) + sum*zmat(i, k)
895 : END DO
896 : END DO
897 54289114 : bsum = zero
898 54289114 : dx = zero
899 162875248 : DO j = 1, n
900 108586134 : sum = zero
901 651652520 : DO i = 1, npt
902 651652520 : sum = sum + w(i)*bmat(i, j)
903 : END DO
904 108586134 : bsum = bsum + sum*d(j)
905 108586134 : jp = npt + j
906 325826260 : DO k = 1, n
907 325826260 : sum = sum + bmat(jp, k)*d(k)
908 : END DO
909 108586134 : vlag(jp) = sum
910 108586134 : bsum = bsum + sum*d(j)
911 162875248 : dx = dx + d(j)*xopt(j)
912 : END DO
913 54289114 : beta = dx*dx + dsq*(xoptsq + dx + dx + half*dsq) + beta - bsum
914 54289114 : 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 54289114 : IF (knew > 0) THEN
921 19195563 : temp = one + alpha*beta/vlag(knew)**2
922 19195563 : 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 54289114 : CALL calc_next_value_of_objective_func
931 54289114 : END SUBROUTINE generate_next_model_step
932 : ! **************************************************************************************************
933 : !> \brief ...
934 : ! **************************************************************************************************
935 54935084 : SUBROUTINE calc_next_value_of_objective_func()
936 164813196 : DO i = 1, n
937 109878112 : xnew(i) = xopt(i) + d(i)
938 164813196 : x(i) = xbase(i) + xnew(i)
939 : END DO
940 54935084 : nf = nf + 1
941 54935084 : CALL check_number_of_steps
942 54935084 : END SUBROUTINE calc_next_value_of_objective_func
943 : ! **************************************************************************************************
944 : !> \brief ...
945 : ! **************************************************************************************************
946 58804619 : SUBROUTINE check_number_of_steps()
947 58804619 : IF (nf > nftest) THEN
948 : ! return to many steps
949 13 : nf = nf - 1
950 13 : opt%state = 3
951 13 : CALL get_state
952 13 : CALL finalize_optimization_after_failure
953 : ELSE
954 :
955 58804606 : CALL get_state
956 :
957 58804606 : opt%state = 2
958 : END IF
959 58804619 : END SUBROUTINE check_number_of_steps
960 : ! **************************************************************************************************
961 : !> \brief ...
962 : ! **************************************************************************************************
963 773858 : SUBROUTINE finalize_optimization_after_failure()
964 773858 : IF (fopt <= f) THEN
965 722384 : DO i = 1, n
966 722384 : x(i) = xbase(i) + xopt(i)
967 : END DO
968 240677 : f = fopt
969 : END IF
970 :
971 773858 : CALL get_state
972 773858 : END SUBROUTINE finalize_optimization_after_failure
973 : ! **************************************************************************************************
974 : !> \brief ...
975 : ! **************************************************************************************************
976 60352322 : SUBROUTINE get_state()
977 60352322 : opt%np = np
978 60352322 : opt%nh = nh
979 60352322 : opt%nptm = nptm
980 60352322 : opt%nftest = nftest
981 60352322 : opt%idz = idz
982 60352322 : opt%itest = itest
983 60352322 : opt%nf = nf
984 60352322 : opt%nfm = nfm
985 60352322 : opt%nfmm = nfmm
986 60352322 : opt%nfsav = nfsav
987 60352322 : opt%knew = knew
988 60352322 : opt%kopt = kopt
989 60352322 : opt%ksave = ksave
990 60352322 : opt%ktemp = ktemp
991 60352322 : opt%rhosq = rhosq
992 60352322 : opt%recip = recip
993 60352322 : opt%reciq = reciq
994 60352322 : opt%fbeg = fbeg
995 60352322 : opt%fopt = fopt
996 60352322 : opt%diffa = diffa
997 60352322 : opt%xoptsq = xoptsq
998 60352322 : opt%rho = rho
999 60352322 : opt%delta = delta
1000 60352322 : opt%dsq = dsq
1001 60352322 : opt%dnorm = dnorm
1002 60352322 : opt%ratio = ratio
1003 60352322 : opt%temp = temp
1004 60352322 : opt%tempq = tempq
1005 60352322 : opt%beta = beta
1006 60352322 : opt%dx = dx
1007 60352322 : opt%vquad = vquad
1008 60352322 : opt%diff = diff
1009 60352322 : opt%diffc = diffc
1010 60352322 : opt%diffb = diffb
1011 60352322 : opt%fsave = fsave
1012 60352322 : opt%detrat = detrat
1013 60352322 : opt%hdiag = hdiag
1014 60352322 : opt%distsq = distsq
1015 60352322 : opt%gisq = gisq
1016 60352322 : opt%gqsq = gqsq
1017 60352322 : opt%f = f
1018 60352322 : opt%bstep = bstep
1019 60352322 : opt%alpha = alpha
1020 60352322 : opt%dstep = dstep
1021 60352322 : END SUBROUTINE get_state
1022 : ! **************************************************************************************************
1023 : !> \brief ...
1024 : ! **************************************************************************************************
1025 58804605 : SUBROUTINE set_state()
1026 58804605 : np = opt%np
1027 58804605 : nh = opt%nh
1028 58804605 : nptm = opt%nptm
1029 58804605 : nftest = opt%nftest
1030 58804605 : idz = opt%idz
1031 58804605 : itest = opt%itest
1032 58804605 : nf = opt%nf
1033 58804605 : nfm = opt%nfm
1034 58804605 : nfmm = opt%nfmm
1035 58804605 : nfsav = opt%nfsav
1036 58804605 : knew = opt%knew
1037 58804605 : kopt = opt%kopt
1038 58804605 : ksave = opt%ksave
1039 58804605 : ktemp = opt%ktemp
1040 58804605 : rhosq = opt%rhosq
1041 58804605 : recip = opt%recip
1042 58804605 : reciq = opt%reciq
1043 58804605 : fbeg = opt%fbeg
1044 58804605 : fopt = opt%fopt
1045 58804605 : diffa = opt%diffa
1046 58804605 : xoptsq = opt%xoptsq
1047 58804605 : rho = opt%rho
1048 58804605 : delta = opt%delta
1049 58804605 : dsq = opt%dsq
1050 58804605 : dnorm = opt%dnorm
1051 58804605 : ratio = opt%ratio
1052 58804605 : temp = opt%temp
1053 58804605 : tempq = opt%tempq
1054 58804605 : beta = opt%beta
1055 58804605 : dx = opt%dx
1056 58804605 : vquad = opt%vquad
1057 58804605 : diff = opt%diff
1058 58804605 : diffc = opt%diffc
1059 58804605 : diffb = opt%diffb
1060 58804605 : fsave = opt%fsave
1061 58804605 : detrat = opt%detrat
1062 58804605 : hdiag = opt%hdiag
1063 58804605 : distsq = opt%distsq
1064 58804605 : gisq = opt%gisq
1065 58804605 : gqsq = opt%gqsq
1066 58804605 : f = opt%f
1067 58804605 : bstep = opt%bstep
1068 58804605 : alpha = opt%alpha
1069 58804605 : dstep = opt%dstep
1070 58804605 : 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 19195563 : 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 19195563 : delsq = delta*delta
1527 19195563 : 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 19195563 : iterc = 0
1533 115179286 : DO k = 1, npt
1534 115179286 : hcol(k) = zero
1535 : END DO
1536 57589643 : DO j = 1, nptm
1537 38394080 : temp = zmat(knew, j)
1538 38394080 : IF (j < idz) temp = -temp
1539 249606987 : DO k = 1, npt
1540 230411424 : hcol(k) = hcol(k) + temp*zmat(k, j)
1541 : END DO
1542 : END DO
1543 19195563 : 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 19195563 : dd = zero
1549 57589643 : DO i = 1, n
1550 38394080 : d(i) = xpt(knew, i) - xopt(i)
1551 38394080 : gc(i) = bmat(knew, i)
1552 38394080 : gd(i) = zero
1553 57589643 : dd = dd + d(i)**2
1554 : END DO
1555 115179286 : DO k = 1, npt
1556 : temp = zero
1557 : sum = zero
1558 288001067 : DO j = 1, n
1559 192017344 : temp = temp + xpt(k, j)*xopt(j)
1560 288001067 : sum = sum + xpt(k, j)*d(j)
1561 : END DO
1562 95983723 : temp = hcol(k)*temp
1563 95983723 : sum = hcol(k)*sum
1564 307196630 : DO i = 1, n
1565 192017344 : gc(i) = gc(i) + temp*xpt(k, i)
1566 288001067 : 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 57589643 : DO i = 1, n
1577 38394080 : gg = gg + gc(i)**2
1578 38394080 : sp = sp + d(i)*gc(i)
1579 57589643 : dhd = dhd + d(i)*gd(i)
1580 : END DO
1581 19195563 : scale = delta/SQRT(dd)
1582 19195563 : IF (sp*dhd < zero) scale = -scale
1583 19195563 : temp = zero
1584 19195563 : IF (sp*sp > 0.99_dp*dd*gg) temp = one
1585 19195563 : tau = scale*(ABS(sp) + half*scale*ABS(dhd))
1586 19195563 : IF (gg*delsq < 0.01_dp*tau*tau) temp = one
1587 57589643 : DO i = 1, n
1588 38394080 : d(i) = scale*d(i)
1589 38394080 : gd(i) = scale*gd(i)
1590 57589643 : 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 27179162 : iterc = iterc + 1
1599 27179162 : dd = zero
1600 27179162 : sp = zero
1601 27179162 : ss = zero
1602 81542496 : DO i = 1, n
1603 54363334 : dd = dd + d(i)**2
1604 54363334 : sp = sp + d(i)*s(i)
1605 81542496 : ss = ss + s(i)**2
1606 : END DO
1607 27179162 : temp = dd*ss - sp*sp
1608 27179162 : IF (temp <= 1.0e-8_dp*dd*ss) EXIT mainloop
1609 25996537 : denom = SQRT(temp)
1610 77994347 : DO i = 1, n
1611 51997810 : s(i) = (dd*s(i) - sp*d(i))/denom
1612 77994347 : 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 155988694 : DO k = 1, npt
1619 : sum = zero
1620 390065395 : DO j = 1, n
1621 390065395 : sum = sum + xpt(k, j)*s(j)
1622 : END DO
1623 129992157 : sum = hcol(k)*sum
1624 416061932 : DO i = 1, n
1625 390065395 : 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 77994347 : DO i = 1, n
1634 51997810 : cf1 = cf1 + s(i)*w(i)
1635 51997810 : cf2 = cf2 + d(i)*gc(i)
1636 51997810 : cf3 = cf3 + s(i)*gc(i)
1637 51997810 : cf4 = cf4 + d(i)*gd(i)
1638 77994347 : cf5 = cf5 + s(i)*gd(i)
1639 : END DO
1640 25996537 : cf1 = half*cf1
1641 25996537 : cf4 = half*cf4 - cf1
1642 : !
1643 : ! Seek the value of the angle that maximizes the modulus of TAU.
1644 : !
1645 25996537 : taubeg = cf1 + cf2 + cf4
1646 25996537 : taumax = taubeg
1647 25996537 : tauold = taubeg
1648 25996537 : isave = 0
1649 25996537 : iu = 49
1650 25996537 : temp = twopi/REAL(iu + 1, DP)
1651 1299826850 : DO i = 1, iu
1652 1273830313 : angle = REAL(i, dp)*temp
1653 1273830313 : cth = COS(angle)
1654 1273830313 : sth = SIN(angle)
1655 1273830313 : tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1656 1273830313 : IF (ABS(tau) > ABS(taumax)) THEN
1657 : taumax = tau
1658 : isave = i
1659 : tempa = tauold
1660 1207285054 : ELSE IF (i == isave + 1) THEN
1661 28265083 : tempb = taU
1662 : END IF
1663 1299826850 : tauold = tau
1664 : END DO
1665 25996537 : IF (isave == 0) tempa = tau
1666 15269864 : IF (isave == iu) tempb = taubeg
1667 25996537 : step = zero
1668 25996537 : IF (tempa /= tempb) THEN
1669 25996537 : tempa = tempa - taumax
1670 25996537 : tempb = tempb - taumax
1671 25996537 : step = half*(tempa - tempb)/(tempa + tempb)
1672 : END IF
1673 25996537 : angle = temp*(REAL(isave, DP) + step)
1674 : !
1675 : ! Calculate the new D and GD. Then test for convergence.
1676 : !
1677 25996537 : cth = COS(angle)
1678 25996537 : sth = SIN(angle)
1679 25996537 : tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1680 77994347 : DO i = 1, n
1681 51997810 : d(i) = cth*d(i) + sth*s(i)
1682 51997810 : gd(i) = cth*gd(i) + sth*w(i)
1683 77994347 : s(i) = gc(i) + gd(i)
1684 : END DO
1685 25996537 : IF (ABS(tau) <= 1.1_dp*ABS(taubeg)) EXIT mainloop
1686 27179162 : IF (iterc >= n) EXIT mainloop
1687 : END DO mainloop
1688 :
1689 19195563 : 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 46387276 : 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 46387276 : delsq = delta*delta
1749 46387276 : iterc = 0
1750 46387276 : itermax = n
1751 46387276 : itersw = itermax
1752 139167609 : DO i = 1, n
1753 139167609 : d(i) = xopt(i)
1754 : END DO
1755 46387276 : CALL updatehd
1756 : !
1757 : ! Prepare for the first line search.
1758 : !
1759 46387276 : qred = zero
1760 46387276 : dd = zero
1761 139167609 : DO i = 1, n
1762 92780333 : step(i) = zero
1763 92780333 : hs(i) = zero
1764 92780333 : g(i) = gq(i) + hd(i)
1765 92780333 : d(i) = -g(i)
1766 139167609 : dd = dd + d(i)**2
1767 : END DO
1768 46387276 : crvmin = zero
1769 46387276 : 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 78303584 : IF (.NOT. jump2) THEN
1781 78302555 : IF (.NOT. jump1) THEN
1782 78302555 : iterc = iterc + 1
1783 78302555 : temp = delsq - ss
1784 78302555 : bstep = temp/(ds + SQRT(ds*ds + dd*temp))
1785 78302555 : CALL updatehd
1786 : END IF
1787 78302555 : jump1 = .FALSE.
1788 78302555 : IF (iterc <= itersw) THEN
1789 78302555 : dhd = zero
1790 234923085 : DO j = 1, n
1791 234923085 : dhd = dhd + d(j)*hd(j)
1792 : END DO
1793 : !
1794 : ! Update CRVMIN and set the step-length ALPHA.
1795 : !
1796 78302555 : alpha = bstep
1797 78302555 : IF (dhd > zero) THEN
1798 69200298 : temp = dhd/dd
1799 69200298 : IF (iterc == 1) crvmin = temp
1800 69200298 : crvmin = MIN(crvmin, temp)
1801 69200298 : alpha = MIN(alpha, gg/dhd)
1802 : END IF
1803 78302555 : qadd = alpha*(gg - half*alpha*dhd)
1804 78302555 : qred = qred + qadd
1805 : !
1806 : ! Update STEP and HS.
1807 : !
1808 78302555 : ggsav = gg
1809 78302555 : gg = zero
1810 234923085 : DO i = 1, n
1811 156620530 : step(i) = step(i) + alpha*d(i)
1812 156620530 : hs(i) = hs(i) + alpha*hd(i)
1813 234923085 : gg = gg + (g(i) + hs(i))**2
1814 : END DO
1815 : !
1816 : ! Begin another conjugate direction iteration if required.
1817 : !
1818 78302555 : IF (alpha < bstep) THEN
1819 51469743 : IF (qadd <= 0.01_dp*qred) EXIT mainloop
1820 50689354 : IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1821 31915484 : IF (iterc == itermax) EXIT mainloop
1822 31915432 : temp = gg/ggsav
1823 31915432 : dd = zero
1824 31915432 : ds = zero
1825 31915432 : ss = zero
1826 95756266 : DO i = 1, n
1827 63840834 : d(i) = temp*d(i) - g(i) - hs(i)
1828 63840834 : dd = dd + d(i)**2
1829 63840834 : ds = ds + d(i)*step(I)
1830 95756266 : ss = ss + step(i)**2
1831 : END DO
1832 31915432 : IF (ds <= zero) EXIT mainloop
1833 31915432 : IF (ss < delsq) CYCLE mainloop
1834 : END IF
1835 26832812 : crvmin = zero
1836 26832812 : itersw = iterc
1837 26832812 : jump2 = .TRUE.
1838 26832812 : 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 26432809 : sg = zero
1849 26432809 : shs = zero
1850 79309512 : DO i = 1, n
1851 52876703 : sg = sg + step(i)*g(i)
1852 79309512 : shs = shs + step(i)*hs(i)
1853 : END DO
1854 26432809 : sgk = sg + shs
1855 26432809 : angtest = sgk/SQRT(gg*delsq)
1856 26432809 : 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 23078979 : iterc = iterc + 1
1862 23078979 : temp = SQRT(delsq*gg - sgk*sgk)
1863 23078979 : tempa = delsq/temp
1864 23078979 : tempb = sgk/temp
1865 69247282 : DO i = 1, n
1866 69247282 : d(i) = tempa*(g(i) + hs(i)) - tempb*step(i)
1867 : END DO
1868 23078979 : CALL updatehd
1869 23078979 : IF (iterc <= itersw) THEN
1870 : jump1 = .TRUE.
1871 : CYCLE mainloop
1872 : END IF
1873 : END IF
1874 23078979 : dg = zero
1875 23078979 : dhd = zero
1876 23078979 : dhs = zero
1877 69247282 : DO i = 1, n
1878 46168303 : dg = dg + d(i)*g(i)
1879 46168303 : dhd = dhd + hd(i)*d(i)
1880 69247282 : dhs = dhs + hd(i)*step(i)
1881 : END DO
1882 : !
1883 : ! Seek the value of the angle that minimizes Q.
1884 : !
1885 23078979 : cf = half*(shs - dhd)
1886 23078979 : qbeg = sg + cf
1887 23078979 : qsav = qbeg
1888 23078979 : qmin = qbeg
1889 23078979 : isave = 0
1890 23078979 : iu = 49
1891 : temp = twopi/REAL(iu + 1, dp)
1892 1153948950 : DO i = 1, iu
1893 1130869971 : angle = REAL(i, dp)*temp
1894 1130869971 : cth = COS(angle)
1895 1130869971 : sth = SIN(angle)
1896 1130869971 : qnew = (sg + cf*cth)*cth + (dg + dhs*cth)*sth
1897 1130869971 : IF (qnew < qmin) THEN
1898 : qmin = qnew
1899 : isave = i
1900 : tempa = qsav
1901 1106848044 : ELSE IF (i == isave + 1) THEN
1902 28789623 : tempb = qnew
1903 : END IF
1904 1153948950 : qsav = qnew
1905 : END DO
1906 23078979 : IF (isave == zero) tempa = qnew
1907 9996550 : IF (isave == iu) tempb = qbeg
1908 23078979 : angle = zero
1909 23078979 : IF (tempa /= tempb) THEN
1910 23078979 : tempa = tempa - qmin
1911 23078979 : tempb = tempb - qmin
1912 23078979 : angle = half*(tempa - tempb)/(tempa + tempb)
1913 : END IF
1914 23078979 : angle = temp*(REAL(isave, DP) + angle)
1915 : !
1916 : ! Calculate the new STEP and HS. Then test for convergence.
1917 : !
1918 23078979 : cth = COS(angle)
1919 23078979 : sth = SIN(angle)
1920 23078979 : reduc = qbeg - (sg + cf*cth)*cth - (dg + dhs*cth)*sth
1921 23078979 : gg = zero
1922 69247282 : DO i = 1, n
1923 46168303 : step(i) = cth*step(i) + sth*d(i)
1924 46168303 : hs(i) = cth*hs(i) + sth*hd(i)
1925 69247282 : gg = gg + (g(i) + hs(i))**2
1926 : END DO
1927 23078979 : qred = qred + reduc
1928 23078979 : ratio = reduc/qred
1929 23078979 : IF (iterc < itermax .AND. ratio > 0.01_dp) THEN
1930 1029 : jump2 = .TRUE.
1931 : ELSE
1932 : EXIT mainloop
1933 : END IF
1934 :
1935 46388152 : IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1936 :
1937 : END DO mainloop
1938 :
1939 : CONTAINS
1940 : ! **************************************************************************************************
1941 : !> \brief ...
1942 : ! **************************************************************************************************
1943 147768810 : SUBROUTINE updatehd
1944 : INTEGER :: i, ih, j, k
1945 :
1946 443337976 : DO i = 1, n
1947 443337976 : hd(i) = zero
1948 : END DO
1949 886675952 : DO k = 1, npt
1950 738907142 : temp = zero
1951 2217319512 : DO j = 1, n
1952 2217319512 : temp = temp + xpt(k, j)*d(j)
1953 : END DO
1954 738907142 : temp = temp*pq(k)
1955 2365088322 : DO i = 1, n
1956 2217319512 : hd(i) = hd(i) + temp*xpt(k, i)
1957 : END DO
1958 : END DO
1959 147768810 : ih = 0
1960 443337976 : DO j = 1, n
1961 886833360 : DO i = 1, j
1962 443495384 : ih = ih + 1
1963 443495384 : IF (i < j) hd(j) = hd(j) + hq(ih)*d(i)
1964 739064550 : hd(i) = hd(i) + hq(ih)*d(j)
1965 : END DO
1966 : END DO
1967 147768810 : 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 53830933 : 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 53830933 : nptm = npt - n - 1
2009 : !
2010 : ! Apply the rotations that put zeros in the KNEW-th row of ZMAT.
2011 : !
2012 53830933 : jl = 1
2013 107669681 : DO j = 2, nptm
2014 107669681 : IF (j == idz) THEN
2015 : jl = idz
2016 53838748 : ELSE IF (zmat(knew, j) /= zero) THEN
2017 52377013 : temp = SQRT(zmat(knew, jl)**2 + zmat(knew, j)**2)
2018 52377013 : tempa = zmat(knew, jl)/temp
2019 52377013 : tempb = zmat(knew, j)/temp
2020 314369390 : DO I = 1, NPT
2021 261992377 : temp = tempa*zmat(i, jl) + tempb*zmat(i, j)
2022 261992377 : zmat(i, j) = tempa*zmat(i, j) - tempb*zmat(i, jl)
2023 314369390 : zmat(i, jl) = temp
2024 : END DO
2025 52377013 : 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 53830933 : tempa = zmat(knew, 1)
2033 53830933 : IF (idz >= 2) tempa = -tempa
2034 53830933 : IF (jl > 1) tempb = zmat(knew, jl)
2035 323001228 : DO i = 1, npt
2036 269170295 : w(i) = tempa*zmat(i, 1)
2037 323001228 : IF (jl > 1) w(i) = w(i) + tempb*zmat(i, jl)
2038 : END DO
2039 53830933 : alpha = w(knew)
2040 53830933 : tau = vlag(knew)
2041 53830933 : tausq = tau*tau
2042 53830933 : denom = alpha*beta + tausq
2043 53830933 : 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 53830933 : iflag = 0
2050 53830933 : IF (JL == 1) THEN
2051 53830933 : temp = SQRT(ABS(denom))
2052 53830933 : tempb = tempa/temp
2053 53830933 : tempa = tau/temp
2054 323001228 : DO i = 1, npt
2055 323001228 : zmat(i, 1) = tempa*zmat(i, 1) - tempb*vlag(i)
2056 : END DO
2057 53830933 : IF (idz == 1 .AND. temp < zero) idz = 2
2058 53830933 : 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 161500614 : DO j = 1, n
2097 107669681 : jp = npt + j
2098 107669681 : w(jp) = bmat(knew, j)
2099 107669681 : tempa = (alpha*vlag(jp) - tau*w(jp))/denom
2100 107669681 : tempb = (-beta*w(jp) - tau*vlag(jp))/denom
2101 861520923 : DO i = 1, jp
2102 700020309 : bmat(i, j) = bmat(i, j) + tempa*vlag(i) + tempb*w(i)
2103 807689990 : IF (i > npt) bmat(jp, i - npt) = bmat(i, j)
2104 : END DO
2105 : END DO
2106 :
2107 53830933 : END SUBROUTINE update
2108 :
2109 0 : END MODULE powell
|
