MODULE cobyla2

! The Fortran 77 version of this code was by Michael Powell
! (M.J.D.Powell @ damtp.cam.ac.uk)

! This Fortran 90 version by Alan Miller
! This version is in a subset of F90 and is compatible with version 3.0a
! of Lahey's free ELF90 compiler.
! Alan.Miller @ vic.cmis.csiro.au
! Latest revision - 28 June 2000

IMPLICIT NONE

INTEGER, PARAMETER, PUBLIC :: dp = SELECTED_REAL_KIND(14, 60)

CONTAINS

!------------------------------------------------------------------------

SUBROUTINE cobyla (n, m, x, rhobeg, rhoend, iprint, maxfun)

INTEGER, INTENT(IN)        :: n
INTEGER, INTENT(IN)        :: m
REAL (dp), INTENT(IN OUT)  :: x(:)
REAL (dp), INTENT(IN OUT)  :: rhobeg
REAL (dp), INTENT(IN OUT)  :: rhoend
INTEGER, INTENT(IN)        :: iprint
INTEGER, INTENT(IN OUT)    :: maxfun

!  This subroutine minimizes an objective function F(X) subject to M
!  inequality constraints on X, where X is a vector of variables that has
!  N components.  The algorithm employs linear approximations to the
!  objective and constraint functions, the approximations being formed by
!  linear interpolation at N+1 points in the space of the variables.
!  We regard these interpolation points as vertices of a simplex.  The
!  parameter RHO controls the size of the simplex and it is reduced
!  automatically from RHOBEG to RHOEND.  For each RHO the subroutine tries
!  to achieve a good vector of variables for the current size, and then
!  RHO is reduced until the value RHOEND is reached.  Therefore RHOBEG and
!  RHOEND should be set to reasonable initial changes to and the required
!  accuracy in the variables respectively, but this accuracy should be
!  viewed as a subject for experimentation because it is not guaranteed.
!  The subroutine has an advantage over many of its competitors, however,
!  which is that it treats each constraint individually when calculating
!  a change to the variables, instead of lumping the constraints together
!  into a single penalty function.  The name of the subroutine is derived
!  from the phrase Constrained Optimization BY Linear Approximations.

!  The user must set the values of N, M, RHOBEG and RHOEND, and must
!  provide an initial vector of variables in X.  Further, the value of
!  IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
!  printing during the calculation. Specifically, there is no output if
!  IPRINT=0 and there is output only at the end of the calculation if
!  IPRINT=1.  Otherwise each new value of RHO and SIGMA is printed.
!  Further, the vector of variables and some function information are
!  given either when RHO is reduced or when each new value of F(X) is
!  computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
!  is a penalty parameter, it being assumed that a change to X is an
!  improvement if it reduces the merit function
!             F(X)+SIGMA*MAX(0.0, - C1(X), - C2(X),..., - CM(X)),
!  where C1,C2,...,CM denote the constraint functions that should become
!  nonnegative eventually, at least to the precision of RHOEND. In the
!  printed output the displayed term that is multiplied by SIGMA is
!  called MAXCV, which stands for 'MAXimum Constraint Violation'.  The
!  argument MAXFUN is an integer variable that must be set by the user to a
!  limit on the number of calls of CALCFC, the purpose of this routine being
!  given below.  The value of MAXFUN will be altered to the number of calls
!  of CALCFC that are made.

!  In order to define the objective and constraint functions, we require
!  a subroutine that has the name and arguments
!             SUBROUTINE CALCFC (N,M,X,F,CON)
!             DIMENSION X(:),CON(:)  .
!  The values of N and M are fixed and have been defined already, while
!  X is now the current vector of variables. The subroutine should return
!  the objective and constraint functions at X in F and CON(1),CON(2),
!  ...,CON(M).  Note that we are trying to adjust X so that F(X) is as
!  small as possible subject to the constraint functions being nonnegative.

!  N.B. If the starting value for any x(i) is set to zero, that value will
!       not be changed.   This can be a useful feature in comparing
!       nested models.   If all the x(i)'s are set to zero, an error message
!       will result.

! Local variable

INTEGER :: mpp

mpp = m + 2
CALL cobylb (n, m, mpp, x, rhobeg, rhoend, iprint, maxfun)

RETURN
END SUBROUTINE cobyla


!------------------------------------------------------------------------------

SUBROUTINE cobylb (n, m, mpp, x, rhobeg, rhoend, iprint, maxfun)

! N.B. Arguments CON, SIM, SIMI, DATMAT, A, VSIG, VETA, SIGBAR, DX, W & IACT
!   have been removed.

INTEGER, INTENT(IN)        :: n
INTEGER, INTENT(IN)        :: m
INTEGER, INTENT(IN)        :: mpp
REAL (dp), INTENT(IN OUT)  :: x(:)
REAL (dp), INTENT(IN)      :: rhobeg
REAL (dp), INTENT(IN)      :: rhoend
INTEGER, INTENT(IN)        :: iprint
INTEGER, INTENT(OUT)       :: maxfun

INTERFACE
  SUBROUTINE calcfc (n, m, x, f, con)
    IMPLICIT NONE
    INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)
    INTEGER, INTENT(IN)    :: n, m
    REAL (dp), INTENT(IN)  :: x(:)
    REAL (dp), INTENT(OUT) :: f
    REAL (dp), INTENT(OUT) :: con(:)
  END SUBROUTINE calcfc
END INTERFACE

!  Set the initial values of some parameters. The last column of SIM holds
!  the optimal vertex of the current simplex, and the preceding N columns
!  hold the displacements from the optimal vertex to the other vertices.
!  Further, SIMI holds the inverse of the matrix that is contained in the
!  first N columns of SIM.

! Local variables

REAL (dp) :: con(mpp), sim(n,n+1), simi(n,n), datmat(mpp,n+1), a(n,m+1),      &
             vsig(n), veta(n), sigbar(n), dx(n), w(n)
REAL (dp) :: alpha, barmu, beta, cmin, cmax, cvmaxm, cvmaxp, delta, denom,    &
             dxsign, edgmax, error, f, gamma, pareta, parmu, parsig, phi,     &
             phimin, prerec, prerem, ratio, resmax, resnew, rho, temp, tempa, &
             total, trured, vmnew, vmold, weta, wsig
INTEGER   :: i, ibrnch, iflag, ifull, iptem, iptemp, j, jdrop, k, l, mp,  &
             nbest, nfvals, np

iptem = MIN(n,5)
iptemp = iptem + 1
np = n + 1
mp = m + 1
alpha = 0.25_dp
beta = 2.1_dp
gamma = 0.5_dp
delta = 1.1_dp
rho = rhobeg
parmu = 0.0_dp
IF (iprint >= 2) WRITE(*, 10) rho
10 FORMAT (/'   The initial value of RHO is', G13.6,   &
           '  and PARMU is set to zero.')
nfvals = 0
temp = 1.0_dp/rho
DO i=1,n
  sim(i,np) = x(i)
  DO j=1,n
    sim(i,j) = 0.0_dp
    simi(i,j) = 0.0_dp
  END DO
  sim(i,i) = rho
  simi(i,i) = temp
END DO
jdrop = np
ibrnch = 0

!  Make the next call of the user-supplied subroutine CALCFC. These
!  instructions are also used for calling CALCFC during the iterations of
!  the algorithm.

40 IF (nfvals >= maxfun .AND. nfvals > 0) THEN
  IF (iprint >= 1) WRITE(*, 50)
  50 FORMAT (/'   Return from subroutine COBYLA because the ',  &
             'MAXFUN limit has been reached.')
  GO TO 600
END IF
nfvals = nfvals + 1
CALL calcfc (n, m, x, f, con)
resmax = 0.0_dp
IF (m > 0) THEN
  DO k=1,m
    resmax = MAX(resmax, - con(k))
  END DO
END IF
IF (nfvals == iprint-1 .OR. iprint == 3) THEN
  WRITE(*, 70) nfvals, f, resmax, x(1:iptem)
  70 FORMAT (/'   NFVALS = ', i5, '   F = ', G13.6, '    MAXCV = ',  &
             G13.6/ ('   X = ', 5G14.6))
  IF (iptem < n) WRITE(*, 80) x(iptemp:n)
  80 FORMAT (G19.6, G15.6)
END IF
con(mp) = f
con(mpp) = resmax
IF (ibrnch == 1) GO TO 440

!  Set the recently calculated function values in a column of DATMAT. This
!  array has a column for each vertex of the current simplex, the entries of
!  each column being the values of the constraint functions (if any)
!  followed by the objective function and the greatest constraint violation
!  at the vertex.

DO k=1,mpp
  datmat(k,jdrop) = con(k)
END DO
IF (nfvals > np) GO TO 130

!  Exchange the new vertex of the initial simplex with the optimal vertex if
!  necessary. Then, if the initial simplex is not complete, pick its next
!  vertex and calculate the function values there.

IF (jdrop <= n) THEN
  IF (datmat(mp,np) <= f) THEN
    x(jdrop) = sim(jdrop,np)
  ELSE
    sim(jdrop,np) = x(jdrop)
    DO k=1,mpp
      datmat(k,jdrop) = datmat(k,np)
      datmat(k,np) = con(k)
    END DO
    DO k=1,jdrop
      sim(jdrop,k) = -rho
      temp = -SUM( simi(k:jdrop, k) )
      simi(jdrop,k) = temp
    END DO
  END IF
END IF
IF (nfvals <= n) THEN
  jdrop = nfvals
  x(jdrop) = x(jdrop) + rho
  GO TO 40
END IF
130 ibrnch = 1

!  Identify the optimal vertex of the current simplex.

140 phimin = datmat(mp,np) + parmu*datmat(mpp,np)
nbest = np
DO j=1,n
  temp = datmat(mp,j) + parmu*datmat(mpp,j)
  IF (temp < phimin) THEN
    nbest = j
    phimin = temp
  ELSE IF (temp == phimin .AND. parmu == 0.0_dp) THEN
    IF (datmat(mpp,j) < datmat(mpp,nbest)) nbest = j
  END IF
END DO

!  Switch the best vertex into pole position if it is not there already,
!  and also update SIM, SIMI and DATMAT.

IF (nbest <= n) THEN
  DO i=1,mpp
    temp = datmat(i,np)
    datmat(i,np) = datmat(i,nbest)
    datmat(i,nbest) = temp
  END DO
  DO i=1,n
    temp = sim(i,nbest)
    sim(i,nbest) = 0.0_dp
    sim(i,np) = sim(i,np) + temp
    tempa = 0.0_dp
    DO k=1,n
      sim(i,k) = sim(i,k) - temp
      tempa = tempa - simi(k,i)
    END DO
    simi(nbest,i) = tempa
  END DO
END IF

!  Make an error return if SIGI is a poor approximation to the inverse of
!  the leading N by N submatrix of SIG.

error = 0.0_dp
DO i=1,n
  DO j=1,n
    temp = 0.0_dp
    IF (i == j) temp = temp - 1.0_dp
    temp = temp + DOT_PRODUCT( simi(i,1:n), sim(1:n,j) )
    error = MAX(error, ABS(temp))
  END DO
END DO
IF (error > 0.1_dp) THEN
  IF (iprint >= 1) WRITE(*, 210)
  210 FORMAT (/'   Return from subroutine COBYLA because ',  &
              'rounding errors are becoming damaging.')
  GO TO 600
END IF

!  Calculate the coefficients of the linear approximations to the objective
!  and constraint functions, placing minus the objective function gradient
!  after the constraint gradients in the array A. The vector W is used for
!  working space.

DO k=1,mp
  con(k) = -datmat(k,np)
  DO j=1,n
    w(j) = datmat(k,j) + con(k)
  END DO
  DO i=1,n
    temp = DOT_PRODUCT( w(1:n), simi(1:n,i) )
    IF (k == mp) temp = -temp
    a(i,k) = temp
  END DO
END DO

!  Calculate the values of sigma and eta, and set IFLAG = 0 if the current
!  simplex is not acceptable.

iflag = 1
parsig = alpha*rho
pareta = beta*rho
DO j=1,n
  wsig = SUM( simi(j,1:n)**2 )
  weta = SUM( sim(1:n,j)**2 )
  vsig(j) = 1.0_dp/SQRT(wsig)
  veta(j) = SQRT(weta)
  IF (vsig(j) < parsig .OR. veta(j) > pareta) iflag = 0
END DO

!  If a new vertex is needed to improve acceptability, then decide which
!  vertex to drop from the simplex.

IF (ibrnch == 1 .OR. iflag == 1) GO TO 370
jdrop = 0
temp = pareta
DO j=1,n
  IF (veta(j) > temp) THEN
    jdrop = j
    temp = veta(j)
  END IF
END DO
IF (jdrop == 0) THEN
  DO j=1,n
    IF (vsig(j) < temp) THEN
      jdrop = j
      temp = vsig(j)
    END IF
  END DO
END IF

!  Calculate the step to the new vertex and its sign.

temp = gamma*rho*vsig(jdrop)
dx(1:n) = temp*simi(jdrop,1:n)
cvmaxp = 0.0_dp
cvmaxm = 0.0_dp
DO k=1,mp
  total = DOT_PRODUCT( a(1:n,k), dx(1:n) )
  IF (k < mp) THEN
    temp = datmat(k,np)
    cvmaxp = MAX(cvmaxp, -total - temp)
    cvmaxm = MAX(cvmaxm, total - temp)
  END IF
END DO
dxsign = 1.0_dp
IF (parmu*(cvmaxp - cvmaxm) > total + total) dxsign = -1.0_dp

!  Update the elements of SIM and SIMI, and set the next X.

temp = 0.0_dp
DO i=1,n
  dx(i) = dxsign*dx(i)
  sim(i,jdrop) = dx(i)
  temp = temp + simi(jdrop,i)*dx(i)
END DO
simi(jdrop,1:n) = simi(jdrop,1:n) / temp
DO j=1,n
  IF (j /= jdrop) THEN
    temp = DOT_PRODUCT( simi(j,1:n), dx(1:n) )
    simi(j,1:n) = simi(j,1:n) - temp*simi(jdrop,1:n)
  END IF
  x(j) = sim(j,np) + dx(j)
END DO
GO TO 40

!  Calculate DX = x(*)-x(0).
!  Branch if the length of DX is less than 0.5*RHO.

370 CALL trstlp (n, m, a, con, rho, dx, ifull)
IF (ifull == 0) THEN
  temp = SUM( dx(1:n)**2 )
  IF (temp < 0.25_dp*rho*rho) THEN
    ibrnch = 1
    GO TO 550
  END IF
END IF

!  Predict the change to F and the new maximum constraint violation if the
!  variables are altered from x(0) to x(0) + DX.

resnew = 0.0_dp
con(mp) = 0.0_dp
DO k=1,mp
  total = con(k) - DOT_PRODUCT( a(1:n,k), dx(1:n) )
  IF (k < mp) resnew = MAX(resnew, total)
END DO

!  Increase PARMU if necessary and branch back if this change alters the
!  optimal vertex. Otherwise PREREM and PREREC will be set to the predicted
!  reductions in the merit function and the maximum constraint violation
!  respectively.

barmu = 0.0_dp
prerec = datmat(mpp,np) - resnew
IF (prerec > 0.0_dp) barmu = total/prerec
IF (parmu < 1.5_dp*barmu) THEN
  parmu = 2.0_dp*barmu
  IF (iprint >= 2) WRITE(*, 410) parmu
  410 FORMAT (/'   Increase in PARMU to', G13.6)
  phi = datmat(mp,np) + parmu*datmat(mpp,np)
  DO j=1,n
    temp = datmat(mp,j) + parmu*datmat(mpp,j)
    IF (temp < phi) GO TO 140
    IF (temp == phi .AND. parmu == 0.0) THEN
      IF (datmat(mpp,j) < datmat(mpp,np)) GO TO 140
    END IF
  END DO
END IF
prerem = parmu*prerec - total

!  Calculate the constraint and objective functions at x(*).
!  Then find the actual reduction in the merit function.

x(1:n) = sim(1:n,np) + dx(1:n)
ibrnch = 1
GO TO 40

440 vmold = datmat(mp,np) + parmu*datmat(mpp,np)
vmnew = f + parmu*resmax
trured = vmold - vmnew
IF (parmu == 0.0_dp .AND. f == datmat(mp,np)) THEN
  prerem = prerec
  trured = datmat(mpp,np) - resmax
END IF

!  Begin the operations that decide whether x(*) should replace one of the
!  vertices of the current simplex, the change being mandatory if TRURED is
!  positive. Firstly, JDROP is set to the index of the vertex that is to be
!  replaced.

ratio = 0.0_dp
IF (trured <= 0.0) ratio = 1.0
jdrop = 0
DO j=1,n
  temp = DOT_PRODUCT( simi(j,1:n), dx(1:n) )
  temp = ABS(temp)
  IF (temp > ratio) THEN
    jdrop = j
    ratio = temp
  END IF
  sigbar(j) = temp*vsig(j)
END DO

!  Calculate the value of ell.

edgmax = delta*rho
l = 0
DO j=1,n
  IF (sigbar(j) >= parsig .OR. sigbar(j) >= vsig(j)) THEN
    temp = veta(j)
    IF (trured > 0.0_dp) THEN
      temp = SUM( (dx(1:n) - sim(1:n,j))**2 )
      temp = SQRT(temp)
    END IF
    IF (temp > edgmax) THEN
      l = j
      edgmax = temp
    END IF
  END IF
END DO
IF (l > 0) jdrop = l
IF (jdrop == 0) GO TO 550

!  Revise the simplex by updating the elements of SIM, SIMI and DATMAT.

temp = 0.0_dp
DO i=1,n
  sim(i,jdrop) = dx(i)
  temp = temp + simi(jdrop,i)*dx(i)
END DO
simi(jdrop,1:n) = simi(jdrop,1:n) / temp
DO j=1,n
  IF (j /= jdrop) THEN
    temp = DOT_PRODUCT( simi(j,1:n), dx(1:n) )
    simi(j,1:n) = simi(j,1:n) - temp*simi(jdrop,1:n)
  END IF
END DO
datmat(1:mpp,jdrop) = con(1:mpp)

!  Branch back for further iterations with the current RHO.

IF (trured > 0.0_dp .AND. trured >= 0.1_dp*prerem) GO TO 140
550 IF (iflag == 0) THEN
  ibrnch = 0
  GO TO 140
END IF

!  Otherwise reduce RHO if it is not at its least value and reset PARMU.

IF (rho > rhoend) THEN
  rho = 0.5_dp*rho
  IF (rho <= 1.5_dp*rhoend) rho = rhoend
  IF (parmu > 0.0_dp) THEN
    denom = 0.0_dp
    DO k=1,mp
      cmin = datmat(k,np)
      cmax = cmin
      DO i=1,n
        cmin = MIN(cmin, datmat(k,i))
        cmax = MAX(cmax, datmat(k,i))
      END DO
      IF (k <= m .AND. cmin < 0.5_dp*cmax) THEN
        temp = MAX(cmax,0.0_dp) - cmin
        IF (denom <= 0.0_dp) THEN
          denom = temp
        ELSE
          denom = MIN(denom,temp)
        END IF
      END IF
    END DO
    IF (denom == 0.0_dp) THEN
      parmu = 0.0_dp
    ELSE IF (cmax - cmin < parmu*denom) THEN
      parmu = (cmax - cmin)/denom
    END IF
  END IF
  IF (iprint >= 2) WRITE(*, 580) rho,parmu
  580 FORMAT (/'   Reduction in RHO to ', G13.6, '  and PARMU = ', G13.6)
  IF (iprint == 2) THEN
    WRITE(*, 70) nfvals, datmat(mp,np), datmat(mpp,np), sim(1:iptem,np)
    IF (iptem < n) WRITE(*, 80) x(iptemp:n)
  END IF
  GO TO 140
END IF

!  Return the best calculated values of the variables.

IF (iprint >= 1) WRITE(*, 590)
590 FORMAT (/'   Normal return from subroutine COBYLA')
IF (ifull == 1) GO TO 620

600 x(1:n) = sim(1:n,np)
f = datmat(mp,np)
resmax = datmat(mpp,np)
620 IF (iprint >= 1) THEN
  WRITE(*, 70) nfvals, f, resmax, x(1:iptem)
  IF (iptem < n) WRITE(*, 80) x(iptemp:n)
END IF
maxfun = nfvals

RETURN
END SUBROUTINE cobylb
!------------------------------------------------------------------------------

SUBROUTINE trstlp (n, m, a, b, rho, dx, ifull)

! N.B. Arguments Z, ZDOTA, VMULTC, SDIRN, DXNEW, VMULTD & IACT have been removed.

INTEGER, INTENT(IN)     :: n
INTEGER, INTENT(IN)     :: m
REAL (dp), INTENT(IN)   :: a(:,:)
REAL (dp), INTENT(IN)   :: b(:)
REAL (dp), INTENT(IN)   :: rho
REAL (dp), INTENT(OUT)  :: dx(:)
INTEGER, INTENT(OUT)    :: ifull

!  This subroutine calculates an N-component vector DX by applying the
!  following two stages. In the first stage, DX is set to the shortest
!  vector that minimizes the greatest violation of the constraints
!    A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K = 2,3,...,M,
!  subject to the Euclidean length of DX being at most RHO. If its length is
!  strictly less than RHO, then we use the resultant freedom in DX to
!  minimize the objective function
!           -A(1,M+1)*DX(1) - A(2,M+1)*DX(2) - ... - A(N,M+1)*DX(N)
!  subject to no increase in any greatest constraint violation. This
!  notation allows the gradient of the objective function to be regarded as
!  the gradient of a constraint. Therefore the two stages are distinguished
!  by MCON .EQ. M and MCON .GT. M respectively. It is possible that a
!  degeneracy may prevent DX from attaining the target length RHO. Then the
!  value IFULL = 0 would be set, but usually IFULL = 1 on return.

!  In general NACT is the number of constraints in the active set and
!  IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT
!  contains a permutation of the remaining constraint indices.  Further, Z
!  is an orthogonal matrix whose first NACT columns can be regarded as the
!  result of Gram-Schmidt applied to the active constraint gradients.  For
!  J = 1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th
!  column of Z with the gradient of the J-th active constraint.  DX is the
!  current vector of variables and here the residuals of the active
!  constraints should be zero. Further, the active constraints have
!  nonnegative Lagrange multipliers that are held at the beginning of
!  VMULTC. The remainder of this vector holds the residuals of the inactive
!  constraints at DX, the ordering of the components of VMULTC being in
!  agreement with the permutation of the indices of the constraints that is
!  in IACT. All these residuals are nonnegative, which is achieved by the
!  shift RESMAX that makes the least residual zero.

!  Initialize Z and some other variables. The value of RESMAX will be
!  appropriate to DX = 0, while ICON will be the index of a most violated
!  constraint if RESMAX is positive. Usually during the first stage the
!  vector SDIRN gives a search direction that reduces all the active
!  constraint violations by one simultaneously.

! Local variables

REAL (dp) :: z(n,n), zdota(m+1), vmultc(m+1), sdirn(n), dxnew(n), vmultd(m+1)
REAL (dp) :: acca, accb, alpha, beta, dd, optnew, optold, ratio, resmax,   &
             resold, sd, sp, spabs, ss, step, stpful, sumabs, temp, tempa, &
             tot, total, vsave, zdotv, zdotw, zdvabs, zdwabs
INTEGER   :: i, iact(m+1), icon, icount, isave, k, kk, kl, kp, kw, mcon,   &
             nact, nactx

ifull = 1
mcon = m
nact = 0
resmax = 0.0_dp
DO i=1,n
  z(i,1:n) = 0.0_dp
  z(i,i) = 1.0_dp
  dx(i) = 0.0_dp
END DO
IF (m >= 1) THEN
  DO k=1,m
    IF (b(k) > resmax) THEN
      resmax = b(k)
      icon = k
    END IF
  END DO
  DO k=1,m
    iact(k) = k
    vmultc(k) = resmax - b(k)
  END DO
END IF
IF (resmax == 0.0_dp) GO TO 480
sdirn(1:n) = 0.0_dp

!  End the current stage of the calculation if 3 consecutive iterations
!  have either failed to reduce the best calculated value of the objective
!  function or to increase the number of active constraints since the best
!  value was calculated. This strategy prevents cycling, but there is a
!  remote possibility that it will cause premature termination.

60 optold = 0.0_dp
icount = 0
70 IF (mcon == m) THEN
  optnew = resmax
ELSE
  optnew = - DOT_PRODUCT( dx(1:n), a(1:n,mcon) )
END IF
IF (icount == 0 .OR. optnew < optold) THEN
  optold = optnew
  nactx = nact
  icount = 3
ELSE IF (nact > nactx) THEN
  nactx = nact
  icount = 3
ELSE
  icount = icount - 1
  IF (icount == 0) GO TO 490
END IF

!  If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to
!  the active set. Apply Givens rotations so that the last N-NACT-1 columns
!  of Z are orthogonal to the gradient of the new constraint, a scalar
!  product being set to zero if its nonzero value could be due to computer
!  rounding errors. The array DXNEW is used for working space.

IF (icon <= nact) GO TO 260
kk = iact(icon)
dxnew(1:n) = a(1:n,kk)
tot = 0.0_dp
k = n
100 IF (k > nact) THEN
  sp = 0.0_dp
  spabs = 0.0_dp
  DO i=1,n
    temp = z(i,k)*dxnew(i)
    sp = sp + temp
    spabs = spabs + ABS(temp)
  END DO
  acca = spabs + 0.1_dp*ABS(sp)
  accb = spabs + 0.2_dp*ABS(sp)
  IF (spabs >= acca .OR. acca >= accb) sp = 0.0_dp
  IF (tot == 0.0_dp) THEN
    tot = sp
  ELSE
    kp = k + 1
    temp = SQRT(sp*sp + tot*tot)
    alpha = sp/temp
    beta = tot/temp
    tot = temp
    DO i=1,n
      temp = alpha*z(i,k) + beta*z(i,kp)
      z(i,kp) = alpha*z(i,kp) - beta*z(i,k)
      z(i,k) = temp
    END DO
  END IF
  k = k - 1
  GO TO 100
END IF

!  Add the new constraint if this can be done without a deletion from the
!  active set.

IF (tot /= 0.0_dp) THEN
  nact = nact + 1
  zdota(nact) = tot
  vmultc(icon) = vmultc(nact)
  vmultc(nact) = 0.0_dp
  GO TO 210
END IF

!  The next instruction is reached if a deletion has to be made from the
!  active set in order to make room for the new active constraint, because
!  the new constraint gradient is a linear combination of the gradients of
!  the old active constraints.  Set the elements of VMULTD to the multipliers
!  of the linear combination.  Further, set IOUT to the index of the
!  constraint to be deleted, but branch if no suitable index can be found.

ratio = -1.0_dp
k = nact
130 zdotv = 0.0_dp
zdvabs = 0.0_dp
DO i=1,n
  temp = z(i,k)*dxnew(i)
  zdotv = zdotv + temp
  zdvabs = zdvabs + ABS(temp)
END DO
acca = zdvabs + 0.1_dp*ABS(zdotv)
accb = zdvabs + 0.2_dp*ABS(zdotv)
IF (zdvabs < acca .AND. acca < accb) THEN
  temp = zdotv/zdota(k)
  IF (temp > 0.0_dp .AND. iact(k) <= m) THEN
    tempa = vmultc(k)/temp
    IF (ratio < 0.0_dp .OR. tempa < ratio) THEN
      ratio = tempa
    END IF
  END IF
  IF (k >= 2) THEN
    kw = iact(k)
    dxnew(1:n) = dxnew(1:n) - temp*a(1:n,kw)
  END IF
  vmultd(k) = temp
ELSE
  vmultd(k) = 0.0_dp
END IF
k = k - 1
IF (k > 0) GO TO 130
IF (ratio < 0.0_dp) GO TO 490

!  Revise the Lagrange multipliers and reorder the active constraints so
!  that the one to be replaced is at the end of the list. Also calculate the
!  new value of ZDOTA(NACT) and branch if it is not acceptable.

DO k=1,nact
  vmultc(k) = MAX(0.0_dp,vmultc(k) - ratio*vmultd(k))
END DO
IF (icon < nact) THEN
  isave = iact(icon)
  vsave = vmultc(icon)
  k = icon
  170 kp = k + 1
  kw = iact(kp)
  sp = DOT_PRODUCT( z(1:n,k), a(1:n,kw) )
  temp = SQRT(sp*sp + zdota(kp)**2)
  alpha = zdota(kp)/temp
  beta = sp/temp
  zdota(kp) = alpha*zdota(k)
  zdota(k) = temp
  DO i=1,n
    temp = alpha*z(i,kp) + beta*z(i,k)
    z(i,kp) = alpha*z(i,k) - beta*z(i,kp)
    z(i,k) = temp
  END DO
  iact(k) = kw
  vmultc(k) = vmultc(kp)
  k = kp
  IF (k < nact) GO TO 170
  iact(k) = isave
  vmultc(k) = vsave
END IF
temp = DOT_PRODUCT( z(1:n,nact), a(1:n,kk) )
IF (temp == 0.0_dp) GO TO 490
zdota(nact) = temp
vmultc(icon) = 0.0_dp
vmultc(nact) = ratio

!  Update IACT and ensure that the objective function continues to be
!  treated as the last active constraint when MCON>M.

210 iact(icon) = iact(nact)
iact(nact) = kk
IF (mcon > m .AND. kk /= mcon) THEN
  k = nact - 1
  sp = DOT_PRODUCT( z(1:n,k), a(1:n,kk) )
  temp = SQRT(sp*sp + zdota(nact)**2)
  alpha = zdota(nact)/temp
  beta = sp/temp
  zdota(nact) = alpha*zdota(k)
  zdota(k) = temp
  DO i=1,n
    temp = alpha*z(i,nact) + beta*z(i,k)
    z(i,nact) = alpha*z(i,k) - beta*z(i,nact)
    z(i,k) = temp
  END DO
  iact(nact) = iact(k)
  iact(k) = kk
  temp = vmultc(k)
  vmultc(k) = vmultc(nact)
  vmultc(nact) = temp
END IF

!  If stage one is in progress, then set SDIRN to the direction of the next
!  change to the current vector of variables.

IF (mcon > m) GO TO 320
kk = iact(nact)
temp = DOT_PRODUCT( sdirn(1:n), a(1:n,kk) )
temp = temp - 1.0_dp
temp = temp/zdota(nact)
sdirn(1:n) = sdirn(1:n) - temp*z(1:n,nact)
GO TO 340

!  Delete the constraint that has the index IACT(ICON) from the active set.

260 IF (icon < nact) THEN
  isave = iact(icon)
  vsave = vmultc(icon)
  k = icon
  DO
    kp = k + 1
    kk = iact(kp)
    sp = DOT_PRODUCT( z(1:n,k), a(1:n,kk) )
    temp = SQRT(sp*sp + zdota(kp)**2)
    alpha = zdota(kp)/temp
    beta = sp/temp
    zdota(kp) = alpha*zdota(k)
    zdota(k) = temp
    DO i=1,n
      temp = alpha*z(i,kp) + beta*z(i,k)
      z(i,kp) = alpha*z(i,k) - beta*z(i,kp)
      z(i,k) = temp
    END DO
    iact(k) = kk
    vmultc(k) = vmultc(kp)
    k = kp
    IF (k >= nact) EXIT
  END DO
  iact(k) = isave
  vmultc(k) = vsave
END IF
nact = nact - 1

!  If stage one is in progress, then set SDIRN to the direction of the next
!  change to the current vector of variables.

IF (mcon > m) GO TO 320
temp = DOT_PRODUCT( sdirn(1:n), z(1:n,nact+1) )
sdirn(1:n) = sdirn(1:n) - temp*z(1:n,nact+1)
GO TO 340

!  Pick the next search direction of stage two.

320 temp = 1.0_dp/zdota(nact)
sdirn(1:n) = temp*z(1:n,nact)

!  Calculate the step to the boundary of the trust region or take the step
!  that reduces RESMAX to zero. The two statements below that include the
!  factor 1.0E-6 prevent some harmless underflows that occurred in a test
!  calculation. Further, we skip the step if it could be zero within a
!  reasonable tolerance for computer rounding errors.

340 dd = rho*rho
sd = 0.0_dp
ss = 0.0_dp
DO i=1,n
  IF (ABS(dx(i)) >= 1.0E-6*rho) dd = dd - dx(i)**2
  sd = sd + dx(i)*sdirn(i)
  ss = ss + sdirn(i)**2
END DO
IF (dd <= 0.0_dp) GO TO 490
temp = SQRT(ss*dd)
IF (ABS(sd) >= 1.0E-6*temp) temp = SQRT(ss*dd + sd*sd)
stpful = dd/(temp + sd)
step = stpful
IF (mcon == m) THEN
  acca = step + 0.1_dp*resmax
  accb = step + 0.2_dp*resmax
  IF (step >= acca .OR. acca >= accb) GO TO 480
  step = MIN(step,resmax)
END IF

!  Set DXNEW to the new variables if STEP is the steplength, and reduce
!  RESMAX to the corresponding maximum residual if stage one is being done.
!  Because DXNEW will be changed during the calculation of some Lagrange
!  multipliers, it will be restored to the following value later.

dxnew(1:n) = dx(1:n) + step*sdirn(1:n)
IF (mcon == m) THEN
  resold = resmax
  resmax = 0.0_dp
  DO k=1,nact
    kk = iact(k)
    temp = b(kk) - DOT_PRODUCT( a(1:n,kk), dxnew(1:n) )
    resmax = MAX(resmax,temp)
  END DO
END IF

!  Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A
!  device is included to force VMULTD(K) = 0.0 if deviations from this value
!  can be attributed to computer rounding errors. First calculate the new
!  Lagrange multipliers.

k = nact
390 zdotw = 0.0_dp
zdwabs = 0.0_dp
DO i=1,n
  temp = z(i,k)*dxnew(i)
  zdotw = zdotw + temp
  zdwabs = zdwabs + ABS(temp)
END DO
acca = zdwabs + 0.1_dp*ABS(zdotw)
accb = zdwabs + 0.2_dp*ABS(zdotw)
IF (zdwabs >= acca .OR. acca >= accb) zdotw = 0.0_dp
vmultd(k) = zdotw / zdota(k)
IF (k >= 2) THEN
  kk = iact(k)
  dxnew(1:n) = dxnew(1:n) - vmultd(k)*a(1:n,kk)
  k = k - 1
  GO TO 390
END IF
IF (mcon > m) vmultd(nact) = MAX(0.0_dp,vmultd(nact))

!  Complete VMULTC by finding the new constraint residuals.

dxnew(1:n) = dx(1:n) + step*sdirn(1:n)
IF (mcon > nact) THEN
  kl = nact + 1
  DO k=kl,mcon
    kk = iact(k)
    total = resmax - b(kk)
    sumabs = resmax + ABS(b(kk))
    DO i=1,n
      temp = a(i,kk)*dxnew(i)
      total = total + temp
      sumabs = sumabs + ABS(temp)
    END DO
    acca = sumabs + 0.1*ABS(total)
    accb = sumabs + 0.2*ABS(total)
    IF (sumabs >= acca .OR. acca >= accb) total = 0.0
    vmultd(k) = total
  END DO
END IF

!  Calculate the fraction of the step from DX to DXNEW that will be taken.

ratio = 1.0_dp
icon = 0
DO k=1,mcon
  IF (vmultd(k) < 0.0_dp) THEN
    temp = vmultc(k)/(vmultc(k) - vmultd(k))
    IF (temp < ratio) THEN
      ratio = temp
      icon = k
    END IF
  END IF
END DO

!  Update DX, VMULTC and RESMAX.

temp = 1.0_dp - ratio
dx(1:n) = temp*dx(1:n) + ratio*dxnew(1:n)
DO k=1,mcon
  vmultc(k) = MAX(0.0_dp,temp*vmultc(k) + ratio*vmultd(k))
END DO
IF (mcon == m) resmax = resold + ratio*(resmax - resold)

!  If the full step is not acceptable then begin another iteration.
!  Otherwise switch to stage two or end the calculation.

IF (icon > 0) GO TO 70
IF (step == stpful) GO TO 500
480 mcon = m + 1
icon = mcon
iact(mcon) = mcon
vmultc(mcon) = 0.0_dp
GO TO 60

!  We employ any freedom that may be available to reduce the objective
!  function before returning a DX whose length is less than RHO.

490 IF (mcon == m) GO TO 480
ifull = 0

500 RETURN
END SUBROUTINE trstlp

END MODULE cobyla2



MODULE common_nprob
IMPLICIT NONE

INTEGER, SAVE, PUBLIC :: nprob

END MODULE common_nprob


!------------------------------------------------------------------------------
!  Main program of test problems in Report DAMTP 1992/NA5.
!------------------------------------------------------------------------------
PROGRAM test_cobyla
USE common_nprob
USE cobyla2
IMPLICIT NONE

INTEGER, PARAMETER :: nn = 10

REAL (dp) :: rhobeg, rhoend, temp, tempa, tempb, tempc, tempd, x(nn),  &
             xopt(nn)
INTEGER   :: i, icase, iprint, m, maxfun, n

INTERFACE
  SUBROUTINE calcfc (n, m, x, f, con)
    IMPLICIT NONE
    INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)
    INTEGER, INTENT(IN)    :: n, m
    REAL (dp), INTENT(IN)  :: x(:)
    REAL (dp), INTENT(OUT) :: f
    REAL (dp), INTENT(OUT) :: con(:)
  END SUBROUTINE calcfc
END INTERFACE

DO nprob=1,10
  IF (nprob == 1) THEN

!  Minimization of a simple quadratic function of two variables.

    WRITE(*, 10)
    10 FORMAT (/'       Output from test problem 1 (Simple quadratic)')
    n = 2
    m = 0
    xopt(1) = -1.0_dp
    xopt(2) = 0.0_dp
  ELSE IF (nprob == 2) THEN

!  Easy two dimensional minimization in unit circle.

    WRITE(*, 20)
    20 FORMAT (/'       Output from test problem 2 (2D unit circle calculation)')
    n = 2
    m = 1
    xopt(1) = SQRT(0.5_dp)
    xopt(2) = -xopt(1)
  ELSE IF (nprob == 3) THEN

!  Easy three dimensional minimization in ellipsoid.

    WRITE(*, 30)
    30 FORMAT (/'       Output from test problem 3 (3D ellipsoid calculation)')
    n = 3
    m = 1
    xopt(1) = 1.0_dp/SQRT(3.0_dp)
    xopt(2) = 1.0_dp/SQRT(6.0_dp)
    xopt(3) = -1.0_dp/3.0_dp
  ELSE IF (nprob == 4) THEN

!  Weak version of Rosenbrock's problem.

    WRITE(*, 40)
    40 FORMAT (/'       Output from test problem 4 (Weak Rosenbrock)')
    n = 2
    m = 0
    xopt(1) = -1.0_dp
    xopt(2) = 1.0_dp
  ELSE IF (nprob == 5) THEN

!  Intermediate version of Rosenbrock's problem.

    WRITE(*, 50)
    50 FORMAT (/'       Output from test problem 5 (Intermediate Rosenbrock)')
    n = 2
    m = 0
    xopt(1) = -1.0_dp
    xopt(2) = 1.0_dp
  ELSE IF (nprob == 6) THEN

!  This problem is taken from Fletcher's book Practical Methods of
!  Optimization and has the equation number (9.1.15).

    WRITE(*, 60)
    60 FORMAT (/'       Output from test problem 6 (Equation ',  &
               '(9.1.15) in Fletcher)')
    n = 2
    m = 2
    xopt(1) = SQRT(0.5_dp)
    xopt(2) = xopt(1)
  ELSE IF (nprob == 7) THEN

!  This problem is taken from Fletcher's book Practical Methods of
!  Optimization and has the equation number (14.4.2).

    WRITE(*, 70)
    70 FORMAT (/'       Output from test problem 7 (Equation ',  &
               '(14.4.2) in Fletcher)')
    n = 3
    m = 3
    xopt(1) = 0.0_dp
    xopt(2) = -3.0_dp
    xopt(3) = -3.0_dp
  ELSE IF (nprob == 8) THEN

!  This problem is taken from page 66 of Hock and Schittkowski's book Test
!  Examples for Nonlinear Programming Codes. It is their test problem Number
!  43, and has the name Rosen-Suzuki.

    WRITE(*, 80)
    80 FORMAT (/'       Output from test problem 8 (Rosen-Suzuki)')
    n=4
    m=3
    xopt(1) = 0.0_dp
    xopt(2) = 1.0_dp
    xopt(3) = 2.0_dp
    xopt(4) = -1.0_dp
  ELSE IF (nprob == 9) THEN

!  This problem is taken from page 111 of Hock and Schittkowski's
!  book Test Examples for Nonlinear Programming Codes. It is their
!  test problem Number 100.

    WRITE(*, 90)
    90 FORMAT (/'       Output from test problem 9 (Hock and Schittkowski 100)')
    n = 7
    m = 4
    xopt(1) = 2.330499_dp
    xopt(2) = 1.951372_dp
    xopt(3) = -0.4775414_dp
    xopt(4) = 4.365726_dp
    xopt(5) = -0.624487_dp
    xopt(6) = 1.038131_dp
    xopt(7) = 1.594227_dp
  ELSE IF (nprob == 10) THEN

!  This problem is taken from page 415 of Luenberger's book Applied
!  Nonlinear Programming. It is to maximize the area of a hexagon of
!  unit diameter.

    WRITE(*, 100)
    100 FORMAT (/'       Output from test problem 10 (Hexagon area)')
    n = 9
    m = 14
  END IF

  DO icase = 1,2
    x(1:n) = 1.0_dp
    rhobeg = 0.5_dp
    rhoend = 1.d-6
    IF (icase == 2) rhoend = 1.d-8
    iprint = 1
    maxfun = 3500
    CALL cobyla (n, m, x, rhobeg, rhoend, iprint, maxfun)
    IF (nprob == 10) THEN
      tempa = x(1) + x(3) + x(5) + x(7)
      tempb = x(2) + x(4) + x(6) + x(8)
      tempc = 0.5_dp/SQRT(tempa*tempa + tempb*tempb)
      tempd = tempc*SQRT(3.0_dp)
      xopt(1) = tempd*tempa + tempc*tempb
      xopt(2) = tempd*tempb - tempc*tempa
      xopt(3) = tempd*tempa - tempc*tempb
      xopt(4) = tempd*tempb + tempc*tempa
      DO i=1,4
        xopt(i+4) = xopt(i)
      END DO
    END IF
    temp = 0.0_dp
    DO i=1,n
      temp = temp + (x(i) - xopt(i))**2
    END DO
    WRITE(*, 150) SQRT(temp)
    150 FORMAT (/'     Least squares error in variables = ', G16.6)
  END DO
  WRITE(*, 170)
  170 FORMAT ('  ----------------------------------------------',  &
              '--------------------')
END DO

STOP
END PROGRAM test_cobyla


!------------------------------------------------------------------------------

SUBROUTINE calcfc (n, m, x, f, con)
USE common_nprob
IMPLICIT NONE

INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)

INTEGER, INTENT(IN)     :: n
INTEGER, INTENT(IN)     :: m
REAL (dp), INTENT(IN)   :: x(:)
REAL (dp), INTENT(OUT)  :: f
REAL (dp), INTENT(OUT)  :: con(:)

IF (nprob == 1) THEN

!  Test problem 1 (Simple quadratic)

  f = 10.0_dp*(x(1) + 1.0_dp)**2 + x(n)**2
ELSE IF (nprob == 2) THEN

!    Test problem 2 (2D unit circle calculation)

  f = x(1)*x(2)
  con(m) = 1.0_dp - x(1)**2 - x(n)**2
ELSE IF (nprob == 3) THEN

!  Test problem 3 (3D ellipsoid calculation)

  f = x(1)*x(2)*x(3)
  con(1) = 1.0_dp - x(1)**2 - 2.0_dp*x(2)**2 - 3.0_dp*x(n)**2
ELSE IF (nprob == 4) THEN

!  Test problem 4 (Weak Rosenbrock)

  f = (x(1)**2 - x(2))**2 + (1.0_dp + x(1))**2
ELSE IF (nprob == 5) THEN

!  Test problem 5 (Intermediate Rosenbrock)

  f = 10.0_dp*(x(1)**2 - x(n))**2 + (1.0_dp + x(1))**2
ELSE IF (nprob == 6) THEN

!  Test problem 6 (Equation (9.1.15) in Fletcher's book)

  f = - x(1) - x(2)
  con(1) = x(2) - x(1)**2
  con(2) = 1.0_dp - x(1)**2 - x(2)**2
ELSE IF (nprob == 7) THEN

!  Test problem 7 (Equation (14.4.2) in Fletcher's book)

  f = x(3)
  con(1) = 5.0_dp*x(1) - x(2) + x(3)
  con(2) = x(3) - x(1)**2 - x(2)**2 - 4.0_dp*x(2)
  con(m) = x(3) - 5.0_dp*x(1) - x(2)
ELSE IF (nprob == 8) THEN

!  Test problem 8 (Rosen-Suzuki)

  f = x(1)**2 + x(2)**2 + 2.0*x(3)**2 + x(4)**2 - 5.0_dp*x(1) - 5.0_dp*x(2)  &
      - 21.0_dp*x(3) + 7.0_dp*x(4)
  con(1) = 8.0_dp - x(1)**2 - x(2)**2 - x(3)**2 - x(4)**2 - x(1) + x(2)  &
           - x(3) + x(4)
  con(2) = 10._dp - x(1)**2 - 2._dp*x(2)**2 - x(3)**2 - 2._dp*x(4)**2 + x(1) + x(4)
  con(m) = 5.0_dp - 2.0*x(1)**2 - x(2)**2 - x(3)**2 - 2.0_dp*x(1) + x(2) + x(4)
ELSE IF (nprob == 9) THEN

!  Test problem 9 (Hock and Schittkowski 100)

  f = (x(1) - 10._dp)**2 + 5._dp*(x(2) - 12._dp)**2 + x(3)**4 + 3._dp*(x(4)  &
       - 11._dp)**2 + 10._dp*x(5)**6 + 7._dp*x(6)**2 + x(7)**4 - 4._dp*x(6)*x(7) &
       - 10._dp*x(6) - 8._dp*x(7)
  con(1) = 127._dp - 2._dp*x(1)**2 - 3._dp*x(2)**4 - x(3) - 4._dp*x(4)**2   &
           - 5._dp*x(5)
  con(2) = 282._dp - 7._dp*x(1) - 3._dp*x(2) - 10._dp*x(3)**2 - x(4) + x(5)
  con(3) = 196._dp - 23._dp*x(1) - x(2)**2 - 6._dp*x(6)**2 + 8._dp*x(7)
  con(4) = - 4._dp*x(1)**2 - x(2)**2 + 3._dp*x(1)*x(2) - 2._dp*x(3)**2 - 5._dp*x(6)  &
           + 11._dp*x(7)
ELSE IF (nprob == 10) THEN

!  Test problem 10 (Hexagon area)

  f = - 0.5_dp*(x(1)*x(4) - x(2)*x(3) + x(3)*x(n) - x(5)*x(n) + x(5)*x(8) &
                - x(6)*x(7))
  con(1) = 1.0_dp - x(3)**2 - x(4)**2
  con(2) = 1.0_dp - x(n)**2
  con(3) = 1.0_dp - x(5)**2 - x(6)**2
  con(4) = 1.0_dp - x(1)**2 - (x(2) - x(n))**2
  con(5) = 1.0_dp - (x(1) - x(5))**2 - (x(2) - x(6))**2
  con(6) = 1.0_dp - (x(1) - x(7))**2 - (x(2) - x(8))**2
  con(7) = 1.0_dp - (x(3) - x(5))**2 - (x(4) - x(6))**2
  con(8) = 1.0_dp - (x(3) - x(7))**2 - (x(4) - x(8))**2
  con(9) = 1.0_dp - x(7)**2 - (x(8) - x(n))**2
  con(10) = x(1)*x(4) - x(2)*x(3)
  con(11) = x(3)*x(n)
  con(12) = - x(5)*x(n)
  con(13) = x(5)*x(8) - x(6)*x(7)
  con(m) = x(n)
END IF

RETURN
END SUBROUTINE calcfc



! The following output was obtained using Lahey's ELF90 compiler
!
!
!   Output from test problem 1 (Simple quadratic)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =    87   F =  0.118975E-10    MAXCV =  0.000000
!  X =   -1.00000      0.262770E-05
!
!    Least squares error in variables =     0.272104E-05
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   116   F =  0.258086E-14    MAXCV =  0.000000
!  X =   -1.00000      0.403932E-07
!
!    Least squares error in variables =     0.415516E-07
! ------------------------------------------------------------------
!
!   Output from test problem 2 (2D unit circle calculation)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =    66   F = -0.500000        MAXCV =  0.199977E-11
!  X =   0.707106     -0.707108
!
!    Least squares error in variables =     0.131782E-05
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =    82   F = -0.500000        MAXCV =  0.305826E-15
!  X =   0.707107     -0.707107
!
!    Least squares error in variables =     0.348402E-08
! ------------------------------------------------------------------
!
!   Output from test problem 3 (3D ellipsoid calculation)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =    94   F = -0.785674E-01    MAXCV =  0.593595E-11
!  X =   0.577351      0.408247     -0.333334
!
!    Least squares error in variables =     0.111077E-05
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   115   F = -0.785674E-01    MAXCV =  0.333853E-15
!  X =   0.577350      0.408248     -0.333333
!
!    Least squares error in variables =     0.693448E-08
! ------------------------------------------------------------------
!
!   Output from test problem 4 (Weak Rosenbrock)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   269   F =  0.703758E-10    MAXCV =  0.000000
!  X =  -0.999992      0.999983
!
!    Least squares error in variables =     0.188469E-04
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   387   F =  0.799591E-14    MAXCV =  0.000000
!  X =   -1.00000       1.00000
!
!    Least squares error in variables =     0.203843E-06
! ------------------------------------------------------------------
!
!   Output from test problem 5 (Intermediate Rosenbrock)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =  2174   F =  0.839835E-08    MAXCV =  0.000000
!  X =  -0.999908      0.999816
!
!    Least squares error in variables =     0.205844E-03
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =  3377   F =  0.948763E-12    MAXCV =  0.000000
!  X =  -0.999999      0.999998
!
!    Least squares error in variables =     0.218951E-05
! ------------------------------------------------------------------
!
!   Output from test problem 6 (Equation (9.1.15) in Fletcher)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =    63   F =  -1.41421        MAXCV =  0.199977E-11
!  X =   0.707106      0.707107
!
!    Least squares error in variables =     0.588093E-06
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =    79   F =  -1.41421        MAXCV =  0.138669E-15
!  X =   0.707107      0.707107
!
!    Least squares error in variables =     0.139562E-08
! ------------------------------------------------------------------
!
!   Output from test problem 7 (Equation (14.4.2) in Fletcher)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =    40   F =  -3.00000        MAXCV =  0.000000
!  X =  -0.644464E-20  -3.00000      -3.00000
!
!    Least squares error in variables =     0.299666E-08
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =    47   F =  -3.00000        MAXCV =  0.424972E-08
!  X =  -0.651433E-20  -3.00000      -3.00000
!
!    Least squares error in variables =     0.200334E-08
! ------------------------------------------------------------------
!
!   Output from test problem 8 (Rosen-Suzuki)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   120   F =  -44.0000        MAXCV =  0.444367E-11
!  X =   0.607440E-06  0.999999       2.00000      -1.00000
!
!    Least squares error in variables =     0.152955E-05
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   150   F =  -44.0000        MAXCV =  0.000000
!  X =   0.193525E-07   1.00000       2.00000      -1.00000
!
!    Least squares error in variables =     0.501895E-07
! ------------------------------------------------------------------
!
!   Output from test problem 9 (Hock and Schittkowski 100)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   422   F =   680.630        MAXCV =  0.595093E-11
!  X =    2.33050       1.95137     -0.477538       4.36573     -0.624488
!    1.03813        1.59423
!
!    Least squares error in variables =     0.416981E-05
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   476   F =   680.630        MAXCV =  0.813238E-14
!  X =    2.33050       1.95137     -0.477541       4.36573     -0.624487
!    1.03813        1.59423
!
!    Least squares error in variables =     0.809366E-06
! ------------------------------------------------------------------
!
!   Output from test problem 10 (Hexagon area)
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   243   F = -0.866025        MAXCV =  0.971976E-23
!  X =   0.687780      0.725919     -0.284774      0.958595      0.687780
!   0.725919      -0.284774
!   0.958595       0.141321E-22
!
!    Least squares error in variables =      1.44346
!
!  Normal return from subroutine COBYLA
!
!  NFVALS =   285   F = -0.866025        MAXCV =  0.220618E-15
!  X =   0.687780      0.725919     -0.284774      0.958595      0.687780
!   0.725919      -0.284774
!   0.958595      -0.814519E-24
!
!    Least squares error in variables =      1.44346
! ------------------------------------------------------------------