(git:0de0cc2)
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 !******************************************************************************
9 MODULE powell
10  USE kinds, ONLY: dp
11  USE mathconstants, ONLY: twopi
12 #include "../base/base_uses.f90"
13 
14  IMPLICIT NONE
15 
16  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'powell'
17 
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
42 
43 !******************************************************************************
44 
45 CONTAINS
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
370 50 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
404 70 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
465 90 nfsav = nf
466  !
467  ! Generate the next trust region step and test its length. Set KNEW
468  ! to -1 if the purpose of the next F will be to improve the model.
469  !
470 100 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  !
491 120 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  !
640 290 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
645 310 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 
659 1000 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  !
760 410 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
852 460 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  !
878 490 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
898 530 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 
2054 END MODULE powell
2055 
2056 !******************************************************************************
pure real(kind=dp) function angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Definition: dumpdcd.F:1008
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public twopi
Definition: powell.F:9
subroutine, public powell_optimize(n, x, optstate)
...
Definition: powell.F:55