comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 14:35:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions
cp2k/tools/scriptmini/scriptmini.f90
Eisuke Kawashima 40b9ef4da9 Fix typo
2020-03-12 16:48:53 +01:00

2754 lines
94 KiB
Fortran
Raw Permalink Blame History

!-----------------------------------------------------------------------
SUBROUTINE WRITE_DOCUMENTATION()
write(6,'(A)') &
" scriptmini : a fortran program to minimise a script ", &
" or external program. ", &
" ", &
" Written by Joost VandeVondele ", &
" ", &
" the script is treated as a black box, that given an ", &
" input vector x returns the function value f. ", &
" ", &
" usage : ", &
" ", &
" ./scriptmini ", &
" ", &
" inputs : ", &
" ", &
" -) the script should be called scriptmini_eval ", &
" -) the script should read the n real variables ", &
" from scriptmini_eval.in ", &
" -) the script should write the function value (1 real number) ", &
" to 'scriptmini_eval.out' ", &
" -) input of scriptmini is 'scriptmini.in' with the format: ", &
" N ", &
" rhobeg rhoend ", &
" maxfun ", &
" iprint ", &
" x[1] x[2] x[3] ... x[N] ", &
" ", &
" where: ", &
" N : integer : is the number of variables ", &
" rhobeg : real : initial trust region radius, +- 10% of the ", &
" largest expected change in the variables ", &
" rhoend : real : final trust region radius, +- the final ", &
" uncertainty in the variables ", &
" maxfun : integer : the maximum number of calls to ", &
" scriptmini_eval [O(10*N**2)] ", &
" iprint : integer : output level (0-3) ", &
" 0 : no output at all (!) ", &
" ... ", &
" 3 : info at every step (recommended) ", &
" x[...] : real : the initial values of the variables "
END SUBROUTINE WRITE_DOCUMENTATION
!-------------------------------------------------------------------------------
MODULE Powell_Optimize
! Code converted using TO_F90 by Alan Miller
! Date: 2002-11-09 Time: 16:58:08
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
PRIVATE
PUBLIC :: uobyqa
CONTAINS
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uobyqa.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE uobyqa(n, x, rhobeg, rhoend, iprint, maxfun)
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: x(:)
REAL (dp), INTENT(IN) :: rhobeg
REAL (dp), INTENT(IN) :: rhoend
INTEGER, INTENT(IN) :: iprint
INTEGER, INTENT(IN) :: maxfun
! This subroutine seeks the least value of a function of many variables,
! by a trust region method that forms quadratic models by interpolation.
! The algorithm is described in "UOBYQA: unconstrained optimization by
! quadratic approximation" by M.J.D. Powell, Report DAMTP 2000/NA14,
! University of Cambridge. The arguments of the subroutine are as follows.
! N must be set to the number of variables and must be at least two.
! Initial values of the variables must be set in X(1),X(2),...,X(N). They
! will be changed to the values that give the least calculated F.
! RHOBEG and RHOEND must be set to the initial and final values of a trust
! region radius, so both must be positive with RHOEND<=RHOBEG. Typically
! RHOBEG should be about one tenth of the greatest expected change to a
! variable, and RHOEND should indicate the accuracy that is required in
! the final values of the variables.
! The value of IPRINT should be set to 0, 1, 2 or 3, which controls the
! amount of printing. Specifically, there is no output if IPRINT=0 and
! there is output only at the return if IPRINT=1. Otherwise, each new
! value of RHO is printed, with the best vector of variables so far and
! the corresponding value of the objective function. Further, each new
! value of F with its variables are output if IPRINT=3.
! MAXFUN must be set to an upper bound on the number of calls of CALFUN.
! The array W will be used for working space. Its length must be at least
! ( N**4 + 8*N**3 + 23*N**2 + 42*N + max [ 2*N**2 + 4, 18*N ] ) / 4.
! SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set F to
! the value of the objective function for the variables X(1),X(2),...,X(N).
INTEGER :: npt
! Partition the working space array, so that different parts of it can be
! treated separately by the subroutine that performs the main calculation.
npt = (n*n + 3*n + 2) / 2
CALL uobyqb(n, x, rhobeg, rhoend, iprint, maxfun, npt)
RETURN
END SUBROUTINE uobyqa
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uobyqb.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE uobyqb(n, x, rhobeg, rhoend, iprint, maxfun, npt)
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: x(:)
REAL (dp), INTENT(IN) :: rhobeg
REAL (dp), INTENT(IN) :: rhoend
INTEGER, INTENT(IN) :: iprint
INTEGER, INTENT(IN) :: maxfun
INTEGER, INTENT(IN) :: npt
INTERFACE
SUBROUTINE calfun(n, x, f)
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN) :: x(:)
REAL (dp), INTENT(OUT) :: f
END SUBROUTINE calfun
END INTERFACE
! The following arrays were previously passed as arguments:
REAL (dp) :: xbase(n), xopt(n), xnew(n), xpt(npt,n), pq(npt-1)
REAL (dp) :: pl(npt,npt-1), h(n,n), g(n), d(n), vlag(npt), w(npt)
! The arguments N, X, RHOBEG, RHOEND, IPRINT and MAXFUN are identical to
! the corresponding arguments in SUBROUTINE UOBYQA.
! NPT is set by UOBYQA to (N*N+3*N+2)/2 for the above dimension statement.
! XBASE will contain a shift of origin that reduces the contributions from
! rounding errors to values of the model and Lagrange functions.
! XOPT will be set to the displacement from XBASE of the vector of
! variables that provides the least calculated F so far.
! XNEW will be set to the displacement from XBASE of the vector of
! variables for the current calculation of F.
! XPT will contain the interpolation point coordinates relative to XBASE.
! PQ will contain the parameters of the quadratic model.
! PL will contain the parameters of the Lagrange functions.
! H will provide the second derivatives that TRSTEP and LAGMAX require.
! G will provide the first derivatives that TRSTEP and LAGMAX require.
! D is reserved for trial steps from XOPT, except that it will contain
! diagonal second derivatives during the initialization procedure.
! VLAG will contain the values of the Lagrange functions at a new point X.
! The array W will be used for working space.
REAL (dp) :: half = 0.5_dp, one = 1.0_dp, tol = 0.01_dp, two = 2.0_dp
REAL (dp) :: zero = 0.0_dp
REAL (dp) :: ddknew, delta, detrat, diff, distest, dnorm, errtol, estim
REAL (dp) :: evalue, f, fbase, fopt, fsave, ratio, rho, rhosq, sixthm
REAL (dp) :: sum, sumg, sumh, temp, tempa, tworsq, vmax, vquad, wmult
INTEGER :: i, ih, ip, iq, iw, j, jswitch, k, knew, kopt, ksave, ktemp
INTEGER :: nf, nftest, nnp, nptm
! Set some constants.
nnp = n + n + 1
nptm = npt - 1
nftest = MAX(maxfun,1)
! Initialization. NF is the number of function calculations so far.
rho = rhobeg
rhosq = rho * rho
nf = 0
DO i = 1, n
xbase(i) = x(i)
xpt(1:npt,i) = zero
END DO
pl(1:npt,1:nptm) = zero
! The branch to label 120 obtains a new value of the objective function
! and then there is a branch back to label 50, because the new function
! value is needed to form the initial quadratic model. The least function
! value so far and its index are noted below.
50 x(1:n) = xbase(1:n) + xpt(nf+1,1:n)
GO TO 150
70 IF (nf == 1) THEN
fopt = f
kopt = nf
fbase = f
j = 0
jswitch = -1
ih = n
ELSE
IF (f < fopt) THEN
fopt = f
kopt = nf
END IF
END IF
! Form the gradient and diagonal second derivatives of the initial
! quadratic model and Lagrange functions.
IF (nf <= nnp) THEN
jswitch = -jswitch
IF (jswitch > 0) THEN
IF (j >= 1) THEN
ih = ih + j
IF (w(j) < zero) THEN
d(j) = (fsave+f-two*fbase) / rhosq
pq(j) = (fsave-f) / (two*rho)
pl(1,ih) = -two / rhosq
pl(nf-1,j) = half / rho
pl(nf-1,ih) = one / rhosq
ELSE
pq(j) = (4.0D0*fsave-3.0D0*fbase-f) / (two*rho)
d(j) = (fbase+f-two*fsave) / rhosq
pl(1,j) = -1.5D0 / rho
pl(1,ih) = one / rhosq
pl(nf-1,j) = two / rho
pl(nf-1,ih) = -two / rhosq
END IF
pq(ih) = d(j)
pl(nf,j) = -half / rho
pl(nf,ih) = one / rhosq
END IF
! Pick the shift from XBASE to the next initial interpolation point
! that provides diagonal second derivatives.
IF (j < n) THEN
j = j + 1
xpt(nf+1,j) = rho
END IF
ELSE
fsave = f
IF (f < fbase) THEN
w(j) = rho
xpt(nf+1,j) = two * rho
ELSE
w(j) = -rho
xpt(nf+1,j) = -rho
END IF
END IF
IF (nf < nnp) GO TO 50
! Form the off-diagonal second derivatives of the initial quadratic model.
ih = n
ip = 1
iq = 2
END IF
ih = ih + 1
IF (nf > nnp) THEN
temp = one / (w(ip)*w(iq))
tempa = f - fbase - w(ip) * pq(ip) - w(iq) * pq(iq)
pq(ih) = (tempa - half*rhosq*(d(ip)+d(iq))) * temp
pl(1,ih) = temp
iw = ip + ip
IF (w(ip) < zero) iw = iw + 1
pl(iw,ih) = -temp
iw = iq + iq
IF (w(iq) < zero) iw = iw + 1
pl(iw,ih) = -temp
pl(nf,ih) = temp
! Pick the shift from XBASE to the next initial interpolation point
! that provides off-diagonal second derivatives.
ip = ip + 1
END IF
IF (ip == iq) THEN
ih = ih + 1
ip = 1
iq = iq + 1
END IF
IF (nf < npt) THEN
xpt(nf+1,ip) = w(ip)
xpt(nf+1,iq) = w(iq)
GO TO 50
END IF
! Set parameters to begin the iterations for the current RHO.
sixthm = zero
delta = rho
80 tworsq = (two*rho) ** 2
rhosq = rho * rho
! Form the gradient of the quadratic model at the trust region centre.
90 knew = 0
ih = n
DO j = 1, n
xopt(j) = xpt(kopt,j)
g(j) = pq(j)
DO i = 1, j
ih = ih + 1
g(i) = g(i) + pq(ih) * xopt(j)
IF (i < j) g(j) = g(j) + pq(ih) * xopt(i)
h(i,j) = pq(ih)
END DO
END DO
! Generate the next trust region step and test its length. Set KNEW
! to -1 if the purpose of the next F will be to improve conditioning,
! and also calculate a lower bound on the Hessian term of the model Q.
CALL trstep(n, g, h, delta, tol, d, evalue)
temp = zero
DO i = 1, n
temp = temp + d(i)**2
END DO
dnorm = MIN(delta,SQRT(temp))
errtol = -one
IF (dnorm < half*rho) THEN
knew = -1
errtol = half * evalue * rho * rho
IF (nf <= npt+9) errtol = zero
GO TO 290
END IF
! Calculate the next value of the objective function.
130 DO i = 1, n
xnew(i) = xopt(i) + d(i)
x(i) = xbase(i) + xnew(i)
END DO
150 IF (nf >= nftest) THEN
IF (iprint > 0) WRITE(*, 5000)
GO TO 420
END IF
nf = nf + 1
CALL calfun(n, x, f)
IF (iprint == 3) THEN
WRITE(*, 5100) nf, f, x(1:n)
END IF
IF (nf <= npt) GO TO 70
IF (knew == -1) GO TO 420
! Use the quadratic model to predict the change in F due to the step D,
! and find the values of the Lagrange functions at the new point.
vquad = zero
ih = n
DO j = 1, n
w(j) = d(j)
vquad = vquad + w(j) * pq(j)
DO i = 1, j
ih = ih + 1
w(ih) = d(i) * xnew(j) + d(j) * xopt(i)
IF (i == j) w(ih) = half * w(ih)
vquad = vquad + w(ih) * pq(ih)
END DO
END DO
DO k = 1, npt
temp = zero
DO j = 1, nptm
temp = temp + w(j) * pl(k,j)
END DO
vlag(k) = temp
END DO
vlag(kopt) = vlag(kopt) + one
! Update SIXTHM, which is a lower bound on one sixth of the greatest
! third derivative of F.
diff = f - fopt - vquad
sum = zero
DO k = 1, npt
temp = zero
DO i = 1, n
temp = temp + (xpt(k,i)-xnew(i)) ** 2
END DO
temp = SQRT(temp)
sum = sum + ABS(temp*temp*temp*vlag(k))
END DO
sixthm = MAX(sixthm, ABS(diff)/sum)
! Update FOPT and XOPT if the new F is the least value of the objective
! function so far. Then branch if D is not a trust region step.
fsave = fopt
IF (f < fopt) THEN
fopt = f
xopt(1:n) = xnew(1:n)
END IF
ksave = knew
IF (knew <= 0) THEN
! Pick the next value of DELTA after a trust region step.
IF (vquad >= zero) THEN
IF (iprint > 0) WRITE(*, 5200)
GO TO 420
END IF
ratio = (f-fsave) / vquad
IF (ratio <= 0.1D0) THEN
delta = half * dnorm
ELSE IF (ratio <= 0.7D0) THEN
delta = MAX(half*delta,dnorm)
ELSE
delta = MAX(delta, 1.25D0*dnorm, dnorm+rho)
END IF
IF (delta <= 1.5D0*rho) delta = rho
! Set KNEW to the index of the next interpolation point to be deleted.
ktemp = 0
detrat = zero
IF (f >= fsave) THEN
ktemp = kopt
detrat = one
END IF
DO k = 1, npt
sum = zero
DO i = 1, n
sum = sum + (xpt(k,i)-xopt(i)) ** 2
END DO
temp = ABS(vlag(k))
IF (sum > rhosq) temp = temp * (sum/rhosq) ** 1.5D0
IF (temp > detrat .AND. k /= ktemp) THEN
detrat = temp
ddknew = sum
knew = k
END IF
END DO
IF (knew == 0) GO TO 290
END IF
! Replace the interpolation point that has index KNEW by the point XNEW,
! and also update the Lagrange functions and the quadratic model.
DO i = 1, n
xpt(knew,i) = xnew(i)
END DO
temp = one / vlag(knew)
DO j = 1, nptm
pl(knew,j) = temp * pl(knew,j)
pq(j) = pq(j) + diff * pl(knew,j)
END DO
DO k = 1, npt
IF (k /= knew) THEN
temp = vlag(k)
DO j = 1, nptm
pl(k,j) = pl(k,j) - temp * pl(knew,j)
END DO
END IF
END DO
! Update KOPT if F is the least calculated value of the objective function.
! Then branch for another trust region calculation. The case KSAVE > 0
! indicates that a model step has just been taken.
IF (f < fsave) THEN
kopt = knew
GO TO 90
END IF
IF (ksave > 0) GO TO 90
IF (dnorm > two*rho) GO TO 90
IF (ddknew > tworsq) GO TO 90
! Alternatively, find out if the interpolation points are close
! enough to the best point so far.
290 DO k = 1, npt
w(k) = zero
DO i = 1, n
w(k) = w(k) + (xpt(k,i)-xopt(i)) ** 2
END DO
END DO
320 knew = -1
distest = tworsq
DO k = 1, npt
IF (w(k) > distest) THEN
knew = k
distest = w(k)
END IF
END DO
! If a point is sufficiently far away, then set the gradient and Hessian
! of its Lagrange function at the centre of the trust region, and find
! half the sum of squares of components of the Hessian.
IF (knew > 0) THEN
ih = n
sumh = zero
DO j = 1, n
g(j) = pl(knew,j)
DO i = 1, j
ih = ih + 1
temp = pl(knew,ih)
g(j) = g(j) + temp * xopt(i)
IF (i < j) THEN
g(i) = g(i) + temp * xopt(j)
sumh = sumh + temp * temp
END IF
h(i,j) = temp
END DO
sumh = sumh + half * temp * temp
END DO
! If ERRTOL is positive, test whether to replace the interpolation point
! with index KNEW, using a bound on the maximum modulus of its Lagrange
! function in the trust region.
IF (errtol > zero) THEN
w(knew) = zero
sumg = zero
DO i = 1, n
sumg = sumg + g(i) ** 2
END DO
estim = rho * (SQRT(sumg)+rho*SQRT(half*sumh))
wmult = sixthm * distest ** 1.5D0
IF (wmult*estim <= errtol) GO TO 320
END IF
! If the KNEW-th point may be replaced, then pick a D that gives a large
! value of the modulus of its Lagrange function within the trust region.
! Here the vector XNEW is used as temporary working space.
CALL lagmax(n, g, h, rho, d, xnew, vmax)
IF (errtol > zero) THEN
IF (wmult*vmax <= errtol) GO TO 320
END IF
GO TO 130
END IF
IF (dnorm > rho) GO TO 90
! Prepare to reduce RHO by shifting XBASE to the best point so far,
! and make the corresponding changes to the gradients of the Lagrange
! functions and the quadratic model.
IF (rho > rhoend) THEN
ih = n
DO j = 1, n
xbase(j) = xbase(j) + xopt(j)
DO k = 1, npt
xpt(k,j) = xpt(k,j) - xopt(j)
END DO
DO i = 1, j
ih = ih + 1
pq(i) = pq(i) + pq(ih) * xopt(j)
IF (i < j) THEN
pq(j) = pq(j) + pq(ih) * xopt(i)
DO k = 1, npt
pl(k,j) = pl(k,j) + pl(k,ih) * xopt(i)
END DO
END IF
DO k = 1, npt
pl(k,i) = pl(k,i) + pl(k,ih) * xopt(j)
END DO
END DO
END DO
! Pick the next values of RHO and DELTA.
delta = half * rho
ratio = rho / rhoend
IF (ratio <= 16.0D0) THEN
rho = rhoend
ELSE IF (ratio <= 250.0D0) THEN
rho = SQRT(ratio) * rhoend
ELSE
rho = 0.1D0 * rho
END IF
delta = MAX(delta,rho)
IF (iprint >= 2) THEN
IF (iprint >= 3) WRITE(*, 5300)
WRITE(*, 5400) rho, nf
WRITE(*, 5500) fopt, xbase(1:n)
END IF
GO TO 80
END IF
! Return from the calculation, after another Newton-Raphson step, if
! it is too short to have been tried before.
IF (errtol >= zero) GO TO 130
420 IF (fopt <= f) THEN
DO i = 1, n
x(i) = xbase(i) + xopt(i)
END DO
f = fopt
END IF
IF (iprint >= 1) THEN
WRITE(*, 5600) nf
WRITE(*, 5500) f, x(1:n)
END IF
RETURN
5000 FORMAT (/T5, 'Return from UOBYQA because CALFUN has been', &
' called MAXFUN times')
5100 FORMAT (/T5, 'Function number',i6,' F =', g18.10, &
' The corresponding X is:'/ (t3, 5g15.6))
5200 FORMAT (/T5, 'Return from UOBYQA because a trust', &
' region step has failed to reduce Q')
5300 FORMAT (' ')
5400 FORMAT (/T5, 'New RHO =', g11.4, ' Number of function values =',i6)
5500 FORMAT (T5, 'Least value of F =', g23.15, &
' The corresponding X is:'/ (t3, 5g15.6))
5600 FORMAT (/T5, 'At the return from UOBYQA', &
' Number of function values =', i6)
END SUBROUTINE uobyqb
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% trstep.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE trstep(n, g, h, delta, tol, d, evalue)
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN) :: g(:)
REAL (dp), INTENT(IN OUT) :: h(:,:)
REAL (dp), INTENT(IN) :: delta
REAL (dp), INTENT(IN) :: tol
REAL (dp), INTENT(OUT) :: d(:)
REAL (dp), INTENT(OUT) :: evalue
! N is the number of variables of a quadratic objective function, Q say.
! G is the gradient of Q at the origin.
! H is the Hessian matrix of Q. Only the upper triangular and diagonal
! parts need be set. The lower triangular part is used to store the
! elements of a Householder similarity transformation.
! DELTA is the trust region radius, and has to be positive.
! TOL is the value of a tolerance from the open interval (0,1).
! D will be set to the calculated vector of variables.
! EVALUE will be set to the least eigenvalue of H if and only if D is a
! Newton-Raphson step. Then EVALUE will be positive, but otherwise it
! will be set to zero.
! Let MAXRED be the maximum of Q(0)-Q(D) subject to ||D|| <= DELTA,
! and let ACTRED be the value of Q(0)-Q(D) that is actually calculated.
! We take the view that any D is acceptable if it has the properties
! ||D|| <= DELTA and ACTRED <= (1-TOL)*MAXRED.
! The calculation of D is done by the method of Section 2 of the paper
! by MJDP in the 1997 Dundee Numerical Analysis Conference Proceedings,
! after transforming H to tridiagonal form.
! The arrays GG, TD, TN, W, PIV and Z will be used for working space.
REAL (dp) :: gg(n), td(n), tn(n), w(n), piv(n), z(n)
REAL (dp) :: delsq, dhd, dnorm, dsq, dtg, dtz, gam, gnorm, gsq, hnorm
REAL (dp) :: par, parl, parlest, paru, paruest, phi, phil, phiu, pivksv
REAL (dp) :: pivot, posdef, scale, shfmax, shfmin, shift, slope, sum
REAL (dp) :: tdmin, temp, tempa, tempb, wsq, wwsq, wz, zsq
INTEGER :: i, iterc, j, jp, k, kp, kpp, ksav, ksave, nm
REAL (dp) :: one = 1.0_dp, two = 2.0_dp, zero = 0.0_dp
! Initialization.
delsq = delta * delta
evalue = zero
nm = n - 1
DO i = 1, n
d(i) = zero
td(i) = h(i,i)
DO j = 1, i
h(i,j) = h(j,i)
END DO
END DO
! Apply Householder transformations to obtain a tridiagonal matrix that
! is similar to H, and put the elements of the Householder vectors in
! the lower triangular part of H. Further, TD and TN will contain the
! diagonal and other nonzero elements of the tridiagonal matrix.
DO k = 1, nm
kp = k + 1
sum = zero
IF (kp < n) THEN
kpp = kp + 1
DO i = kpp, n
sum = sum + h(i,k) ** 2
END DO
END IF
IF (sum == zero) THEN
tn(k) = h(kp,k)
h(kp,k) = zero
ELSE
temp = h(kp,k)
tn(k) = SIGN(SQRT(sum+temp*temp),temp)
h(kp,k) = -sum / (temp+tn(k))
temp = SQRT(two/(sum+h(kp,k)**2))
DO i = kp, n
w(i) = temp * h(i,k)
h(i,k) = w(i)
z(i) = td(i) * w(i)
END DO
wz = zero
DO j = kp, nm
jp = j + 1
DO i = jp, n
z(i) = z(i) + h(i,j) * w(j)
z(j) = z(j) + h(i,j) * w(i)
END DO
wz = wz + w(j) * z(j)
END DO
wz = wz + w(n) * z(n)
DO j = kp, n
td(j) = td(j) + w(j) * (wz*w(j)-two*z(j))
IF (j < n) THEN
jp = j + 1
DO i = jp, n
h(i,j) = h(i,j) - w(i) * z(j) - w(j) * (z(i)-wz*w(i))
END DO
END IF
END DO
END IF
END DO
! Form GG by applying the similarity transformation to G.
gsq = zero
DO i = 1, n
gg(i) = g(i)
gsq = gsq + g(i) ** 2
END DO
gnorm = SQRT(gsq)
DO k = 1, nm
kp = k + 1
sum = zero
DO i = kp, n
sum = sum + gg(i) * h(i,k)
END DO
DO i = kp, n
gg(i) = gg(i) - sum * h(i,k)
END DO
END DO
! Begin the trust region calculation with a tridiagonal matrix by
! calculating the norm of H. Then treat the case when H is zero.
hnorm = ABS(td(1)) + ABS(tn(1))
tdmin = td(1)
tn(n) = zero
DO i = 2, n
temp = ABS(tn(i-1)) + ABS(td(i)) + ABS(tn(i))
hnorm = MAX(hnorm,temp)
tdmin = MIN(tdmin,td(i))
END DO
IF (hnorm == zero) THEN
IF (gnorm == zero) GO TO 420
scale = delta / gnorm
DO i = 1, n
d(i) = -scale * gg(i)
END DO
GO TO 380
END IF
! Set the initial values of PAR and its bounds.
parl = MAX(zero, -tdmin, gnorm/delta-hnorm)
parlest = parl
par = parl
paru = zero
paruest = zero
posdef = zero
iterc = 0
! Calculate the pivots of the Cholesky factorization of (H+PAR*I).
160 iterc = iterc + 1
ksav = 0
piv(1) = td(1) + par
k = 1
170 IF (piv(k) > zero) THEN
piv(k+1) = td(k+1) + par - tn(k) ** 2 / piv(k)
ELSE
IF (piv(k) < zero .OR. tn(k) /= zero) GO TO 180
ksav = k
piv(k+1) = td(k+1) + par
END IF
k = k + 1
IF (k < n) GO TO 170
IF (piv(k) >= zero) THEN
IF (piv(k) == zero) ksav = k
! Branch if all the pivots are positive, allowing for the case when
! G is zero.
IF (ksav == 0 .AND. gsq > zero) GO TO 250
IF (gsq == zero) THEN
IF (par == zero) GO TO 380
paru = par
paruest = par
IF (ksav == 0) GO TO 210
END IF
k = ksav
END IF
! Set D to a direction of nonpositive curvature of the given tridiagonal
! matrix, and thus revise PARLEST.
180 d(k) = one
IF (ABS(tn(k)) <= ABS(piv(k))) THEN
dsq = one
dhd = piv(k)
ELSE
temp = td(k+1) + par
IF (temp <= ABS(piv(k))) THEN
d(k+1) = SIGN(one,-tn(k))
dhd = piv(k) + temp - two * ABS(tn(k))
ELSE
d(k+1) = -tn(k) / temp
dhd = piv(k) + tn(k) * d(k+1)
END IF
dsq = one + d(k+1) ** 2
END IF
190 IF (k > 1) THEN
k = k - 1
IF (tn(k) /= zero) THEN
d(k) = -tn(k) * d(k+1) / piv(k)
dsq = dsq + d(k) ** 2
GO TO 190
END IF
d(1:k) = zero
END IF
parl = par
parlest = par - dhd / dsq
! Terminate with D set to a multiple of the current D if the following
! test suggests that it suitable to do so.
210 temp = paruest
IF (gsq == zero) temp = temp * (one-tol)
IF (paruest > zero .AND. parlest >= temp) THEN
dtg = DOT_PRODUCT( d(1:n), gg(1:n) )
scale = -SIGN(delta/SQRT(dsq),dtg)
d(1:n) = scale * d(1:n)
GO TO 380
END IF
! Pick the value of PAR for the next iteration.
240 IF (paru == zero) THEN
par = two * parlest + gnorm / delta
ELSE
par = 0.5D0 * (parl+paru)
par = MAX(par,parlest)
END IF
IF (paruest > zero) par = MIN(par,paruest)
GO TO 160
! Calculate D for the current PAR in the positive definite case.
250 w(1) = -gg(1) / piv(1)
DO i = 2, n
w(i) = (-gg(i)-tn(i-1)*w(i-1)) / piv(i)
END DO
d(n) = w(n)
DO i = nm, 1, -1
d(i) = w(i) - tn(i) * d(i+1) / piv(i)
END DO
! Branch if a Newton-Raphson step is acceptable.
dsq = zero
wsq = zero
DO i = 1, n
dsq = dsq + d(i) ** 2
wsq = wsq + piv(i) * w(i) ** 2
END DO
IF (par /= zero .OR. dsq > delsq) THEN
! Make the usual test for acceptability of a full trust region step.
dnorm = SQRT(dsq)
phi = one / dnorm - one / delta
temp = tol * (one+par*dsq/wsq) - dsq * phi * phi
IF (temp >= zero) THEN
scale = delta / dnorm
DO i = 1, n
d(i) = scale * d(i)
END DO
GO TO 380
END IF
IF (iterc >= 2 .AND. par <= parl) GO TO 380
IF (paru > zero .AND. par >= paru) GO TO 380
! Complete the iteration when PHI is negative.
IF (phi < zero) THEN
parlest = par
IF (posdef == one) THEN
IF (phi <= phil) GO TO 380
slope = (phi-phil) / (par-parl)
parlest = par - phi / slope
END IF
slope = one / gnorm
IF (paru > zero) slope = (phiu-phi) / (paru-par)
temp = par - phi / slope
IF (paruest > zero) temp = MIN(temp,paruest)
paruest = temp
posdef = one
parl = par
phil = phi
GO TO 240
END IF
! If required, calculate Z for the alternative test for convergence.
IF (posdef == zero) THEN
w(1) = one / piv(1)
DO i = 2, n
temp = -tn(i-1) * w(i-1)
w(i) = (SIGN(one,temp)+temp) / piv(i)
END DO
z(n) = w(n)
DO i = nm, 1, -1
z(i) = w(i) - tn(i) * z(i+1) / piv(i)
END DO
wwsq = zero
zsq = zero
dtz = zero
DO i = 1, n
wwsq = wwsq + piv(i) * w(i) ** 2
zsq = zsq + z(i) ** 2
dtz = dtz + d(i) * z(i)
END DO
! Apply the alternative test for convergence.
tempa = ABS(delsq-dsq)
tempb = SQRT(dtz*dtz+tempa*zsq)
gam = tempa / (SIGN(tempb,dtz)+dtz)
temp = tol * (wsq+par*delsq) - gam * gam * wwsq
IF (temp >= zero) THEN
DO i = 1, n
d(i) = d(i) + gam * z(i)
END DO
GO TO 380
END IF
parlest = MAX(parlest,par-wwsq/zsq)
END IF
! Complete the iteration when PHI is positive.
slope = one / gnorm
IF (paru > zero) THEN
IF (phi >= phiu) GO TO 380
slope = (phiu-phi) / (paru-par)
END IF
parlest = MAX(parlest,par-phi/slope)
paruest = par
IF (posdef == one) THEN
slope = (phi-phil) / (par-parl)
paruest = par - phi / slope
END IF
paru = par
phiu = phi
GO TO 240
END IF
! Set EVALUE to the least eigenvalue of the second derivative matrix if
! D is a Newton-Raphson step. SHFMAX will be an upper bound on EVALUE.
shfmin = zero
pivot = td(1)
shfmax = pivot
DO k = 2, n
pivot = td(k) - tn(k-1) ** 2 / pivot
shfmax = MIN(shfmax,pivot)
END DO
! Find EVALUE by a bisection method, but occasionally SHFMAX may be
! adjusted by the rule of false position.
ksave = 0
350 shift = 0.5D0 * (shfmin+shfmax)
k = 1
temp = td(1) - shift
360 IF (temp > zero) THEN
piv(k) = temp
IF (k < n) THEN
temp = td(k+1) - shift - tn(k) ** 2 / temp
k = k + 1
GO TO 360
END IF
shfmin = shift
ELSE
IF (k < ksave) GO TO 370
IF (k == ksave) THEN
IF (pivksv == zero) GO TO 370
IF (piv(k)-temp < temp-pivksv) THEN
pivksv = temp
shfmax = shift
ELSE
pivksv = zero
shfmax = (shift*piv(k) - shfmin*temp) / (piv(k)-temp)
END IF
ELSE
ksave = k
pivksv = temp
shfmax = shift
END IF
END IF
IF (shfmin <= 0.99D0*shfmax) GO TO 350
370 evalue = shfmin
! Apply the inverse Householder transformations to D.
380 nm = n - 1
DO k = nm, 1, -1
kp = k + 1
sum = zero
DO i = kp, n
sum = sum + d(i) * h(i,k)
END DO
DO i = kp, n
d(i) = d(i) - sum * h(i,k)
END DO
END DO
! Return from the subroutine.
420 RETURN
END SUBROUTINE trstep
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lagmax.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE lagmax(n, g, h, rho, d, v, vmax)
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN) :: g(:)
REAL (dp), INTENT(OUT) :: h(:,:)
REAL (dp), INTENT(IN) :: rho
REAL (dp), INTENT(OUT) :: d(:)
REAL (dp), INTENT(OUT) :: v(:)
REAL (dp), INTENT(OUT) :: vmax
! N is the number of variables of a quadratic objective function, Q say.
! G is the gradient of Q at the origin.
! H is the symmetric Hessian matrix of Q. Only the upper triangular and
! diagonal parts need be set.
! RHO is the trust region radius, and has to be positive.
! D will be set to the calculated vector of variables.
! The array V will be used for working space.
! VMAX will be set to |Q(0)-Q(D)|.
! Calculating the D that maximizes |Q(0)-Q(D)| subject to ||D|| <= RHO
! requires of order N**3 operations, but sometimes it is adequate if
! |Q(0)-Q(D)| is within about 0.9 of its greatest possible value. This
! subroutine provides such a solution in only of order N**2 operations,
! where the claim of accuracy has been tested by numerical experiments.
REAL (dp) :: half = 0.5_dp, one = 1.0_dp, zero = 0.0_dp
REAL (dp) :: dd, dhd, dlin, dsq, gd, gg, ghg, gnorm, halfrt, hmax, ratio
REAL (dp) :: scale, sum, sumv, temp, tempa, tempb, tempc, tempd, tempv
REAL (dp) :: vhg, vhv, vhw, vlin, vmu, vnorm, vsq, vv, wcos, whw, wsin, wsq
INTEGER :: i, j, k
! Preliminary calculations.
halfrt = SQRT(half)
! Pick V such that ||HV|| / ||V|| is large.
hmax = zero
DO i = 1, n
sum = zero
DO j = 1, n
h(j,i) = h(i,j)
sum = sum + h(i,j) ** 2
END DO
IF (sum > hmax) THEN
hmax = sum
k = i
END IF
END DO
DO j = 1, n
v(j) = h(k,j)
END DO
! Set D to a vector in the subspace spanned by V and HV that maximizes
! |(D,HD)|/(D,D), except that we set D=HV if V and HV are nearly parallel.
! The vector that has the name D at label 60 used to be the vector W.
vsq = zero
vhv = zero
dsq = zero
DO i = 1, n
vsq = vsq + v(i) ** 2
d(i) = DOT_PRODUCT( h(i,1:n), v(1:n) )
vhv = vhv + v(i) * d(i)
dsq = dsq + d(i) ** 2
END DO
IF (vhv*vhv <= 0.9999D0*dsq*vsq) THEN
temp = vhv / vsq
wsq = zero
DO i = 1, n
d(i) = d(i) - temp * v(i)
wsq = wsq + d(i) ** 2
END DO
whw = zero
ratio = SQRT(wsq/vsq)
DO i = 1, n
temp = DOT_PRODUCT( h(i,1:n), d(1:n) )
whw = whw + temp * d(i)
v(i) = ratio * v(i)
END DO
vhv = ratio * ratio * vhv
vhw = ratio * wsq
temp = half * (whw-vhv)
temp = temp + SIGN(SQRT(temp**2+vhw**2),whw+vhv)
DO i = 1, n
d(i) = vhw * v(i) + temp * d(i)
END DO
END IF
! We now turn our attention to the subspace spanned by G and D. A multiple
! of the current D is returned if that choice seems to be adequate.
gg = zero
gd = zero
dd = zero
dhd = zero
DO i = 1, n
gg = gg + g(i) ** 2
gd = gd + g(i) * d(i)
dd = dd + d(i) ** 2
sum = DOT_PRODUCT( h(i,1:n), d(1:n) )
dhd = dhd + sum * d(i)
END DO
temp = gd / gg
vv = zero
scale = SIGN(rho/SQRT(dd),gd*dhd)
DO i = 1, n
v(i) = d(i) - temp * g(i)
vv = vv + v(i) ** 2
d(i) = scale * d(i)
END DO
gnorm = SQRT(gg)
IF (gnorm*dd <= 0.5D-2*rho*ABS(dhd) .OR. vv/dd <= 1.0D-4) THEN
vmax = ABS(scale*(gd + half*scale*dhd))
GO TO 170
END IF
! G and V are now orthogonal in the subspace spanned by G and D. Hence
! we generate an orthonormal basis of this subspace such that (D,HV) is
! negligible or zero, where D and V will be the basis vectors.
ghg = zero
vhg = zero
vhv = zero
DO i = 1, n
sum = DOT_PRODUCT( h(i,1:n), g(1:n) )
sumv = DOT_PRODUCT( h(i,1:n), v(1:n) )
ghg = ghg + sum * g(i)
vhg = vhg + sumv * g(i)
vhv = vhv + sumv * v(i)
END DO
vnorm = SQRT(vv)
ghg = ghg / gg
vhg = vhg / (vnorm*gnorm)
vhv = vhv / vv
IF (ABS(vhg) <= 0.01D0*MAX(ABS(ghg),ABS(vhv))) THEN
vmu = ghg - vhv
wcos = one
wsin = zero
ELSE
temp = half * (ghg-vhv)
vmu = temp + SIGN(SQRT(temp**2+vhg**2),temp)
temp = SQRT(vmu**2+vhg**2)
wcos = vmu / temp
wsin = vhg / temp
END IF
tempa = wcos / gnorm
tempb = wsin / vnorm
tempc = wcos / vnorm
tempd = wsin / gnorm
DO i = 1, n
d(i) = tempa * g(i) + tempb * v(i)
v(i) = tempc * v(i) - tempd * g(i)
END DO
! The final D is a multiple of the current D, V, D+V or D-V. We make the
! choice from these possibilities that is optimal.
dlin = wcos * gnorm / rho
vlin = -wsin * gnorm / rho
tempa = ABS(dlin) + half * ABS(vmu+vhv)
tempb = ABS(vlin) + half * ABS(ghg-vmu)
tempc = halfrt * (ABS(dlin)+ABS(vlin)) + 0.25D0 * ABS(ghg+vhv)
IF (tempa >= tempb .AND. tempa >= tempc) THEN
tempd = SIGN(rho,dlin*(vmu+vhv))
tempv = zero
ELSE IF (tempb >= tempc) THEN
tempd = zero
tempv = SIGN(rho,vlin*(ghg-vmu))
ELSE
tempd = SIGN(halfrt*rho,dlin*(ghg+vhv))
tempv = SIGN(halfrt*rho,vlin*(ghg+vhv))
END IF
DO i = 1, n
d(i) = tempd * d(i) + tempv * v(i)
END DO
vmax = rho * rho * MAX(tempa,tempb,tempc)
170 RETURN
END SUBROUTINE lagmax
END MODULE Powell_Optimize
!-------------------------------------------------------------------------------
!
! Main program scriptmini
!
! reads input and starts optimisation
!
!-------------------------------------------------------------------------------
PROGRAM scriptmini
USE Powell_Optimize
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
REAL (dp) :: rhobeg, rhoend
REAL(dp), DIMENSION(:), ALLOCATABLE :: x
INTEGER :: iprint, maxfun, n,istat
OPEN(UNIT=17,FILE="scriptmini.in",FORM="FORMATTED",STATUS="OLD",IOSTAT=ISTAT)
IF (ISTAT.NE.0) THEN
CALL write_documentation()
STOP " Unable to open scriptmini.in "
END IF
READ(17,*,IOSTAT=ISTAT) N
IF (ISTAT.NE.0) THEN
CALL write_documentation()
STOP " Unable to read N in scriptmini.in "
END IF
ALLOCATE(x(N))
READ(17,*,IOSTAT=ISTAT) rhobeg,rhoend
IF (ISTAT.NE.0) THEN
CALL write_documentation()
STOP " Unable to read rhobeg,rhoend in scriptmini.in "
END IF
READ(17,*,IOSTAT=ISTAT) maxfun
IF (ISTAT.NE.0) THEN
CALL write_documentation()
STOP " Unable to read maxfun in scriptmini.in "
END IF
READ(17,*,IOSTAT=ISTAT) iprint
IF (ISTAT.NE.0) THEN
CALL write_documentation()
STOP " Unable to read iprint in scriptmini.in "
END IF
READ(17,*,IOSTAT=ISTAT) x
IF (ISTAT.NE.0) THEN
CALL write_documentation()
STOP " Unable to read x in scriptmini.in "
END IF
IF (.FALSE.) THEN
CALL uobyqa (n, x, rhobeg, rhoend, iprint, maxfun)
ELSE
CALL newuoa (n, x, rhobeg, rhoend, iprint, maxfun)
ENDIF
DEALLOCATE(x)
END PROGRAM scriptmini
!-------------------------------------------------------------------------------
!
! calfun: the actual evaluation of the script
!
!
!-------------------------------------------------------------------------------
SUBROUTINE calfun(n, x, f)
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN) :: x(:)
REAL (dp), INTENT(OUT) :: f
character(LEN=40) :: format
! write variables on a single line in the file
WRITE(format,'(A1,I4.4,A12)') '(',n,'(1X,F30.20))'
OPEN(UNIT=17,FILE="scriptmini_eval.in")
WRITE(UNIT=17,FMT=format) x(1:n)
CLOSE(UNIT=17)
! execute scriptmini_eval
CALL system("./scriptmini_eval")
! read value of the energy back
OPEN(UNIT=17,FILE="scriptmini_eval.out")
READ(UNIT=17,FMT=*) f
CLOSE(UNIT=17)
END SUBROUTINE calfun
!
!
!
!
!
SUBROUTINE NEWUOA (N,X,RHOBEG,RHOEND,IPRINT,MAXFUN)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(*)
REAL*8, DIMENSION(:), ALLOCATABLE :: W
NPT=2*N+1
ALLOCATE(W((NPT+13)*(NPT+N)+3*N*(N+3)/2))
!
! This subroutine seeks the least value of a function of many variab
! by a trust region method that forms quadratic models by interpolat
! There can be some freedom in the interpolation conditions, which i
! taken up by minimizing the Frobenius norm of the change to the sec
! derivative of the quadratic model, beginning with a zero matrix. T
! arguments of the subroutine are as follows.
!
! N must be set to the number of variables and must be at least two.
! NPT is the number of interpolation conditions. Its value must be i
! interval [N+2,(N+1)(N+2)/2].
! Initial values of the variables must be set in X(1),X(2),...,X(N).
! will be changed to the values that give the least calculated F.
! RHOBEG and RHOEND must be set to the initial and final values of a
! region radius, so both must be positive with RHOEND<=RHOBEG. Typ
! RHOBEG should be about one tenth of the greatest expected change
! variable, and RHOEND should indicate the accuracy that is requir
! the final values of the variables.
! The value of IPRINT should be set to 0, 1, 2 or 3, which controls
! amount of printing. Specifically, there is no output if IPRINT=0
! there is output only at the return if IPRINT=1. Otherwise, each
! value of RHO is printed, with the best vector of variables so fa
! the corresponding value of the objective function. Further, each
! value of F with its variables are output if IPRINT=3.
! MAXFUN must be set to an upper bound on the number of calls of CAL
! The array W will be used for working space. Its length must be at
! (NPT+13)*(NPT+N)+3*N*(N+3)/2.
!
! SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must se
! the value of the objective function for the variables X(1),X(2),..
!
! Partition the working space array, so that different parts of it c
! treated separately by the subroutine that performs the main calcul
!
NP=N+1
NPTM=NPT-NP
IF (NPT .LT. N+2 .OR. NPT .GT. ((N+2)*NP)/2) THEN
PRINT 10
10 FORMAT (/4X,'Return from NEWUOA because NPT is not in', &
& ' the required interval')
GO TO 20
END IF
NDIM=NPT+N
IXB=1
IXO=IXB+N
IXN=IXO+N
IXP=IXN+N
IFV=IXP+N*NPT
IGQ=IFV+NPT
IHQ=IGQ+N
IPQ=IHQ+(N*NP)/2
IBMAT=IPQ+NPT
IZMAT=IBMAT+NDIM*N
ID=IZMAT+NPT*NPTM
IVL=ID+N
IW=IVL+NDIM
!
! The above settings provide a partition of W for subroutine NEWUOB.
! The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements
! W plus the space that is needed by the last array of NEWUOB.
!
CALL NEWUOB (N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W(IXB), &
& W(IXO),W(IXN),W(IXP),W(IFV),W(IGQ),W(IHQ),W(IPQ),W(IBMAT), &
& W(IZMAT),NDIM,W(ID),W(IVL),W(IW))
20 RETURN
END
SUBROUTINE NEWUOB (N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,XBASE, &
& XOPT,XNEW,XPT,FVAL,GQ,HQ,PQ,BMAT,ZMAT,NDIM,D,VLAG,W)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(1:N),XBASE(*),XOPT(*),XNEW(*),XPT(NPT,*),FVAL(*), &
& GQ(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),D(*),VLAG(*),W(*)
!
! The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are ide
! to the corresponding arguments in SUBROUTINE NEWUOA.
! XBASE will hold a shift of origin that should reduce the contribut
! from rounding errors to values of the model and Lagrange functio
! XOPT will be set to the displacement from XBASE of the vector of
! variables that provides the least calculated F so far.
! XNEW will be set to the displacement from XBASE of the vector of
! variables for the current calculation of F.
! XPT will contain the interpolation point coordinates relative to X
! FVAL will hold the values of F at the interpolation points.
! GQ will hold the gradient of the quadratic model at XBASE.
! HQ will hold the explicit second derivatives of the quadratic mode
! PQ will contain the parameters of the implicit second derivatives
! the quadratic model.
! BMAT will hold the last N columns of H.
! ZMAT will hold the factorization of the leading NPT by NPT submatr
! H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, wh
! the elements of DZ are plus or minus one, as specified by IDZ.
! NDIM is the first dimension of BMAT and has the value NPT+N.
! D is reserved for trial steps from XOPT.
! VLAG will contain the values of the Lagrange functions at a new po
! They are part of a product that requires VLAG to be of length ND
! The array W will be used for working space. Its length must be at
! 10*NDIM = 10*(NPT+N).
!
! Set some constants.
!
INTERFACE
SUBROUTINE calfun(n, x, f)
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN) :: x(:)
REAL (dp), INTENT(OUT) :: f
END SUBROUTINE calfun
END INTERFACE
HALF=0.5D0
ONE=1.0D0
TENTH=0.1D0
ZERO=0.0D0
NP=N+1
NH=(N*NP)/2
NPTM=NPT-NP
NFTEST=MAX0(MAXFUN,1)
!
! Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero.
!
DO 20 J=1,N
XBASE(J)=X(J)
DO 10 K=1,NPT
10 XPT(K,J)=ZERO
DO 20 I=1,NDIM
20 BMAT(I,J)=ZERO
DO 30 IH=1,NH
30 HQ(IH)=ZERO
DO 40 K=1,NPT
PQ(K)=ZERO
DO 40 J=1,NPTM
40 ZMAT(K,J)=ZERO
!
! Begin the initialization procedure. NF becomes one more than the n
! of function values so far. The coordinates of the displacement of
! next initial interpolation point from XBASE are set in XPT(NF,.).
!
RHOSQ=RHOBEG*RHOBEG
RECIP=ONE/RHOSQ
RECIQ=DSQRT(HALF)/RHOSQ
NF=0
50 NFM=NF
NFMM=NF-N
NF=NF+1
IF (NFM .LE. 2*N) THEN
IF (NFM .GE. 1 .AND. NFM .LE. N) THEN
XPT(NF,NFM)=RHOBEG
ELSE IF (NFM .GT. N) THEN
XPT(NF,NFMM)=-RHOBEG
END IF
ELSE
ITEMP=(NFMM-1)/N
JPT=NFM-ITEMP*N-N
IPT=JPT+ITEMP
IF (IPT .GT. N) THEN
ITEMP=JPT
JPT=IPT-N
IPT=ITEMP
END IF
XIPT=RHOBEG
IF (FVAL(IPT+NP) .LT. FVAL(IPT+1)) XIPT=-XIPT
XJPT=RHOBEG
IF (FVAL(JPT+NP) .LT. FVAL(JPT+1)) XJPT=-XJPT
XPT(NF,IPT)=XIPT
XPT(NF,JPT)=XJPT
END IF
!
! Calculate the next value of F, label 70 being reached immediately
! after this calculation. The least function value so far and its in
! are required.
!
DO 60 J=1,N
60 X(J)=XPT(NF,J)+XBASE(J)
GOTO 310
70 FVAL(NF)=F
IF (NF .EQ. 1) THEN
FBEG=F
FOPT=F
KOPT=1
ELSE IF (F .LT. FOPT) THEN
FOPT=F
KOPT=NF
END IF
!
! Set the nonzero initial elements of BMAT and the quadratic model i
! the cases when NF is at most 2*N+1.
!
IF (NFM .LE. 2*N) THEN
IF (NFM .GE. 1 .AND. NFM .LE. N) THEN
GQ(NFM)=(F-FBEG)/RHOBEG
IF (NPT .LT. NF+N) THEN
BMAT(1,NFM)=-ONE/RHOBEG
BMAT(NF,NFM)=ONE/RHOBEG
BMAT(NPT+NFM,NFM)=-HALF*RHOSQ
END IF
ELSE IF (NFM .GT. N) THEN
BMAT(NF-N,NFMM)=HALF/RHOBEG
BMAT(NF,NFMM)=-HALF/RHOBEG
ZMAT(1,NFMM)=-RECIQ-RECIQ
ZMAT(NF-N,NFMM)=RECIQ
ZMAT(NF,NFMM)=RECIQ
IH=(NFMM*(NFMM+1))/2
TEMP=(FBEG-F)/RHOBEG
HQ(IH)=(GQ(NFMM)-TEMP)/RHOBEG
GQ(NFMM)=HALF*(GQ(NFMM)+TEMP)
END IF
!
! Set the off-diagonal second derivatives of the Lagrange functions
! the initial quadratic model.
!
ELSE
IH=(IPT*(IPT-1))/2+JPT
IF (XIPT .LT. ZERO) IPT=IPT+N
IF (XJPT .LT. ZERO) JPT=JPT+N
ZMAT(1,NFMM)=RECIP
ZMAT(NF,NFMM)=RECIP
ZMAT(IPT+1,NFMM)=-RECIP
ZMAT(JPT+1,NFMM)=-RECIP
HQ(IH)=(FBEG-FVAL(IPT+1)-FVAL(JPT+1)+F)/(XIPT*XJPT)
END IF
IF (NF .LT. NPT) GOTO 50
!
! Begin the iterative procedure, because the initial model is comple
!
RHO=RHOBEG
DELTA=RHO
IDZ=1
DIFFA=ZERO
DIFFB=ZERO
ITEST=0
XOPTSQ=ZERO
DO 80 I=1,N
XOPT(I)=XPT(KOPT,I)
80 XOPTSQ=XOPTSQ+XOPT(I)**2
90 NFSAV=NF
!
! Generate the next trust region step and test its length. Set KNEW
! to -1 if the purpose of the next F will be to improve the model.
!
100 KNEW=0
CALL TRSAPP (N,NPT,XOPT,XPT,GQ,HQ,PQ,DELTA,D,W,W(NP), &
& W(NP+N),W(NP+2*N),CRVMIN)
DSQ=ZERO
DO 110 I=1,N
110 DSQ=DSQ+D(I)**2
DNORM=DMIN1(DELTA,DSQRT(DSQ))
IF (DNORM .LT. HALF*RHO) THEN
KNEW=-1
DELTA=TENTH*DELTA
RATIO=-1.0D0
IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
IF (NF .LE. NFSAV+2) GOTO 460
TEMP=0.125D0*CRVMIN*RHO*RHO
IF (TEMP .LE. DMAX1(DIFFA,DIFFB,DIFFC)) GOTO 460
GOTO 490
END IF
!
! Shift XBASE if XOPT may be too far from XBASE. First make the chan
! to BMAT that do not depend on ZMAT.
!
120 IF (DSQ .LE. 1.0D-3*XOPTSQ) THEN
TEMPQ=0.25D0*XOPTSQ
DO 140 K=1,NPT
SUM=ZERO
DO 130 I=1,N
130 SUM=SUM+XPT(K,I)*XOPT(I)
TEMP=PQ(K)*SUM
SUM=SUM-HALF*XOPTSQ
W(NPT+K)=SUM
DO 140 I=1,N
GQ(I)=GQ(I)+TEMP*XPT(K,I)
XPT(K,I)=XPT(K,I)-HALF*XOPT(I)
VLAG(I)=BMAT(K,I)
W(I)=SUM*XPT(K,I)+TEMPQ*XOPT(I)
IP=NPT+I
DO 140 J=1,I
140 BMAT(IP,J)=BMAT(IP,J)+VLAG(I)*W(J)+W(I)*VLAG(J)
!
! Then the revisions of BMAT that depend on ZMAT are calculated.
!
DO 180 K=1,NPTM
SUMZ=ZERO
DO 150 I=1,NPT
SUMZ=SUMZ+ZMAT(I,K)
150 W(I)=W(NPT+I)*ZMAT(I,K)
DO 170 J=1,N
SUM=TEMPQ*SUMZ*XOPT(J)
DO 160 I=1,NPT
160 SUM=SUM+W(I)*XPT(I,J)
VLAG(J)=SUM
IF (K .LT. IDZ) SUM=-SUM
DO 170 I=1,NPT
170 BMAT(I,J)=BMAT(I,J)+SUM*ZMAT(I,K)
DO 180 I=1,N
IP=I+NPT
TEMP=VLAG(I)
IF (K .LT. IDZ) TEMP=-TEMP
DO 180 J=1,I
180 BMAT(IP,J)=BMAT(IP,J)+TEMP*VLAG(J)
!
! The following instructions complete the shift of XBASE, including
! the changes to the parameters of the quadratic model.
!
IH=0
DO 200 J=1,N
W(J)=ZERO
DO 190 K=1,NPT
W(J)=W(J)+PQ(K)*XPT(K,J)
190 XPT(K,J)=XPT(K,J)-HALF*XOPT(J)
DO 200 I=1,J
IH=IH+1
IF (I .LT. J) GQ(J)=GQ(J)+HQ(IH)*XOPT(I)
GQ(I)=GQ(I)+HQ(IH)*XOPT(J)
HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J)
200 BMAT(NPT+I,J)=BMAT(NPT+J,I)
DO 210 J=1,N
XBASE(J)=XBASE(J)+XOPT(J)
210 XOPT(J)=ZERO
XOPTSQ=ZERO
END IF
!
! Pick the model step if KNEW is positive. A different choice of D
! may be made later, if the choice of D by BIGLAG causes substantial
! cancellation in DENOM.
!
IF (KNEW .GT. 0) THEN
CALL BIGLAG (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KNEW,DSTEP, &
& D,ALPHA,VLAG,VLAG(NPT+1),W,W(NP),W(NP+N))
END IF
!
! Calculate VLAG and BETA for the current choice of D. The first NPT
! components of W_check will be held in W.
!
DO 230 K=1,NPT
SUMA=ZERO
SUMB=ZERO
SUM=ZERO
DO 220 J=1,N
SUMA=SUMA+XPT(K,J)*D(J)
SUMB=SUMB+XPT(K,J)*XOPT(J)
220 SUM=SUM+BMAT(K,J)*D(J)
W(K)=SUMA*(HALF*SUMA+SUMB)
230 VLAG(K)=SUM
BETA=ZERO
DO 250 K=1,NPTM
SUM=ZERO
DO 240 I=1,NPT
240 SUM=SUM+ZMAT(I,K)*W(I)
IF (K .LT. IDZ) THEN
BETA=BETA+SUM*SUM
SUM=-SUM
ELSE
BETA=BETA-SUM*SUM
END IF
DO 250 I=1,NPT
250 VLAG(I)=VLAG(I)+SUM*ZMAT(I,K)
BSUM=ZERO
DX=ZERO
DO 280 J=1,N
SUM=ZERO
DO 260 I=1,NPT
260 SUM=SUM+W(I)*BMAT(I,J)
BSUM=BSUM+SUM*D(J)
JP=NPT+J
DO 270 K=1,N
270 SUM=SUM+BMAT(JP,K)*D(K)
VLAG(JP)=SUM
BSUM=BSUM+SUM*D(J)
280 DX=DX+D(J)*XOPT(J)
BETA=DX*DX+DSQ*(XOPTSQ+DX+DX+HALF*DSQ)+BETA-BSUM
VLAG(KOPT)=VLAG(KOPT)+ONE
!
! If KNEW is positive and if the cancellation in DENOM is unacceptab
! then BIGDEN calculates an alternative model step, XNEW being used
! working space.
!
IF (KNEW .GT. 0) THEN
TEMP=ONE+ALPHA*BETA/VLAG(KNEW)**2
IF (DABS(TEMP) .LE. 0.8D0) THEN
CALL BIGDEN (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KOPT, &
& KNEW,D,W,VLAG,BETA,XNEW,W(NDIM+1),W(6*NDIM+1))
END IF
END IF
!
! Calculate the next value of the objective function.
!
290 DO 300 I=1,N
XNEW(I)=XOPT(I)+D(I)
300 X(I)=XBASE(I)+XNEW(I)
NF=NF+1
310 IF (NF .GT. NFTEST) THEN
NF=NF-1
IF (IPRINT .GT. 0) PRINT 320
320 FORMAT (/4X,'Return from NEWUOA because CALFUN has been', &
& ' called MAXFUN times.')
GOTO 530
END IF
CALL CALFUN (N,X,F)
IF (IPRINT .EQ. 3) THEN
PRINT 330, NF,F,(X(I),I=1,N)
330 FORMAT (/4X,'Function number',I6,' F =',1PD18.10, &
& ' The corresponding X is:'/(2X,5D15.6))
END IF
IF (NF .LE. NPT) GOTO 70
IF (KNEW .EQ. -1) GOTO 530
!
! Use the quadratic model to predict the change in F due to the step
! and set DIFF to the error of this prediction.
!
VQUAD=ZERO
IH=0
DO 340 J=1,N
VQUAD=VQUAD+D(J)*GQ(J)
DO 340 I=1,J
IH=IH+1
TEMP=D(I)*XNEW(J)+D(J)*XOPT(I)
IF (I .EQ. J) TEMP=HALF*TEMP
340 VQUAD=VQUAD+TEMP*HQ(IH)
DO 350 K=1,NPT
350 VQUAD=VQUAD+PQ(K)*W(K)
DIFF=F-FOPT-VQUAD
DIFFC=DIFFB
DIFFB=DIFFA
DIFFA=DABS(DIFF)
IF (DNORM .GT. RHO) NFSAV=NF
!
! Update FOPT and XOPT if the new F is the least value of the object
! function so far. The branch when KNEW is positive occurs if D is n
! a trust region step.
!
FSAVE=FOPT
IF (F .LT. FOPT) THEN
FOPT=F
XOPTSQ=ZERO
DO 360 I=1,N
XOPT(I)=XNEW(I)
360 XOPTSQ=XOPTSQ+XOPT(I)**2
END IF
KSAVE=KNEW
IF (KNEW .GT. 0) GOTO 410
!
! Pick the next value of DELTA after a trust region step.
!
IF (VQUAD .GE. ZERO) THEN
IF (IPRINT .GT. 0) PRINT 370
370 FORMAT (/4X,'Return from NEWUOA because a trust', &
& ' region step has failed to reduce Q.')
GOTO 530
END IF
RATIO=(F-FSAVE)/VQUAD
IF (RATIO .LE. TENTH) THEN
DELTA=HALF*DNORM
ELSE IF (RATIO .LE. 0.7D0) THEN
DELTA=DMAX1(HALF*DELTA,DNORM)
ELSE
DELTA=DMAX1(HALF*DELTA,DNORM+DNORM)
END IF
IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
!
! Set KNEW to the index of the next interpolation point to be delete
!
RHOSQ=DMAX1(TENTH*DELTA,RHO)**2
KTEMP=0
DETRAT=ZERO
IF (F .GE. FSAVE) THEN
KTEMP=KOPT
DETRAT=ONE
END IF
DO 400 K=1,NPT
HDIAG=ZERO
DO 380 J=1,NPTM
TEMP=ONE
IF (J .LT. IDZ) TEMP=-ONE
380 HDIAG=HDIAG+TEMP*ZMAT(K,J)**2
TEMP=DABS(BETA*HDIAG+VLAG(K)**2)
DISTSQ=ZERO
DO 390 J=1,N
390 DISTSQ=DISTSQ+(XPT(K,J)-XOPT(J))**2
IF (DISTSQ .GT. RHOSQ) TEMP=TEMP*(DISTSQ/RHOSQ)**3
IF (TEMP .GT. DETRAT .AND. K .NE. KTEMP) THEN
DETRAT=TEMP
KNEW=K
END IF
400 END DO
IF (KNEW .EQ. 0) GOTO 460
!
! Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point
! can be moved. Begin the updating of the quadratic model, starting
! with the explicit second derivative term.
!
410 CALL UPDATE (N,NPT,BMAT,ZMAT,IDZ,NDIM,VLAG,BETA,KNEW,W)
FVAL(KNEW)=F
IH=0
DO 420 I=1,N
TEMP=PQ(KNEW)*XPT(KNEW,I)
DO 420 J=1,I
IH=IH+1
420 HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J)
PQ(KNEW)=ZERO
!
! Update the other second derivative parameters, and then the gradie
! vector of the model. Also include the new interpolation point.
!
DO 440 J=1,NPTM
TEMP=DIFF*ZMAT(KNEW,J)
IF (J .LT. IDZ) TEMP=-TEMP
DO 440 K=1,NPT
440 PQ(K)=PQ(K)+TEMP*ZMAT(K,J)
GQSQ=ZERO
DO 450 I=1,N
GQ(I)=GQ(I)+DIFF*BMAT(KNEW,I)
GQSQ=GQSQ+GQ(I)**2
450 XPT(KNEW,I)=XNEW(I)
!
! If a trust region step makes a small change to the objective funct
! then calculate the gradient of the least Frobenius norm interpolan
! XBASE, and store it in W, using VLAG for a vector of right hand si
!
IF (KSAVE .EQ. 0 .AND. DELTA .EQ. RHO) THEN
IF (DABS(RATIO) .GT. 1.0D-2) THEN
ITEST=0
ELSE
DO 700 K=1,NPT
700 VLAG(K)=FVAL(K)-FVAL(KOPT)
GISQ=ZERO
DO 720 I=1,N
SUM=ZERO
DO 710 K=1,NPT
710 SUM=SUM+BMAT(K,I)*VLAG(K)
GISQ=GISQ+SUM*SUM
720 W(I)=SUM
!
! Test whether to replace the new quadratic model by the least Frobe
! norm interpolant, making the replacement if the test is satisfied.
!
ITEST=ITEST+1
IF (GQSQ .LT. 1.0D2*GISQ) ITEST=0
IF (ITEST .GE. 3) THEN
DO 730 I=1,N
730 GQ(I)=W(I)
DO 740 IH=1,NH
740 HQ(IH)=ZERO
DO 760 J=1,NPTM
W(J)=ZERO
DO 750 K=1,NPT
750 W(J)=W(J)+VLAG(K)*ZMAT(K,J)
760 IF (J .LT. IDZ) W(J)=-W(J)
DO 770 K=1,NPT
PQ(K)=ZERO
DO 770 J=1,NPTM
770 PQ(K)=PQ(K)+ZMAT(K,J)*W(J)
ITEST=0
END IF
END IF
END IF
IF (F .LT. FSAVE) KOPT=KNEW
!
! If a trust region step has provided a sufficient decrease in F, th
! branch for another trust region calculation. The case KSAVE>0 occu
! when the new function value was calculated by a model step.
!
IF (F .LE. FSAVE+TENTH*VQUAD) GOTO 100
IF (KSAVE .GT. 0) GOTO 100
!
! Alternatively, find out if the interpolation points are close enou
! to the best point so far.
!
KNEW=0
460 DISTSQ=4.0D0*DELTA*DELTA
DO 480 K=1,NPT
SUM=ZERO
DO 470 J=1,N
470 SUM=SUM+(XPT(K,J)-XOPT(J))**2
IF (SUM .GT. DISTSQ) THEN
KNEW=K
DISTSQ=SUM
END IF
480 END DO
!
! If KNEW is positive, then set DSTEP, and branch back for the next
! iteration, which will generate a "model step".
!
IF (KNEW .GT. 0) THEN
DSTEP=DMAX1(DMIN1(TENTH*DSQRT(DISTSQ),HALF*DELTA),RHO)
DSQ=DSTEP*DSTEP
GOTO 120
END IF
IF (RATIO .GT. ZERO) GOTO 100
IF (DMAX1(DELTA,DNORM) .GT. RHO) GOTO 100
!
! The calculations with the current value of RHO are complete. Pick
! next values of RHO and DELTA.
!
490 IF (RHO .GT. RHOEND) THEN
DELTA=HALF*RHO
RATIO=RHO/RHOEND
IF (RATIO .LE. 16.0D0) THEN
RHO=RHOEND
ELSE IF (RATIO .LE. 250.0D0) THEN
RHO=DSQRT(RATIO)*RHOEND
ELSE
RHO=TENTH*RHO
END IF
DELTA=DMAX1(DELTA,RHO)
IF (IPRINT .GE. 2) THEN
IF (IPRINT .GE. 3) PRINT 500
500 FORMAT (5X)
PRINT 510, RHO,NF
510 FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of', &
& ' function values =',I6)
PRINT 520, FOPT,(XBASE(I)+XOPT(I),I=1,N)
520 FORMAT (4X,'Least value of F =',1PD23.15,9X, &
& 'The corresponding X is:'/(2X,5D15.6))
END IF
GOTO 90
END IF
!
! Return from the calculation, after another Newton-Raphson step, if
! it is too short to have been tried before.
!
IF (KNEW .EQ. -1) GOTO 290
530 IF (FOPT .LE. F) THEN
DO 540 I=1,N
540 X(I)=XBASE(I)+XOPT(I)
F=FOPT
END IF
IF (IPRINT .GE. 1) THEN
PRINT 550, NF
550 FORMAT (/4X,'At the return from NEWUOA',5X, &
& 'Number of function values =',I6)
PRINT 520, F,(X(I),I=1,N)
END IF
RETURN
END
SUBROUTINE BIGDEN (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KOPT, &
& KNEW,D,W,VLAG,BETA,S,WVEC,PROD)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XOPT(*),XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),D(*), &
& W(*),VLAG(*),S(*),WVEC(NDIM,*),PROD(NDIM,*)
DIMENSION DEN(9),DENEX(9),PAR(9)
!
! N is the number of variables.
! NPT is the number of interpolation equations.
! XOPT is the best interpolation point so far.
! XPT contains the coordinates of the current interpolation points.
! BMAT provides the last N columns of H.
! ZMAT and IDZ give a factorization of the first NPT by NPT submatri
! NDIM is the first dimension of BMAT and has the value NPT+N.
! KOPT is the index of the optimal interpolation point.
! KNEW is the index of the interpolation point that is going to be m
! D will be set to the step from XOPT to the new point, and on entry
! should be the D that was calculated by the last call of BIGLAG.
! length of the initial D provides a trust region bound on the fin
! W will be set to Wcheck for the final choice of D.
! VLAG will be set to Theta*Wcheck+e_b for the final choice of D.
! BETA will be set to the value that will occur in the updating form
! when the KNEW-th interpolation point is moved to its new positio
! S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be us
! for working space.
!
! D is calculated in a way that should provide a denominator with a
! modulus in the updating formula when the KNEW-th interpolation poi
! shifted to the new position XOPT+D.
!
! Set some constants.
!
HALF=0.5D0
ONE=1.0D0
QUART=0.25D0
TWO=2.0D0
ZERO=0.0D0
TWOPI=8.0D0*DATAN(ONE)
NPTM=NPT-N-1
!
! Store the first NPT elements of the KNEW-th column of H in W(N+1)
! to W(N+NPT).
!
DO 10 K=1,NPT
10 W(N+K)=ZERO
DO 20 J=1,NPTM
TEMP=ZMAT(KNEW,J)
IF (J .LT. IDZ) TEMP=-TEMP
DO 20 K=1,NPT
20 W(N+K)=W(N+K)+TEMP*ZMAT(K,J)
ALPHA=W(N+KNEW)
!
! The initial search direction D is taken from the last call of BIGL
! and the initial S is set below, usually to the direction from X_OP
! to X_KNEW, but a different direction to an interpolation point may
! be chosen, in order to prevent S from being nearly parallel to D.
!
DD=ZERO
DS=ZERO
SS=ZERO
XOPTSQ=ZERO
DO 30 I=1,N
DD=DD+D(I)**2
S(I)=XPT(KNEW,I)-XOPT(I)
DS=DS+D(I)*S(I)
SS=SS+S(I)**2
30 XOPTSQ=XOPTSQ+XOPT(I)**2
IF (DS*DS .GT. 0.99D0*DD*SS) THEN
KSAV=KNEW
DTEST=DS*DS/SS
DO 50 K=1,NPT
IF (K .NE. KOPT) THEN
DSTEMP=ZERO
SSTEMP=ZERO
DO 40 I=1,N
DIFF=XPT(K,I)-XOPT(I)
DSTEMP=DSTEMP+D(I)*DIFF
40 SSTEMP=SSTEMP+DIFF*DIFF
IF (DSTEMP*DSTEMP/SSTEMP .LT. DTEST) THEN
KSAV=K
DTEST=DSTEMP*DSTEMP/SSTEMP
DS=DSTEMP
SS=SSTEMP
END IF
END IF
50 CONTINUE
DO 60 I=1,N
60 S(I)=XPT(KSAV,I)-XOPT(I)
END IF
SSDEN=DD*SS-DS*DS
ITERC=0
DENSAV=ZERO
!
! Begin the iteration by overwriting S with a vector that has the
! required length and direction.
!
70 ITERC=ITERC+1
TEMP=ONE/DSQRT(SSDEN)
XOPTD=ZERO
XOPTS=ZERO
DO 80 I=1,N
S(I)=TEMP*(DD*S(I)-DS*D(I))
XOPTD=XOPTD+XOPT(I)*D(I)
80 XOPTS=XOPTS+XOPT(I)*S(I)
!
! Set the coefficients of the first two terms of BETA.
!
TEMPA=HALF*XOPTD*XOPTD
TEMPB=HALF*XOPTS*XOPTS
DEN(1)=DD*(XOPTSQ+HALF*DD)+TEMPA+TEMPB
DEN(2)=TWO*XOPTD*DD
DEN(3)=TWO*XOPTS*DD
DEN(4)=TEMPA-TEMPB
DEN(5)=XOPTD*XOPTS
DO 90 I=6,9
90 DEN(I)=ZERO
!
! Put the coefficients of Wcheck in WVEC.
!
DO 110 K=1,NPT
TEMPA=ZERO
TEMPB=ZERO
TEMPC=ZERO
DO 100 I=1,N
TEMPA=TEMPA+XPT(K,I)*D(I)
TEMPB=TEMPB+XPT(K,I)*S(I)
100 TEMPC=TEMPC+XPT(K,I)*XOPT(I)
WVEC(K,1)=QUART*(TEMPA*TEMPA+TEMPB*TEMPB)
WVEC(K,2)=TEMPA*TEMPC
WVEC(K,3)=TEMPB*TEMPC
WVEC(K,4)=QUART*(TEMPA*TEMPA-TEMPB*TEMPB)
110 WVEC(K,5)=HALF*TEMPA*TEMPB
DO 120 I=1,N
IP=I+NPT
WVEC(IP,1)=ZERO
WVEC(IP,2)=D(I)
WVEC(IP,3)=S(I)
WVEC(IP,4)=ZERO
120 WVEC(IP,5)=ZERO
!
! Put the coefficients of THETA*Wcheck in PROD.
!
DO 190 JC=1,5
NW=NPT
IF (JC .EQ. 2 .OR. JC .EQ. 3) NW=NDIM
DO 130 K=1,NPT
130 PROD(K,JC)=ZERO
DO 150 J=1,NPTM
SUM=ZERO
DO 140 K=1,NPT
140 SUM=SUM+ZMAT(K,J)*WVEC(K,JC)
IF (J .LT. IDZ) SUM=-SUM
DO 150 K=1,NPT
150 PROD(K,JC)=PROD(K,JC)+SUM*ZMAT(K,J)
IF (NW .EQ. NDIM) THEN
DO 170 K=1,NPT
SUM=ZERO
DO 160 J=1,N
160 SUM=SUM+BMAT(K,J)*WVEC(NPT+J,JC)
170 PROD(K,JC)=PROD(K,JC)+SUM
END IF
DO 190 J=1,N
SUM=ZERO
DO 180 I=1,NW
180 SUM=SUM+BMAT(I,J)*WVEC(I,JC)
190 PROD(NPT+J,JC)=SUM
!
! Include in DEN the part of BETA that depends on THETA.
!
DO 210 K=1,NDIM
SUM=ZERO
DO 200 I=1,5
PAR(I)=HALF*PROD(K,I)*WVEC(K,I)
200 SUM=SUM+PAR(I)
DEN(1)=DEN(1)-PAR(1)-SUM
TEMPA=PROD(K,1)*WVEC(K,2)+PROD(K,2)*WVEC(K,1)
TEMPB=PROD(K,2)*WVEC(K,4)+PROD(K,4)*WVEC(K,2)
TEMPC=PROD(K,3)*WVEC(K,5)+PROD(K,5)*WVEC(K,3)
DEN(2)=DEN(2)-TEMPA-HALF*(TEMPB+TEMPC)
DEN(6)=DEN(6)-HALF*(TEMPB-TEMPC)
TEMPA=PROD(K,1)*WVEC(K,3)+PROD(K,3)*WVEC(K,1)
TEMPB=PROD(K,2)*WVEC(K,5)+PROD(K,5)*WVEC(K,2)
TEMPC=PROD(K,3)*WVEC(K,4)+PROD(K,4)*WVEC(K,3)
DEN(3)=DEN(3)-TEMPA-HALF*(TEMPB-TEMPC)
DEN(7)=DEN(7)-HALF*(TEMPB+TEMPC)
TEMPA=PROD(K,1)*WVEC(K,4)+PROD(K,4)*WVEC(K,1)
DEN(4)=DEN(4)-TEMPA-PAR(2)+PAR(3)
TEMPA=PROD(K,1)*WVEC(K,5)+PROD(K,5)*WVEC(K,1)
TEMPB=PROD(K,2)*WVEC(K,3)+PROD(K,3)*WVEC(K,2)
DEN(5)=DEN(5)-TEMPA-HALF*TEMPB
DEN(8)=DEN(8)-PAR(4)+PAR(5)
TEMPA=PROD(K,4)*WVEC(K,5)+PROD(K,5)*WVEC(K,4)
210 DEN(9)=DEN(9)-HALF*TEMPA
!
! Extend DEN so that it holds all the coefficients of DENOM.
!
SUM=ZERO
DO 220 I=1,5
PAR(I)=HALF*PROD(KNEW,I)**2
220 SUM=SUM+PAR(I)
DENEX(1)=ALPHA*DEN(1)+PAR(1)+SUM
TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,2)
TEMPB=PROD(KNEW,2)*PROD(KNEW,4)
TEMPC=PROD(KNEW,3)*PROD(KNEW,5)
DENEX(2)=ALPHA*DEN(2)+TEMPA+TEMPB+TEMPC
DENEX(6)=ALPHA*DEN(6)+TEMPB-TEMPC
TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,3)
TEMPB=PROD(KNEW,2)*PROD(KNEW,5)
TEMPC=PROD(KNEW,3)*PROD(KNEW,4)
DENEX(3)=ALPHA*DEN(3)+TEMPA+TEMPB-TEMPC
DENEX(7)=ALPHA*DEN(7)+TEMPB+TEMPC
TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,4)
DENEX(4)=ALPHA*DEN(4)+TEMPA+PAR(2)-PAR(3)
TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,5)
DENEX(5)=ALPHA*DEN(5)+TEMPA+PROD(KNEW,2)*PROD(KNEW,3)
DENEX(8)=ALPHA*DEN(8)+PAR(4)-PAR(5)
DENEX(9)=ALPHA*DEN(9)+PROD(KNEW,4)*PROD(KNEW,5)
!
! Seek the value of the angle that maximizes the modulus of DENOM.
!
SUM=DENEX(1)+DENEX(2)+DENEX(4)+DENEX(6)+DENEX(8)
DENOLD=SUM
DENMAX=SUM
ISAVE=0
IU=49
TEMP=TWOPI/DBLE(IU+1)
PAR(1)=ONE
DO 250 I=1,IU
ANGLE=DBLE(I)*TEMP
PAR(2)=DCOS(ANGLE)
PAR(3)=DSIN(ANGLE)
DO 230 J=4,8,2
PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1)
230 PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2)
SUMOLD=SUM
SUM=ZERO
DO 240 J=1,9
240 SUM=SUM+DENEX(J)*PAR(J)
IF (DABS(SUM) .GT. DABS(DENMAX)) THEN
DENMAX=SUM
ISAVE=I
TEMPA=SUMOLD
ELSE IF (I .EQ. ISAVE+1) THEN
TEMPB=SUM
END IF
250 END DO
IF (ISAVE .EQ. 0) TEMPA=SUM
IF (ISAVE .EQ. IU) TEMPB=DENOLD
STEP=ZERO
IF (TEMPA .NE. TEMPB) THEN
TEMPA=TEMPA-DENMAX
TEMPB=TEMPB-DENMAX
STEP=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
END IF
ANGLE=TEMP*(DBLE(ISAVE)+STEP)
!
! Calculate the new parameters of the denominator, the new VLAG vect
! and the new D. Then test for convergence.
!
PAR(2)=DCOS(ANGLE)
PAR(3)=DSIN(ANGLE)
DO 260 J=4,8,2
PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1)
260 PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2)
BETA=ZERO
DENMAX=ZERO
DO 270 J=1,9
BETA=BETA+DEN(J)*PAR(J)
270 DENMAX=DENMAX+DENEX(J)*PAR(J)
DO 280 K=1,NDIM
VLAG(K)=ZERO
DO 280 J=1,5
280 VLAG(K)=VLAG(K)+PROD(K,J)*PAR(J)
TAU=VLAG(KNEW)
DD=ZERO
TEMPA=ZERO
TEMPB=ZERO
DO 290 I=1,N
D(I)=PAR(2)*D(I)+PAR(3)*S(I)
W(I)=XOPT(I)+D(I)
DD=DD+D(I)**2
TEMPA=TEMPA+D(I)*W(I)
290 TEMPB=TEMPB+W(I)*W(I)
IF (ITERC .GE. N) GOTO 340
IF (ITERC .GT. 1) DENSAV=DMAX1(DENSAV,DENOLD)
IF (DABS(DENMAX) .LE. 1.1D0*DABS(DENSAV)) GOTO 340
DENSAV=DENMAX
!
! Set S to half the gradient of the denominator with respect to D.
! Then branch for the next iteration.
!
DO 300 I=1,N
TEMP=TEMPA*XOPT(I)+TEMPB*D(I)-VLAG(NPT+I)
300 S(I)=TAU*BMAT(KNEW,I)+ALPHA*TEMP
DO 320 K=1,NPT
SUM=ZERO
DO 310 J=1,N
310 SUM=SUM+XPT(K,J)*W(J)
TEMP=(TAU*W(N+K)-ALPHA*VLAG(K))*SUM
DO 320 I=1,N
320 S(I)=S(I)+TEMP*XPT(K,I)
SS=ZERO
DS=ZERO
DO 330 I=1,N
SS=SS+S(I)**2
330 DS=DS+D(I)*S(I)
SSDEN=DD*SS-DS*DS
IF (SSDEN .GE. 1.0D-8*DD*SS) GOTO 70
!
! Set the vector W before the RETURN from the subroutine.
!
340 DO 350 K=1,NDIM
W(K)=ZERO
DO 350 J=1,5
350 W(K)=W(K)+WVEC(K,J)*PAR(J)
VLAG(KOPT)=VLAG(KOPT)+ONE
RETURN
END
SUBROUTINE BIGLAG (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KNEW, &
& DELTA,D,ALPHA,HCOL,GC,GD,S,W)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XOPT(*),XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),D(*), &
& HCOL(*),GC(*),GD(*),S(*),W(*)
!
! N is the number of variables.
! NPT is the number of interpolation equations.
! XOPT is the best interpolation point so far.
! XPT contains the coordinates of the current interpolation points.
! BMAT provides the last N columns of H.
! ZMAT and IDZ give a factorization of the first NPT by NPT submatri
! NDIM is the first dimension of BMAT and has the value NPT+N.
! KNEW is the index of the interpolation point that is going to be m
! DELTA is the current trust region bound.
! D will be set to the step from XOPT to the new point.
! ALPHA will be set to the KNEW-th diagonal element of the H matrix.
! HCOL, GC, GD, S and W will be used for working space.
!
! The step D is calculated in a way that attempts to maximize the mo
! of LFUNC(XOPT+D), subject to the bound ||D|| .LE. DELTA, where LFU
! the KNEW-th Lagrange function.
!
! Set some constants.
!
HALF=0.5D0
ONE=1.0D0
ZERO=0.0D0
TWOPI=8.0D0*DATAN(ONE)
DELSQ=DELTA*DELTA
NPTM=NPT-N-1
!
! Set the first NPT components of HCOL to the leading elements of th
! KNEW-th column of H.
!
ITERC=0
DO 10 K=1,NPT
10 HCOL(K)=ZERO
DO 20 J=1,NPTM
TEMP=ZMAT(KNEW,J)
IF (J .LT. IDZ) TEMP=-TEMP
DO 20 K=1,NPT
20 HCOL(K)=HCOL(K)+TEMP*ZMAT(K,J)
ALPHA=HCOL(KNEW)
!
! Set the unscaled initial direction D. Form the gradient of LFUNC a
! XOPT, and multiply D by the second derivative matrix of LFUNC.
!
DD=ZERO
DO 30 I=1,N
D(I)=XPT(KNEW,I)-XOPT(I)
GC(I)=BMAT(KNEW,I)
GD(I)=ZERO
30 DD=DD+D(I)**2
DO 50 K=1,NPT
TEMP=ZERO
SUM=ZERO
DO 40 J=1,N
TEMP=TEMP+XPT(K,J)*XOPT(J)
40 SUM=SUM+XPT(K,J)*D(J)
TEMP=HCOL(K)*TEMP
SUM=HCOL(K)*SUM
DO 50 I=1,N
GC(I)=GC(I)+TEMP*XPT(K,I)
50 GD(I)=GD(I)+SUM*XPT(K,I)
!
! Scale D and GD, with a sign change if required. Set S to another
! vector in the initial two dimensional subspace.
!
GG=ZERO
SP=ZERO
DHD=ZERO
DO 60 I=1,N
GG=GG+GC(I)**2
SP=SP+D(I)*GC(I)
60 DHD=DHD+D(I)*GD(I)
SCALE=DELTA/DSQRT(DD)
IF (SP*DHD .LT. ZERO) SCALE=-SCALE
TEMP=ZERO
IF (SP*SP .GT. 0.99D0*DD*GG) TEMP=ONE
TAU=SCALE*(DABS(SP)+HALF*SCALE*DABS(DHD))
IF (GG*DELSQ .LT. 0.01D0*TAU*TAU) TEMP=ONE
DO 70 I=1,N
D(I)=SCALE*D(I)
GD(I)=SCALE*GD(I)
70 S(I)=GC(I)+TEMP*GD(I)
!
! Begin the iteration by overwriting S with a vector that has the
! required length and direction, except that termination occurs if
! the given D and S are nearly parallel.
!
80 ITERC=ITERC+1
DD=ZERO
SP=ZERO
SS=ZERO
DO 90 I=1,N
DD=DD+D(I)**2
SP=SP+D(I)*S(I)
90 SS=SS+S(I)**2
TEMP=DD*SS-SP*SP
IF (TEMP .LE. 1.0D-8*DD*SS) GOTO 160
DENOM=DSQRT(TEMP)
DO 100 I=1,N
S(I)=(DD*S(I)-SP*D(I))/DENOM
100 W(I)=ZERO
!
! Calculate the coefficients of the objective function on the circle
! beginning with the multiplication of S by the second derivative ma
!
DO 120 K=1,NPT
SUM=ZERO
DO 110 J=1,N
110 SUM=SUM+XPT(K,J)*S(J)
SUM=HCOL(K)*SUM
DO 120 I=1,N
120 W(I)=W(I)+SUM*XPT(K,I)
CF1=ZERO
CF2=ZERO
CF3=ZERO
CF4=ZERO
CF5=ZERO
DO 130 I=1,N
CF1=CF1+S(I)*W(I)
CF2=CF2+D(I)*GC(I)
CF3=CF3+S(I)*GC(I)
CF4=CF4+D(I)*GD(I)
130 CF5=CF5+S(I)*GD(I)
CF1=HALF*CF1
CF4=HALF*CF4-CF1
!
! Seek the value of the angle that maximizes the modulus of TAU.
!
TAUBEG=CF1+CF2+CF4
TAUMAX=TAUBEG
TAUOLD=TAUBEG
ISAVE=0
IU=49
TEMP=TWOPI/DBLE(IU+1)
DO 140 I=1,IU
ANGLE=DBLE(I)*TEMP
CTH=DCOS(ANGLE)
STH=DSIN(ANGLE)
TAU=CF1+(CF2+CF4*CTH)*CTH+(CF3+CF5*CTH)*STH
IF (DABS(TAU) .GT. DABS(TAUMAX)) THEN
TAUMAX=TAU
ISAVE=I
TEMPA=TAUOLD
ELSE IF (I .EQ. ISAVE+1) THEN
TEMPB=TAU
END IF
140 TAUOLD=TAU
IF (ISAVE .EQ. 0) TEMPA=TAU
IF (ISAVE .EQ. IU) TEMPB=TAUBEG
STEP=ZERO
IF (TEMPA .NE. TEMPB) THEN
TEMPA=TEMPA-TAUMAX
TEMPB=TEMPB-TAUMAX
STEP=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
END IF
ANGLE=TEMP*(DBLE(ISAVE)+STEP)
!
! Calculate the new D and GD. Then test for convergence.
!
CTH=DCOS(ANGLE)
STH=DSIN(ANGLE)
TAU=CF1+(CF2+CF4*CTH)*CTH+(CF3+CF5*CTH)*STH
DO 150 I=1,N
D(I)=CTH*D(I)+STH*S(I)
GD(I)=CTH*GD(I)+STH*W(I)
150 S(I)=GC(I)+GD(I)
IF (DABS(TAU) .LE. 1.1D0*DABS(TAUBEG)) GOTO 160
IF (ITERC .LT. N) GOTO 80
160 RETURN
END
SUBROUTINE TRSAPP (N,NPT,XOPT,XPT,GQ,HQ,PQ,DELTA,STEP, &
& D,G,HD,HS,CRVMIN)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XOPT(*),XPT(NPT,*),GQ(*),HQ(*),PQ(*),STEP(*), &
& D(*),G(*),HD(*),HS(*)
!
! N is the number of variables of a quadratic objective function, Q
! The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meani
! in order to define the current quadratic model Q.
! DELTA is the trust region radius, and has to be positive.
! STEP will be set to the calculated trial step.
! The arrays D, G, HD and HS will be used for working space.
! CRVMIN will be set to the least curvature of H along the conjugate
! directions that occur, except that it is set to zero if STEP goe
! all the way to the trust region boundary.
!
! The calculation of STEP begins with the truncated conjugate gradie
! method. If the boundary of the trust region is reached, then furth
! changes to STEP may be made, each one being in the 2D space spanne
! by the current STEP and the corresponding gradient of Q. Thus STEP
! should provide a substantial reduction to Q within the trust regio
!
! Initialization, which includes setting HD to H times XOPT.
!
HALF=0.5D0
ZERO=0.0D0
TWOPI=8.0D0*DATAN(1.0D0)
DELSQ=DELTA*DELTA
ITERC=0
ITERMAX=N
ITERSW=ITERMAX
DO 10 I=1,N
10 D(I)=XOPT(I)
GOTO 170
!
! Prepare for the first line search.
!
20 QRED=ZERO
DD=ZERO
DO 30 I=1,N
STEP(I)=ZERO
HS(I)=ZERO
G(I)=GQ(I)+HD(I)
D(I)=-G(I)
30 DD=DD+D(I)**2
CRVMIN=ZERO
IF (DD .EQ. ZERO) GOTO 160
DS=ZERO
SS=ZERO
GG=DD
GGBEG=GG
!
! Calculate the step to the trust region boundary and the product HD
!
40 ITERC=ITERC+1
TEMP=DELSQ-SS
BSTEP=TEMP/(DS+DSQRT(DS*DS+DD*TEMP))
GOTO 170
50 DHD=ZERO
DO 60 J=1,N
60 DHD=DHD+D(J)*HD(J)
!
! Update CRVMIN and set the step-length ALPHA.
!
ALPHA=BSTEP
IF (DHD .GT. ZERO) THEN
TEMP=DHD/DD
IF (ITERC .EQ. 1) CRVMIN=TEMP
CRVMIN=DMIN1(CRVMIN,TEMP)
ALPHA=DMIN1(ALPHA,GG/DHD)
END IF
QADD=ALPHA*(GG-HALF*ALPHA*DHD)
QRED=QRED+QADD
!
! Update STEP and HS.
!
GGSAV=GG
GG=ZERO
DO 70 I=1,N
STEP(I)=STEP(I)+ALPHA*D(I)
HS(I)=HS(I)+ALPHA*HD(I)
70 GG=GG+(G(I)+HS(I))**2
!
! Begin another conjugate direction iteration if required.
!
IF (ALPHA .LT. BSTEP) THEN
IF (QADD .LE. 0.01D0*QRED) GOTO 160
IF (GG .LE. 1.0D-4*GGBEG) GOTO 160
IF (ITERC .EQ. ITERMAX) GOTO 160
TEMP=GG/GGSAV
DD=ZERO
DS=ZERO
SS=ZERO
DO 80 I=1,N
D(I)=TEMP*D(I)-G(I)-HS(I)
DD=DD+D(I)**2
DS=DS+D(I)*STEP(I)
80 SS=SS+STEP(I)**2
IF (DS .LE. ZERO) GOTO 160
IF (SS .LT. DELSQ) GOTO 40
END IF
CRVMIN=ZERO
ITERSW=ITERC
!
! Test whether an alternative iteration is required.
!
90 IF (GG .LE. 1.0D-4*GGBEG) GOTO 160
SG=ZERO
SHS=ZERO
DO 100 I=1,N
SG=SG+STEP(I)*G(I)
100 SHS=SHS+STEP(I)*HS(I)
SGK=SG+SHS
ANGTEST=SGK/DSQRT(GG*DELSQ)
IF (ANGTEST .LE. -0.99D0) GOTO 160
!
! Begin the alternative iteration by calculating D and HD and some
! scalar products.
!
ITERC=ITERC+1
TEMP=DSQRT(DELSQ*GG-SGK*SGK)
TEMPA=DELSQ/TEMP
TEMPB=SGK/TEMP
DO 110 I=1,N
110 D(I)=TEMPA*(G(I)+HS(I))-TEMPB*STEP(I)
GOTO 170
120 DG=ZERO
DHD=ZERO
DHS=ZERO
DO 130 I=1,N
DG=DG+D(I)*G(I)
DHD=DHD+HD(I)*D(I)
130 DHS=DHS+HD(I)*STEP(I)
!
! Seek the value of the angle that minimizes Q.
!
CF=HALF*(SHS-DHD)
QBEG=SG+CF
QSAV=QBEG
QMIN=QBEG
ISAVE=0
IU=49
TEMP=TWOPI/DBLE(IU+1)
DO 140 I=1,IU
ANGLE=DBLE(I)*TEMP
CTH=DCOS(ANGLE)
STH=DSIN(ANGLE)
QNEW=(SG+CF*CTH)*CTH+(DG+DHS*CTH)*STH
IF (QNEW .LT. QMIN) THEN
QMIN=QNEW
ISAVE=I
TEMPA=QSAV
ELSE IF (I .EQ. ISAVE+1) THEN
TEMPB=QNEW
END IF
140 QSAV=QNEW
IF (ISAVE .EQ. ZERO) TEMPA=QNEW
IF (ISAVE .EQ. IU) TEMPB=QBEG
ANGLE=ZERO
IF (TEMPA .NE. TEMPB) THEN
TEMPA=TEMPA-QMIN
TEMPB=TEMPB-QMIN
ANGLE=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
END IF
ANGLE=TEMP*(DBLE(ISAVE)+ANGLE)
!
! Calculate the new STEP and HS. Then test for convergence.
!
CTH=DCOS(ANGLE)
STH=DSIN(ANGLE)
REDUC=QBEG-(SG+CF*CTH)*CTH-(DG+DHS*CTH)*STH
GG=ZERO
DO 150 I=1,N
STEP(I)=CTH*STEP(I)+STH*D(I)
HS(I)=CTH*HS(I)+STH*HD(I)
150 GG=GG+(G(I)+HS(I))**2
QRED=QRED+REDUC
RATIO=REDUC/QRED
IF (ITERC .LT. ITERMAX .AND. RATIO .GT. 0.01D0) GOTO 90
160 RETURN
!
! The following instructions act as a subroutine for setting the vec
! HD to the vector D multiplied by the second derivative matrix of Q
! They are called from three different places, which are distinguish
! by the value of ITERC.
!
170 DO 180 I=1,N
180 HD(I)=ZERO
DO 200 K=1,NPT
TEMP=ZERO
DO 190 J=1,N
190 TEMP=TEMP+XPT(K,J)*D(J)
TEMP=TEMP*PQ(K)
DO 200 I=1,N
200 HD(I)=HD(I)+TEMP*XPT(K,I)
IH=0
DO 210 J=1,N
DO 210 I=1,J
IH=IH+1
IF (I .LT. J) HD(J)=HD(J)+HQ(IH)*D(I)
210 HD(I)=HD(I)+HQ(IH)*D(J)
IF (ITERC .EQ. 0) GOTO 20
IF (ITERC .LE. ITERSW) GOTO 50
GOTO 120
END
SUBROUTINE UPDATE (N,NPT,BMAT,ZMAT,IDZ,NDIM,VLAG,BETA,KNEW,W)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION BMAT(NDIM,*),ZMAT(NPT,*),VLAG(*),W(*)
!
! The arrays BMAT and ZMAT with IDZ are updated, in order to shift t
! interpolation point that has index KNEW. On entry, VLAG contains t
! components of the vector Theta*Wcheck+e_b of the updating formula
! (6.11), and BETA holds the value of the parameter that has this na
! The vector W is used for working space.
!
! Set some constants.
!
ONE=1.0D0
ZERO=0.0D0
NPTM=NPT-N-1
!
! Apply the rotations that put zeros in the KNEW-th row of ZMAT.
!
JL=1
DO 20 J=2,NPTM
IF (J .EQ. IDZ) THEN
JL=IDZ
ELSE IF (ZMAT(KNEW,J) .NE. ZERO) THEN
TEMP=DSQRT(ZMAT(KNEW,JL)**2+ZMAT(KNEW,J)**2)
TEMPA=ZMAT(KNEW,JL)/TEMP
TEMPB=ZMAT(KNEW,J)/TEMP
DO 10 I=1,NPT
TEMP=TEMPA*ZMAT(I,JL)+TEMPB*ZMAT(I,J)
ZMAT(I,J)=TEMPA*ZMAT(I,J)-TEMPB*ZMAT(I,JL)
10 ZMAT(I,JL)=TEMP
ZMAT(KNEW,J)=ZERO
END IF
20 END DO
!
! Put the first NPT components of the KNEW-th column of HLAG into W,
! and calculate the parameters of the updating formula.
!
TEMPA=ZMAT(KNEW,1)
IF (IDZ .GE. 2) TEMPA=-TEMPA
IF (JL .GT. 1) TEMPB=ZMAT(KNEW,JL)
DO 30 I=1,NPT
W(I)=TEMPA*ZMAT(I,1)
IF (JL .GT. 1) W(I)=W(I)+TEMPB*ZMAT(I,JL)
30 END DO
ALPHA=W(KNEW)
TAU=VLAG(KNEW)
TAUSQ=TAU*TAU
DENOM=ALPHA*BETA+TAUSQ
VLAG(KNEW)=VLAG(KNEW)-ONE
!
! Complete the updating of ZMAT when there is only one nonzero eleme
! in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to
! then the first column of ZMAT will be exchanged with another one l
!
IFLAG=0
IF (JL .EQ. 1) THEN
TEMP=DSQRT(DABS(DENOM))
TEMPB=TEMPA/TEMP
TEMPA=TAU/TEMP
DO 40 I=1,NPT
40 ZMAT(I,1)=TEMPA*ZMAT(I,1)-TEMPB*VLAG(I)
IF (IDZ .EQ. 1 .AND. TEMP .LT. ZERO) IDZ=2
IF (IDZ .GE. 2 .AND. TEMP .GE. ZERO) IFLAG=1
ELSE
!
! Complete the updating of ZMAT in the alternative case.
!
JA=1
IF (BETA .GE. ZERO) JA=JL
JB=JL+1-JA
TEMP=ZMAT(KNEW,JB)/DENOM
TEMPA=TEMP*BETA
TEMPB=TEMP*TAU
TEMP=ZMAT(KNEW,JA)
SCALA=ONE/DSQRT(DABS(BETA)*TEMP*TEMP+TAUSQ)
SCALB=SCALA*DSQRT(DABS(DENOM))
DO 50 I=1,NPT
ZMAT(I,JA)=SCALA*(TAU*ZMAT(I,JA)-TEMP*VLAG(I))
50 ZMAT(I,JB)=SCALB*(ZMAT(I,JB)-TEMPA*W(I)-TEMPB*VLAG(I))
IF (DENOM .LE. ZERO) THEN
IF (BETA .LT. ZERO) IDZ=IDZ+1
IF (BETA .GE. ZERO) IFLAG=1
END IF
END IF
!
! IDZ is reduced in the following case, and usually the first column
! of ZMAT is exchanged with a later one.
!
IF (IFLAG .EQ. 1) THEN
IDZ=IDZ-1
DO 60 I=1,NPT
TEMP=ZMAT(I,1)
ZMAT(I,1)=ZMAT(I,IDZ)
60 ZMAT(I,IDZ)=TEMP
END IF
!
! Finally, update the matrix BMAT.
!
DO 70 J=1,N
JP=NPT+J
W(JP)=BMAT(KNEW,J)
TEMPA=(ALPHA*VLAG(JP)-TAU*W(JP))/DENOM
TEMPB=(-BETA*W(JP)-TAU*VLAG(JP))/DENOM
DO 70 I=1,JP
BMAT(I,J)=BMAT(I,J)+TEMPA*VLAG(I)+TEMPB*W(I)
IF (I .GT. NPT) BMAT(JP,I-NPT)=BMAT(I,J)
70 CONTINUE
RETURN
END
Reference in a new issue View git blame Copy permalink