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