(git:374b731)
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-2024 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
43!******************************************************************************
44
45CONTAINS
46
47!******************************************************************************
48! **************************************************************************************************
49!> \brief ...
50!> \param n ...
51!> \param x ...
52!> \param optstate ...
53! **************************************************************************************************
54 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 CALL timeset(routinen, handle)
64
65 SELECT CASE (optstate%state)
66 CASE (0)
67 npt = 2*n + 1
68 ALLOCATE (optstate%w((npt + 13)*(npt + n) + 3*n*(n + 3)/2))
69 ALLOCATE (optstate%xopt(n))
70 ! Initialize w
71 optstate%w = 0.0_dp
72 optstate%state = 1
73 CALL newuoa(n, x, optstate)
74 CASE (1, 2)
75 CALL newuoa(n, x, optstate)
76 CASE (3)
77 IF (optstate%unit > 0) THEN
78 WRITE (optstate%unit, *) "POWELL| Exceeding maximum number of steps"
79 END IF
80 optstate%state = -1
81 CASE (4)
82 IF (optstate%unit > 0) THEN
83 WRITE (optstate%unit, *) "POWELL| Error in trust region"
84 END IF
85 optstate%state = -1
86 CASE (5)
87 IF (optstate%unit > 0) THEN
88 WRITE (optstate%unit, *) "POWELL| N out of range"
89 END IF
90 optstate%state = -1
91 CASE (6, 7)
92 optstate%state = -1
93 CASE (8)
94 x(1:n) = optstate%xopt(1:n)
95 DEALLOCATE (optstate%w)
96 DEALLOCATE (optstate%xopt)
97 optstate%state = -1
98 CASE DEFAULT
99 cpabort("")
100 END SELECT
101
102 CALL timestop(handle)
103
104 END SUBROUTINE powell_optimize
105!******************************************************************************
106! **************************************************************************************************
107!> \brief ...
108!> \param n ...
109!> \param x ...
110!> \param optstate ...
111! **************************************************************************************************
112 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 maxfun = optstate%maxfun
124 rhobeg = optstate%rhobeg
125 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 np = n + 1
162 npt = 2*n + 1
163 nptm = npt - np
164 IF (npt < n + 2 .OR. npt > ((n + 2)*np)/2) THEN
165 optstate%state = 5
166 RETURN
167 END IF
168 ndim = npt + n
169 ixb = 1
170 ixo = ixb + n
171 ixn = ixo + n
172 ixp = ixn + n
173 ifv = ixp + n*npt
174 igq = ifv + npt
175 ihq = igq + n
176 ipq = ihq + (n*np)/2
177 ibmat = ipq + npt
178 izmat = ibmat + ndim*n
179 id = izmat + npt*nptm
180 ivl = id + n
181 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 optstate%w(ivl:), optstate%w(iw:), optstate)
191
192 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 SUBROUTINE newuob(n, npt, x, rhobeg, rhoend, maxfun, xbase, &
222 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 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 idz = 0
283 itest = 0
284 nf = 0
285 nfm = 0
286 nfmm = 0
287 nfsav = 0
288 knew = 0
289 kopt = 0
290 ksave = 0
291 ktemp = 0
292 rhosq = 0._dp
293 recip = 0._dp
294 reciq = 0._dp
295 fbeg = 0._dp
296 fopt = 0._dp
297 diffa = 0._dp
298 xoptsq = 0._dp
299 rho = 0._dp
300 delta = 0._dp
301 dsq = 0._dp
302 dnorm = 0._dp
303 ratio = 0._dp
304 temp = 0._dp
305 tempq = 0._dp
306 beta = 0._dp
307 dx = 0._dp
308 vquad = 0._dp
309 diff = 0._dp
310 diffc = 0._dp
311 diffb = 0._dp
312 fsave = 0._dp
313 detrat = 0._dp
314 hdiag = 0._dp
315 distsq = 0._dp
316 gisq = 0._dp
317 gqsq = 0._dp
318 f = 0._dp
319 bstep = 0._dp
320 alpha = 0._dp
321 dstep = 0._dp
322 !
323 END IF
324
325 ipt = 0
326 jpt = 0
327 xipt = 0._dp
328 xjpt = 0._dp
329
330 half = 0.5_dp
331 one = 1.0_dp
332 tenth = 0.1_dp
333 zero = 0.0_dp
334 np = n + 1
335 nh = (n*np)/2
336 nptm = npt - np
337 nftest = max(maxfun, 1)
338
339 IF (opt%state == 2) GOTO 1000
340 !
341 ! Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero.
342 !
343 DO j = 1, n
344 xbase(j) = x(j)
345 DO k = 1, npt
346 xpt(k, j) = zero
347 END DO
348 DO i = 1, ndim
349 bmat(i, j) = zero
350 END DO
351 END DO
352 DO ih = 1, nh
353 hq(ih) = zero
354 END DO
355 DO k = 1, npt
356 pq(k) = zero
357 DO j = 1, nptm
358 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 rhosq = rhobeg*rhobeg
367 recip = one/rhosq
368 reciq = sqrt(half)/rhosq
369 nf = 0
37050 nfm = nf
371 nfmm = nf - n
372 nf = nf + 1
373 IF (nfm <= 2*n) THEN
374 IF (nfm >= 1 .AND. nfm <= n) THEN
375 xpt(nf, nfm) = rhobeg
376 ELSE IF (nfm > n) THEN
377 xpt(nf, nfmm) = -rhobeg
378 END IF
379 ELSE
380 itemp = (nfmm - 1)/n
381 jpt = nfm - itemp*n - n
382 ipt = jpt + itemp
383 IF (ipt > n) THEN
384 itemp = jpt
385 jpt = ipt - n
386 ipt = itemp
387 END IF
388 xipt = rhobeg
389 IF (fval(ipt + np) < fval(ipt + 1)) xipt = -xipt
390 xjpt = rhobeg
391 IF (fval(jpt + np) < fval(jpt + 1)) xjpt = -xjpt
392 xpt(nf, ipt) = xipt
393 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 DO j = 1, n
401 x(j) = xpt(nf, j) + xbase(j)
402 END DO
403 GOTO 310
40470 fval(nf) = f
405 IF (nf == 1) THEN
406 fbeg = f
407 fopt = f
408 kopt = 1
409 ELSE IF (f < fopt) THEN
410 fopt = f
411 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 IF (nfm <= 2*n) THEN
418 IF (nfm >= 1 .AND. nfm <= n) THEN
419 gq(nfm) = (f - fbeg)/rhobeg
420 IF (npt < nf + n) THEN
421 bmat(1, nfm) = -one/rhobeg
422 bmat(nf, nfm) = one/rhobeg
423 bmat(npt + nfm, nfm) = -half*rhosq
424 END IF
425 ELSE IF (nfm > n) THEN
426 bmat(nf - n, nfmm) = half/rhobeg
427 bmat(nf, nfmm) = -half/rhobeg
428 zmat(1, nfmm) = -reciq - reciq
429 zmat(nf - n, nfmm) = reciq
430 zmat(nf, nfmm) = reciq
431 ih = (nfmm*(nfmm + 1))/2
432 temp = (fbeg - f)/rhobeg
433 hq(ih) = (gq(nfmm) - temp)/rhobeg
434 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 ih = (ipt*(ipt - 1))/2 + jpt
442 IF (xipt < zero) ipt = ipt + n
443 IF (xjpt < zero) jpt = jpt + n
444 zmat(1, nfmm) = recip
445 zmat(nf, nfmm) = recip
446 zmat(ipt + 1, nfmm) = -recip
447 zmat(jpt + 1, nfmm) = -recip
448 hq(ih) = (fbeg - fval(ipt + 1) - fval(jpt + 1) + f)/(xipt*xjpt)
449 END IF
450 IF (nf < npt) GOTO 50
451 !
452 ! Begin the iterative procedure, because the initial model is comple
453 !
454 rho = rhobeg
455 delta = rho
456 idz = 1
457 diffa = zero
458 diffb = zero
459 itest = 0
460 xoptsq = zero
461 DO i = 1, n
462 xopt(i) = xpt(kopt, i)
463 xoptsq = xoptsq + xopt(i)**2
464 END DO
46590 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 !
470100 knew = 0
471 CALL trsapp(n, npt, xopt, xpt, gq, hq, pq, delta, d, w, w(np), w(np + n), w(np + 2*n), crvmin)
472 dsq = zero
473 DO i = 1, n
474 dsq = dsq + d(i)**2
475 END DO
476 dnorm = min(delta, sqrt(dsq))
477 IF (dnorm < half*rho) THEN
478 knew = -1
479 delta = tenth*delta
480 ratio = -1.0_dp
481 IF (delta <= 1.5_dp*rho) delta = rho
482 IF (nf <= nfsav + 2) GOTO 460
483 temp = 0.125_dp*crvmin*rho*rho
484 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 !
491120 IF (dsq <= 1.0e-3_dp*xoptsq) THEN
492 tempq = 0.25_dp*xoptsq
493 DO k = 1, npt
494 sum = zero
495 DO i = 1, n
496 sum = sum + xpt(k, i)*xopt(i)
497 END DO
498 temp = pq(k)*sum
499 sum = sum - half*xoptsq
500 w(npt + k) = sum
501 DO i = 1, n
502 gq(i) = gq(i) + temp*xpt(k, i)
503 xpt(k, i) = xpt(k, i) - half*xopt(i)
504 vlag(i) = bmat(k, i)
505 w(i) = sum*xpt(k, i) + tempq*xopt(i)
506 ip = npt + i
507 DO j = 1, i
508 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 DO k = 1, nptm
516 sumz = zero
517 DO i = 1, npt
518 sumz = sumz + zmat(i, k)
519 w(i) = w(npt + i)*zmat(i, k)
520 END DO
521 DO j = 1, n
522 sum = tempq*sumz*xopt(j)
523 DO i = 1, npt
524 sum = sum + w(i)*xpt(i, j)
525 vlag(j) = sum
526 IF (k < idz) sum = -sum
527 END DO
528 DO i = 1, npt
529 bmat(i, j) = bmat(i, j) + sum*zmat(i, k)
530 END DO
531 END DO
532 DO i = 1, n
533 ip = i + npt
534 temp = vlag(i)
535 IF (k < idz) temp = -temp
536 DO j = 1, i
537 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 DO j = 1, n
547 w(j) = zero
548 DO k = 1, npt
549 w(j) = w(j) + pq(k)*xpt(k, j)
550 xpt(k, j) = xpt(k, j) - half*xopt(j)
551 END DO
552 DO i = 1, j
553 ih = ih + 1
554 IF (i < j) gq(j) = gq(j) + hq(ih)*xopt(i)
555 gq(i) = gq(i) + hq(ih)*xopt(j)
556 hq(ih) = hq(ih) + w(i)*xopt(j) + xopt(i)*w(j)
557 bmat(npt + i, j) = bmat(npt + j, i)
558 END DO
559 END DO
560 DO j = 1, n
561 xbase(j) = xbase(j) + xopt(j)
562 xopt(j) = zero
563 END DO
564 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 IF (knew > 0) THEN
572 CALL biglag(n, npt, xopt, xpt, bmat, zmat, idz, ndim, knew, dstep, &
573 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 DO k = 1, npt
580 suma = zero
581 sumb = zero
582 sum = zero
583 DO j = 1, n
584 suma = suma + xpt(k, j)*d(j)
585 sumb = sumb + xpt(k, j)*xopt(j)
586 sum = sum + bmat(k, j)*d(j)
587 END DO
588 w(k) = suma*(half*suma + sumb)
589 vlag(k) = sum
590 END DO
591 beta = zero
592 DO k = 1, nptm
593 sum = zero
594 DO i = 1, npt
595 sum = sum + zmat(i, k)*w(i)
596 END DO
597 IF (k < idz) THEN
598 beta = beta + sum*sum
599 sum = -sum
600 ELSE
601 beta = beta - sum*sum
602 END IF
603 DO i = 1, npt
604 vlag(i) = vlag(i) + sum*zmat(i, k)
605 END DO
606 END DO
607 bsum = zero
608 dx = zero
609 DO j = 1, n
610 sum = zero
611 DO i = 1, npt
612 sum = sum + w(i)*bmat(i, j)
613 END DO
614 bsum = bsum + sum*d(j)
615 jp = npt + j
616 DO k = 1, n
617 sum = sum + bmat(jp, k)*d(k)
618 END DO
619 vlag(jp) = sum
620 bsum = bsum + sum*d(j)
621 dx = dx + d(j)*xopt(j)
622 END DO
623 beta = dx*dx + dsq*(xoptsq + dx + dx + half*dsq) + beta - bsum
624 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 IF (knew > 0) THEN
631 temp = one + alpha*beta/vlag(knew)**2
632 IF (abs(temp) <= 0.8_dp) THEN
633 CALL bigden(n, npt, xopt, xpt, bmat, zmat, idz, ndim, kopt, &
634 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 !
640290 DO i = 1, n
641 xnew(i) = xopt(i) + d(i)
642 x(i) = xbase(i) + xnew(i)
643 END DO
644 nf = nf + 1
645310 IF (nf > nftest) THEN
646 ! return to many steps
647 nf = nf - 1
648 opt%state = 3
649 CALL get_state
650 GOTO 530
651 END IF
652
653 CALL get_state
654
655 opt%state = 2
656
657 RETURN
658
6591000 CONTINUE
660
661 CALL set_state
662
663 IF (nf <= npt) GOTO 70
664 IF (knew == -1) THEN
665 opt%state = 6
666 CALL get_state
667 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 vquad = zero
674 ih = 0
675 DO j = 1, n
676 vquad = vquad + d(j)*gq(j)
677 DO i = 1, j
678 ih = ih + 1
679 temp = d(i)*xnew(j) + d(j)*xopt(i)
680 IF (i == j) temp = half*temp
681 vquad = vquad + temp*hq(ih)
682 END DO
683 END DO
684 DO k = 1, npt
685 vquad = vquad + pq(k)*w(k)
686 END DO
687 diff = f - fopt - vquad
688 diffc = diffb
689 diffb = diffa
690 diffa = abs(diff)
691 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 fsave = fopt
698 IF (f < fopt) THEN
699 fopt = f
700 xoptsq = zero
701 DO i = 1, n
702 xopt(i) = xnew(i)
703 xoptsq = xoptsq + xopt(i)**2
704 END DO
705 END IF
706 ksave = knew
707 IF (knew > 0) GOTO 410
708 !
709 ! Pick the next value of DELTA after a trust region step.
710 !
711 IF (vquad >= zero) THEN
712 ! Return because a trust region step has failed to reduce Q
713 opt%state = 4
714 CALL get_state
715 GOTO 530
716 END IF
717 ratio = (f - fsave)/vquad
718 IF (ratio <= tenth) THEN
719 delta = half*dnorm
720 ELSE IF (ratio <= 0.7_dp) THEN
721 delta = max(half*delta, dnorm)
722 ELSE
723 delta = max(half*delta, dnorm + dnorm)
724 END IF
725 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 rhosq = max(tenth*delta, rho)**2
730 ktemp = 0
731 detrat = zero
732 IF (f >= fsave) THEN
733 ktemp = kopt
734 detrat = one
735 END IF
736 DO k = 1, npt
737 hdiag = zero
738 DO j = 1, nptm
739 temp = one
740 IF (j < idz) temp = -one
741 hdiag = hdiag + temp*zmat(k, j)**2
742 END DO
743 temp = abs(beta*hdiag + vlag(k)**2)
744 distsq = zero
745 DO j = 1, n
746 distsq = distsq + (xpt(k, j) - xopt(j))**2
747 END DO
748 IF (distsq > rhosq) temp = temp*(distsq/rhosq)**3
749 IF (temp > detrat .AND. k /= ktemp) THEN
750 detrat = temp
751 knew = k
752 END IF
753 END DO
754 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 !
760410 CALL update(n, npt, bmat, zmat, idz, ndim, vlag, beta, knew, w)
761 fval(knew) = f
762 ih = 0
763 DO i = 1, n
764 temp = pq(knew)*xpt(knew, i)
765 DO j = 1, i
766 ih = ih + 1
767 hq(ih) = hq(ih) + temp*xpt(knew, j)
768 END DO
769 END DO
770 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 DO j = 1, nptm
776 temp = diff*zmat(knew, j)
777 IF (j < idz) temp = -temp
778 DO k = 1, npt
779 pq(k) = pq(k) + temp*zmat(k, j)
780 END DO
781 END DO
782 gqsq = zero
783 DO i = 1, n
784 gq(i) = gq(i) + diff*bmat(knew, i)
785 gqsq = gqsq + gq(i)**2
786 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 IF (ksave == 0 .AND. delta == rho) THEN
794 IF (abs(ratio) > 1.0e-2_dp) THEN
795 itest = 0
796 ELSE
797 DO k = 1, npt
798 vlag(k) = fval(k) - fval(kopt)
799 END DO
800 gisq = zero
801 DO i = 1, n
802 sum = zero
803 DO k = 1, npt
804 sum = sum + bmat(k, i)*vlag(k)
805 END DO
806 gisq = gisq + sum*sum
807 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 itest = itest + 1
814 IF (gqsq < 1.0e2_dp*gisq) itest = 0
815 IF (itest >= 3) THEN
816 DO i = 1, n
817 gq(i) = w(i)
818 END DO
819 DO ih = 1, nh
820 hq(ih) = zero
821 END DO
822 DO j = 1, nptm
823 w(j) = zero
824 DO k = 1, npt
825 w(j) = w(j) + vlag(k)*zmat(k, j)
826 END DO
827 IF (j < idz) w(j) = -w(j)
828 END DO
829 DO k = 1, npt
830 pq(k) = zero
831 DO j = 1, nptm
832 pq(k) = pq(k) + zmat(k, j)*w(j)
833 END DO
834 END DO
835 itest = 0
836 END IF
837 END IF
838 END IF
839 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 IF (f <= fsave + tenth*vquad) GOTO 100
846 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 knew = 0
852460 distsq = 4.0_dp*delta*delta
853 DO k = 1, npt
854 sum = zero
855 DO j = 1, n
856 sum = sum + (xpt(k, j) - xopt(j))**2
857 END DO
858 IF (sum > distsq) THEN
859 knew = k
860 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 IF (knew > 0) THEN
868 dstep = max(min(tenth*sqrt(distsq), half*delta), rho)
869 dsq = dstep*dstep
870 GOTO 120
871 END IF
872 IF (ratio > zero) GOTO 100
873 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 !
878490 IF (rho > rhoend) THEN
879 delta = half*rho
880 ratio = rho/rhoend
881 IF (ratio <= 16.0_dp) THEN
882 rho = rhoend
883 ELSE IF (ratio <= 250.0_dp) THEN
884 rho = sqrt(ratio)*rhoend
885 ELSE
886 rho = tenth*rho
887 END IF
888 delta = max(delta, rho)
889 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 IF (knew == -1) GOTO 290
896 opt%state = 7
897 CALL get_state
898530 IF (fopt <= f) THEN
899 DO i = 1, n
900 x(i) = xbase(i) + xopt(i)
901 END DO
902 f = fopt
903 END IF
904
905 CALL get_state
906
907 !******************************************************************************
908 CONTAINS
909 !******************************************************************************
910! **************************************************************************************************
911!> \brief ...
912! **************************************************************************************************
913 SUBROUTINE get_state()
914 opt%np = np
915 opt%nh = nh
916 opt%nptm = nptm
917 opt%nftest = nftest
918 opt%idz = idz
919 opt%itest = itest
920 opt%nf = nf
921 opt%nfm = nfm
922 opt%nfmm = nfmm
923 opt%nfsav = nfsav
924 opt%knew = knew
925 opt%kopt = kopt
926 opt%ksave = ksave
927 opt%ktemp = ktemp
928 opt%rhosq = rhosq
929 opt%recip = recip
930 opt%reciq = reciq
931 opt%fbeg = fbeg
932 opt%fopt = fopt
933 opt%diffa = diffa
934 opt%xoptsq = xoptsq
935 opt%rho = rho
936 opt%delta = delta
937 opt%dsq = dsq
938 opt%dnorm = dnorm
939 opt%ratio = ratio
940 opt%temp = temp
941 opt%tempq = tempq
942 opt%beta = beta
943 opt%dx = dx
944 opt%vquad = vquad
945 opt%diff = diff
946 opt%diffc = diffc
947 opt%diffb = diffb
948 opt%fsave = fsave
949 opt%detrat = detrat
950 opt%hdiag = hdiag
951 opt%distsq = distsq
952 opt%gisq = gisq
953 opt%gqsq = gqsq
954 opt%f = f
955 opt%bstep = bstep
956 opt%alpha = alpha
957 opt%dstep = dstep
958 END SUBROUTINE get_state
959
960 !******************************************************************************
961! **************************************************************************************************
962!> \brief ...
963! **************************************************************************************************
964 SUBROUTINE set_state()
965 np = opt%np
966 nh = opt%nh
967 nptm = opt%nptm
968 nftest = opt%nftest
969 idz = opt%idz
970 itest = opt%itest
971 nf = opt%nf
972 nfm = opt%nfm
973 nfmm = opt%nfmm
974 nfsav = opt%nfsav
975 knew = opt%knew
976 kopt = opt%kopt
977 ksave = opt%ksave
978 ktemp = opt%ktemp
979 rhosq = opt%rhosq
980 recip = opt%recip
981 reciq = opt%reciq
982 fbeg = opt%fbeg
983 fopt = opt%fopt
984 diffa = opt%diffa
985 xoptsq = opt%xoptsq
986 rho = opt%rho
987 delta = opt%delta
988 dsq = opt%dsq
989 dnorm = opt%dnorm
990 ratio = opt%ratio
991 temp = opt%temp
992 tempq = opt%tempq
993 beta = opt%beta
994 dx = opt%dx
995 vquad = opt%vquad
996 diff = opt%diff
997 diffc = opt%diffc
998 diffb = opt%diffb
999 fsave = opt%fsave
1000 detrat = opt%detrat
1001 hdiag = opt%hdiag
1002 distsq = opt%distsq
1003 gisq = opt%gisq
1004 gqsq = opt%gqsq
1005 f = opt%f
1006 bstep = opt%bstep
1007 alpha = opt%alpha
1008 dstep = opt%dstep
1009 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 SUBROUTINE bigden(n, npt, xopt, xpt, bmat, zmat, idz, ndim, kopt, &
1036 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 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 DO k = 1, npt
1091 w(n + k) = zero
1092 END DO
1093 DO j = 1, nptm
1094 temp = zmat(knew, j)
1095 IF (j < idz) temp = -temp
1096 DO k = 1, npt
1097 w(n + k) = w(n + k) + temp*zmat(k, j)
1098 END DO
1099 END DO
1100 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 dd = zero
1108 ds = zero
1109 ss = zero
1110 xoptsq = zero
1111 DO i = 1, n
1112 dd = dd + d(i)**2
1113 s(i) = xpt(knew, i) - xopt(i)
1114 ds = ds + d(i)*s(i)
1115 ss = ss + s(i)**2
1116 xoptsq = xoptsq + xopt(i)**2
1117 END DO
1118 IF (ds*ds > 0.99_dp*dd*ss) THEN
1119 ksav = knew
1120 dtest = ds*ds/ss
1121 DO k = 1, npt
1122 IF (k /= kopt) THEN
1123 dstemp = zero
1124 sstemp = zero
1125 DO i = 1, n
1126 diff = xpt(k, i) - xopt(i)
1127 dstemp = dstemp + d(i)*diff
1128 sstemp = sstemp + diff*diff
1129 END DO
1130 IF (dstemp*dstemp/sstemp < dtest) THEN
1131 ksav = k
1132 dtest = dstemp*dstemp/sstemp
1133 ds = dstemp
1134 ss = sstemp
1135 END IF
1136 END IF
1137 END DO
1138 DO i = 1, n
1139 s(i) = xpt(ksav, i) - xopt(i)
1140 END DO
1141 END IF
1142 ssden = dd*ss - ds*ds
1143 iterc = 0
1144 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 iterc = iterc + 1
1151 temp = one/sqrt(ssden)
1152 xoptd = zero
1153 xopts = zero
1154 DO i = 1, n
1155 s(i) = temp*(dd*s(i) - ds*d(i))
1156 xoptd = xoptd + xopt(i)*d(i)
1157 xopts = xopts + xopt(i)*s(i)
1158 END DO
1159 !
1160 ! Set the coefficients of the first two terms of BETA.
1161 !
1162 tempa = half*xoptd*xoptd
1163 tempb = half*xopts*xopts
1164 den(1) = dd*(xoptsq + half*dd) + tempa + tempb
1165 den(2) = two*xoptd*dd
1166 den(3) = two*xopts*dd
1167 den(4) = tempa - tempb
1168 den(5) = xoptd*xopts
1169 DO i = 6, 9
1170 den(i) = zero
1171 END DO
1172 !
1173 ! Put the coefficients of Wcheck in WVEC.
1174 !
1175 DO k = 1, npt
1176 tempa = zero
1177 tempb = zero
1178 tempc = zero
1179 DO i = 1, n
1180 tempa = tempa + xpt(k, i)*d(i)
1181 tempb = tempb + xpt(k, i)*s(i)
1182 tempc = tempc + xpt(k, i)*xopt(i)
1183 END DO
1184 wvec(k, 1) = quart*(tempa*tempa + tempb*tempb)
1185 wvec(k, 2) = tempa*tempc
1186 wvec(k, 3) = tempb*tempc
1187 wvec(k, 4) = quart*(tempa*tempa - tempb*tempb)
1188 wvec(k, 5) = half*tempa*tempb
1189 END DO
1190 DO i = 1, n
1191 ip = i + npt
1192 wvec(ip, 1) = zero
1193 wvec(ip, 2) = d(i)
1194 wvec(ip, 3) = s(i)
1195 wvec(ip, 4) = zero
1196 wvec(ip, 5) = zero
1197 END DO
1198 !
1199 ! Put the coefficients of THETA*Wcheck in PROD.
1200 !
1201 DO jc = 1, 5
1202 nw = npt
1203 IF (jc == 2 .OR. jc == 3) nw = ndim
1204 DO k = 1, npt
1205 prod(k, jc) = zero
1206 END DO
1207 DO j = 1, nptm
1208 sum = zero
1209 DO k = 1, npt
1210 sum = sum + zmat(k, j)*wvec(k, jc)
1211 END DO
1212 IF (j < idz) sum = -sum
1213 DO k = 1, npt
1214 prod(k, jc) = prod(k, jc) + sum*zmat(k, j)
1215 END DO
1216 END DO
1217 IF (nw == ndim) THEN
1218 DO k = 1, npt
1219 sum = zero
1220 DO j = 1, n
1221 sum = sum + bmat(k, j)*wvec(npt + j, jc)
1222 END DO
1223 prod(k, jc) = prod(k, jc) + sum
1224 END DO
1225 END IF
1226 DO j = 1, n
1227 sum = zero
1228 DO i = 1, nw
1229 sum = sum + bmat(i, j)*wvec(i, jc)
1230 END DO
1231 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 DO k = 1, ndim
1238 sum = zero
1239 DO i = 1, 5
1240 par(i) = half*prod(k, i)*wvec(k, i)
1241 sum = sum + par(i)
1242 END DO
1243 den(1) = den(1) - par(1) - sum
1244 tempa = prod(k, 1)*wvec(k, 2) + prod(k, 2)*wvec(k, 1)
1245 tempb = prod(k, 2)*wvec(k, 4) + prod(k, 4)*wvec(k, 2)
1246 tempc = prod(k, 3)*wvec(k, 5) + prod(k, 5)*wvec(k, 3)
1247 den(2) = den(2) - tempa - half*(tempb + tempc)
1248 den(6) = den(6) - half*(tempb - tempc)
1249 tempa = prod(k, 1)*wvec(k, 3) + prod(k, 3)*wvec(k, 1)
1250 tempb = prod(k, 2)*wvec(k, 5) + prod(k, 5)*wvec(k, 2)
1251 tempc = prod(k, 3)*wvec(k, 4) + prod(k, 4)*wvec(k, 3)
1252 den(3) = den(3) - tempa - half*(tempb - tempc)
1253 den(7) = den(7) - half*(tempb + tempc)
1254 tempa = prod(k, 1)*wvec(k, 4) + prod(k, 4)*wvec(k, 1)
1255 den(4) = den(4) - tempa - par(2) + par(3)
1256 tempa = prod(k, 1)*wvec(k, 5) + prod(k, 5)*wvec(k, 1)
1257 tempb = prod(k, 2)*wvec(k, 3) + prod(k, 3)*wvec(k, 2)
1258 den(5) = den(5) - tempa - half*tempb
1259 den(8) = den(8) - par(4) + par(5)
1260 tempa = prod(k, 4)*wvec(k, 5) + prod(k, 5)*wvec(k, 4)
1261 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 DO i = 1, 5
1268 par(i) = half*prod(knew, i)**2
1269 sum = sum + par(i)
1270 END DO
1271 denex(1) = alpha*den(1) + par(1) + sum
1272 tempa = two*prod(knew, 1)*prod(knew, 2)
1273 tempb = prod(knew, 2)*prod(knew, 4)
1274 tempc = prod(knew, 3)*prod(knew, 5)
1275 denex(2) = alpha*den(2) + tempa + tempb + tempc
1276 denex(6) = alpha*den(6) + tempb - tempc
1277 tempa = two*prod(knew, 1)*prod(knew, 3)
1278 tempb = prod(knew, 2)*prod(knew, 5)
1279 tempc = prod(knew, 3)*prod(knew, 4)
1280 denex(3) = alpha*den(3) + tempa + tempb - tempc
1281 denex(7) = alpha*den(7) + tempb + tempc
1282 tempa = two*prod(knew, 1)*prod(knew, 4)
1283 denex(4) = alpha*den(4) + tempa + par(2) - par(3)
1284 tempa = two*prod(knew, 1)*prod(knew, 5)
1285 denex(5) = alpha*den(5) + tempa + prod(knew, 2)*prod(knew, 3)
1286 denex(8) = alpha*den(8) + par(4) - par(5)
1287 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 sum = denex(1) + denex(2) + denex(4) + denex(6) + denex(8)
1292 denold = sum
1293 denmax = sum
1294 isave = 0
1295 iu = 49
1296 temp = twopi/real(iu + 1, dp)
1297 par(1) = one
1298 DO i = 1, iu
1299 angle = real(i, dp)*temp
1300 par(2) = cos(angle)
1301 par(3) = sin(angle)
1302 DO j = 4, 8, 2
1303 par(j) = par(2)*par(j - 2) - par(3)*par(j - 1)
1304 par(j + 1) = par(2)*par(j - 1) + par(3)*par(j - 2)
1305 END DO
1306 sumold = sum
1307 sum = zero
1308 DO j = 1, 9
1309 sum = sum + denex(j)*par(j)
1310 END DO
1311 IF (abs(sum) > abs(denmax)) THEN
1312 denmax = sum
1313 isave = i
1314 tempa = sumold
1315 ELSE IF (i == isave + 1) THEN
1316 tempb = sum
1317 END IF
1318 END DO
1319 IF (isave == 0) tempa = sum
1320 IF (isave == iu) tempb = denold
1321 step = zero
1322 IF (tempa /= tempb) THEN
1323 tempa = tempa - denmax
1324 tempb = tempb - denmax
1325 step = half*(tempa - tempb)/(tempa + tempb)
1326 END IF
1327 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 par(2) = cos(angle)
1333 par(3) = sin(angle)
1334 DO j = 4, 8, 2
1335 par(j) = par(2)*par(j - 2) - par(3)*par(j - 1)
1336 par(j + 1) = par(2)*par(j - 1) + par(3)*par(j - 2)
1337 END DO
1338 beta = zero
1339 denmax = zero
1340 DO j = 1, 9
1341 beta = beta + den(j)*par(j)
1342 denmax = denmax + denex(j)*par(j)
1343 END DO
1344 DO k = 1, ndim
1345 vlag(k) = zero
1346 DO j = 1, 5
1347 vlag(k) = vlag(k) + prod(k, j)*par(j)
1348 END DO
1349 END DO
1350 tau = vlag(knew)
1351 dd = zero
1352 tempa = zero
1353 tempb = zero
1354 DO i = 1, n
1355 d(i) = par(2)*d(i) + par(3)*s(i)
1356 w(i) = xopt(i) + d(i)
1357 dd = dd + d(i)**2
1358 tempa = tempa + d(i)*w(i)
1359 tempb = tempb + w(i)*w(i)
1360 END DO
1361 IF (iterc >= n) EXIT mainloop
1362 IF (iterc >= 1) densav = max(densav, denold)
1363 IF (abs(denmax) <= 1.1_dp*abs(densav)) EXIT mainloop
1364 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 DO i = 1, n
1370 temp = tempa*xopt(i) + tempb*d(i) - vlag(npt + i)
1371 s(i) = tau*bmat(knew, i) + alpha*temp
1372 END DO
1373 DO k = 1, npt
1374 sum = zero
1375 DO j = 1, n
1376 sum = sum + xpt(k, j)*w(j)
1377 END DO
1378 temp = (tau*w(n + k) - alpha*vlag(k))*sum
1379 DO i = 1, n
1380 s(i) = s(i) + temp*xpt(k, i)
1381 END DO
1382 END DO
1383 ss = zero
1384 ds = zero
1385 DO i = 1, n
1386 ss = ss + s(i)**2
1387 ds = ds + d(i)*s(i)
1388 END DO
1389 ssden = dd*ss - ds*ds
1390 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 DO k = 1, ndim
1396 w(k) = zero
1397 DO j = 1, 5
1398 w(k) = w(k) + wvec(k, j)*par(j)
1399 END DO
1400 END DO
1401 vlag(kopt) = vlag(kopt) + one
1402
1403 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 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 delsq = delta*delta
1469 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 iterc = 0
1475 DO k = 1, npt
1476 hcol(k) = zero
1477 END DO
1478 DO j = 1, nptm
1479 temp = zmat(knew, j)
1480 IF (j < idz) temp = -temp
1481 DO k = 1, npt
1482 hcol(k) = hcol(k) + temp*zmat(k, j)
1483 END DO
1484 END DO
1485 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 dd = zero
1491 DO i = 1, n
1492 d(i) = xpt(knew, i) - xopt(i)
1493 gc(i) = bmat(knew, i)
1494 gd(i) = zero
1495 dd = dd + d(i)**2
1496 END DO
1497 DO k = 1, npt
1498 temp = zero
1499 sum = zero
1500 DO j = 1, n
1501 temp = temp + xpt(k, j)*xopt(j)
1502 sum = sum + xpt(k, j)*d(j)
1503 END DO
1504 temp = hcol(k)*temp
1505 sum = hcol(k)*sum
1506 DO i = 1, n
1507 gc(i) = gc(i) + temp*xpt(k, i)
1508 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 DO i = 1, n
1519 gg = gg + gc(i)**2
1520 sp = sp + d(i)*gc(i)
1521 dhd = dhd + d(i)*gd(i)
1522 END DO
1523 scale = delta/sqrt(dd)
1524 IF (sp*dhd < zero) scale = -scale
1525 temp = zero
1526 IF (sp*sp > 0.99_dp*dd*gg) temp = one
1527 tau = scale*(abs(sp) + half*scale*abs(dhd))
1528 IF (gg*delsq < 0.01_dp*tau*tau) temp = one
1529 DO i = 1, n
1530 d(i) = scale*d(i)
1531 gd(i) = scale*gd(i)
1532 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 iterc = iterc + 1
1541 dd = zero
1542 sp = zero
1543 ss = zero
1544 DO i = 1, n
1545 dd = dd + d(i)**2
1546 sp = sp + d(i)*s(i)
1547 ss = ss + s(i)**2
1548 END DO
1549 temp = dd*ss - sp*sp
1550 IF (temp <= 1.0e-8_dp*dd*ss) EXIT mainloop
1551 denom = sqrt(temp)
1552 DO i = 1, n
1553 s(i) = (dd*s(i) - sp*d(i))/denom
1554 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 DO k = 1, npt
1561 sum = zero
1562 DO j = 1, n
1563 sum = sum + xpt(k, j)*s(j)
1564 END DO
1565 sum = hcol(k)*sum
1566 DO i = 1, n
1567 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 DO i = 1, n
1576 cf1 = cf1 + s(i)*w(i)
1577 cf2 = cf2 + d(i)*gc(i)
1578 cf3 = cf3 + s(i)*gc(i)
1579 cf4 = cf4 + d(i)*gd(i)
1580 cf5 = cf5 + s(i)*gd(i)
1581 END DO
1582 cf1 = half*cf1
1583 cf4 = half*cf4 - cf1
1584 !
1585 ! Seek the value of the angle that maximizes the modulus of TAU.
1586 !
1587 taubeg = cf1 + cf2 + cf4
1588 taumax = taubeg
1589 tauold = taubeg
1590 isave = 0
1591 iu = 49
1592 temp = twopi/real(iu + 1, dp)
1593 DO i = 1, iu
1594 angle = real(i, dp)*temp
1595 cth = cos(angle)
1596 sth = sin(angle)
1597 tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1598 IF (abs(tau) > abs(taumax)) THEN
1599 taumax = tau
1600 isave = i
1601 tempa = tauold
1602 ELSE IF (i == isave + 1) THEN
1603 tempb = tau
1604 END IF
1605 tauold = tau
1606 END DO
1607 IF (isave == 0) tempa = tau
1608 IF (isave == iu) tempb = taubeg
1609 step = zero
1610 IF (tempa /= tempb) THEN
1611 tempa = tempa - taumax
1612 tempb = tempb - taumax
1613 step = half*(tempa - tempb)/(tempa + tempb)
1614 END IF
1615 angle = temp*(real(isave, dp) + step)
1616 !
1617 ! Calculate the new D and GD. Then test for convergence.
1618 !
1619 cth = cos(angle)
1620 sth = sin(angle)
1621 tau = cf1 + (cf2 + cf4*cth)*cth + (cf3 + cf5*cth)*sth
1622 DO i = 1, n
1623 d(i) = cth*d(i) + sth*s(i)
1624 gd(i) = cth*gd(i) + sth*w(i)
1625 s(i) = gc(i) + gd(i)
1626 END DO
1627 IF (abs(tau) <= 1.1_dp*abs(taubeg)) EXIT mainloop
1628 IF (iterc >= n) EXIT mainloop
1629 END DO mainloop
1630
1631 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 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 delsq = delta*delta
1692 iterc = 0
1693 itermax = n
1694 itersw = itermax
1695 DO i = 1, n
1696 d(i) = xopt(i)
1697 END DO
1698 CALL updatehd
1699 !
1700 ! Prepare for the first line search.
1701 !
1702 qred = zero
1703 dd = zero
1704 DO i = 1, n
1705 step(i) = zero
1706 hs(i) = zero
1707 g(i) = gq(i) + hd(i)
1708 d(i) = -g(i)
1709 dd = dd + d(i)**2
1710 END DO
1711 crvmin = zero
1712 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 IF (.NOT. jump2) THEN
1724 IF (.NOT. jump1) THEN
1725 iterc = iterc + 1
1726 temp = delsq - ss
1727 bstep = temp/(ds + sqrt(ds*ds + dd*temp))
1728 CALL updatehd
1729 END IF
1730 jump1 = .false.
1731 IF (iterc <= itersw) THEN
1732 dhd = zero
1733 DO j = 1, n
1734 dhd = dhd + d(j)*hd(j)
1735 END DO
1736 !
1737 ! Update CRVMIN and set the step-length ALPHA.
1738 !
1739 alpha = bstep
1740 IF (dhd > zero) THEN
1741 temp = dhd/dd
1742 IF (iterc == 1) crvmin = temp
1743 crvmin = min(crvmin, temp)
1744 alpha = min(alpha, gg/dhd)
1745 END IF
1746 qadd = alpha*(gg - half*alpha*dhd)
1747 qred = qred + qadd
1748 !
1749 ! Update STEP and HS.
1750 !
1751 ggsav = gg
1752 gg = zero
1753 DO i = 1, n
1754 step(i) = step(i) + alpha*d(i)
1755 hs(i) = hs(i) + alpha*hd(i)
1756 gg = gg + (g(i) + hs(i))**2
1757 END DO
1758 !
1759 ! Begin another conjugate direction iteration if required.
1760 !
1761 IF (alpha < bstep) THEN
1762 IF (qadd <= 0.01_dp*qred) EXIT mainloop
1763 IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1764 IF (iterc == itermax) EXIT mainloop
1765 temp = gg/ggsav
1766 dd = zero
1767 ds = zero
1768 ss = zero
1769 DO i = 1, n
1770 d(i) = temp*d(i) - g(i) - hs(i)
1771 dd = dd + d(i)**2
1772 ds = ds + d(i)*step(i)
1773 ss = ss + step(i)**2
1774 END DO
1775 IF (ds <= zero) EXIT mainloop
1776 IF (ss < delsq) cycle mainloop
1777 END IF
1778 crvmin = zero
1779 itersw = iterc
1780 jump2 = .true.
1781 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 IF (jump2) THEN
1791 sg = zero
1792 shs = zero
1793 DO i = 1, n
1794 sg = sg + step(i)*g(i)
1795 shs = shs + step(i)*hs(i)
1796 END DO
1797 sgk = sg + shs
1798 angtest = sgk/sqrt(gg*delsq)
1799 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 iterc = iterc + 1
1805 temp = sqrt(delsq*gg - sgk*sgk)
1806 tempa = delsq/temp
1807 tempb = sgk/temp
1808 DO i = 1, n
1809 d(i) = tempa*(g(i) + hs(i)) - tempb*step(i)
1810 END DO
1811 CALL updatehd
1812 IF (iterc <= itersw) THEN
1813 jump1 = .true.
1814 cycle mainloop
1815 END IF
1816 END IF
1817 dg = zero
1818 dhd = zero
1819 dhs = zero
1820 DO i = 1, n
1821 dg = dg + d(i)*g(i)
1822 dhd = dhd + hd(i)*d(i)
1823 dhs = dhs + hd(i)*step(i)
1824 END DO
1825 !
1826 ! Seek the value of the angle that minimizes Q.
1827 !
1828 cf = half*(shs - dhd)
1829 qbeg = sg + cf
1830 qsav = qbeg
1831 qmin = qbeg
1832 isave = 0
1833 iu = 49
1834 temp = twopi/real(iu + 1, dp)
1835 DO i = 1, iu
1836 angle = real(i, dp)*temp
1837 cth = cos(angle)
1838 sth = sin(angle)
1839 qnew = (sg + cf*cth)*cth + (dg + dhs*cth)*sth
1840 IF (qnew < qmin) THEN
1841 qmin = qnew
1842 isave = i
1843 tempa = qsav
1844 ELSE IF (i == isave + 1) THEN
1845 tempb = qnew
1846 END IF
1847 qsav = qnew
1848 END DO
1849 IF (isave == zero) tempa = qnew
1850 IF (isave == iu) tempb = qbeg
1851 angle = zero
1852 IF (tempa /= tempb) THEN
1853 tempa = tempa - qmin
1854 tempb = tempb - qmin
1855 angle = half*(tempa - tempb)/(tempa + tempb)
1856 END IF
1857 angle = temp*(real(isave, dp) + angle)
1858 !
1859 ! Calculate the new STEP and HS. Then test for convergence.
1860 !
1861 cth = cos(angle)
1862 sth = sin(angle)
1863 reduc = qbeg - (sg + cf*cth)*cth - (dg + dhs*cth)*sth
1864 gg = zero
1865 DO i = 1, n
1866 step(i) = cth*step(i) + sth*d(i)
1867 hs(i) = cth*hs(i) + sth*hd(i)
1868 gg = gg + (g(i) + hs(i))**2
1869 END DO
1870 qred = qred + reduc
1871 ratio = reduc/qred
1872 IF (iterc < itermax .AND. ratio > 0.01_dp) THEN
1873 jump2 = .true.
1874 ELSE
1875 EXIT mainloop
1876 END IF
1877
1878 IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop
1879
1880 END DO mainloop
1881
1882 !*******************************************************************************
1883 CONTAINS
1884! **************************************************************************************************
1885!> \brief ...
1886! **************************************************************************************************
1887 SUBROUTINE updatehd
1888 INTEGER :: i, ih, j, k
1889
1890 DO i = 1, n
1891 hd(i) = zero
1892 END DO
1893 DO k = 1, npt
1894 temp = zero
1895 DO j = 1, n
1896 temp = temp + xpt(k, j)*d(j)
1897 END DO
1898 temp = temp*pq(k)
1899 DO i = 1, n
1900 hd(i) = hd(i) + temp*xpt(k, i)
1901 END DO
1902 END DO
1903 ih = 0
1904 DO j = 1, n
1905 DO i = 1, j
1906 ih = ih + 1
1907 IF (i < j) hd(j) = hd(j) + hq(ih)*d(i)
1908 hd(i) = hd(i) + hq(ih)*d(j)
1909 END DO
1910 END DO
1911 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 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 nptm = npt - n - 1
1954 !
1955 ! Apply the rotations that put zeros in the KNEW-th row of ZMAT.
1956 !
1957 jl = 1
1958 DO j = 2, nptm
1959 IF (j == idz) THEN
1960 jl = idz
1961 ELSE IF (zmat(knew, j) /= zero) THEN
1962 temp = sqrt(zmat(knew, jl)**2 + zmat(knew, j)**2)
1963 tempa = zmat(knew, jl)/temp
1964 tempb = zmat(knew, j)/temp
1965 DO i = 1, npt
1966 temp = tempa*zmat(i, jl) + tempb*zmat(i, j)
1967 zmat(i, j) = tempa*zmat(i, j) - tempb*zmat(i, jl)
1968 zmat(i, jl) = temp
1969 END DO
1970 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 tempa = zmat(knew, 1)
1978 IF (idz >= 2) tempa = -tempa
1979 IF (jl > 1) tempb = zmat(knew, jl)
1980 DO i = 1, npt
1981 w(i) = tempa*zmat(i, 1)
1982 IF (jl > 1) w(i) = w(i) + tempb*zmat(i, jl)
1983 END DO
1984 alpha = w(knew)
1985 tau = vlag(knew)
1986 tausq = tau*tau
1987 denom = alpha*beta + tausq
1988 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 iflag = 0
1995 IF (jl == 1) THEN
1996 temp = sqrt(abs(denom))
1997 tempb = tempa/temp
1998 tempa = tau/temp
1999 DO i = 1, npt
2000 zmat(i, 1) = tempa*zmat(i, 1) - tempb*vlag(i)
2001 END DO
2002 IF (idz == 1 .AND. temp < zero) idz = 2
2003 IF (idz >= 2 .AND. temp >= zero) iflag = 1
2004 ELSE
2005 !
2006 ! Complete the updating of ZMAT in the alternative case.
2007 !
2008 ja = 1
2009 IF (beta >= zero) ja = jl
2010 jb = jl + 1 - ja
2011 temp = zmat(knew, jb)/denom
2012 tempa = temp*beta
2013 tempb = temp*tau
2014 temp = zmat(knew, ja)
2015 scala = one/sqrt(abs(beta)*temp*temp + tausq)
2016 scalb = scala*sqrt(abs(denom))
2017 DO i = 1, npt
2018 zmat(i, ja) = scala*(tau*zmat(i, ja) - temp*vlag(i))
2019 zmat(i, jb) = scalb*(zmat(i, jb) - tempa*w(i) - tempb*vlag(i))
2020 END DO
2021 IF (denom <= zero) THEN
2022 IF (beta < zero) idz = idz + 1
2023 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 idz = idz - 1
2032 DO i = 1, npt
2033 temp = zmat(i, 1)
2034 zmat(i, 1) = zmat(i, idz)
2035 zmat(i, idz) = temp
2036 END DO
2037 END IF
2038 !
2039 ! Finally, update the matrix BMAT.
2040 !
2041 DO j = 1, n
2042 jp = npt + j
2043 w(jp) = bmat(knew, j)
2044 tempa = (alpha*vlag(jp) - tau*w(jp))/denom
2045 tempb = (-beta*w(jp) - tau*vlag(jp))/denom
2046 DO i = 1, jp
2047 bmat(i, j) = bmat(i, j) + tempa*vlag(i) + tempb*w(i)
2048 IF (i > npt) bmat(jp, i - npt) = bmat(i, j)
2049 END DO
2050 END DO
2051
2052 END SUBROUTINE update
2053
2054END MODULE powell
2055
2056!******************************************************************************
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:55
powell::opt_state_type
Definition powell.F:18