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