!-----------------------------------------------------------------------
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