scbd/cobyla/email.txt

                                 COBYLA
                                 ~~~~~~

     Here is a single-precision Fortran implementation of the algorithm for
constrained optimization that is the subject of the report I have written on
"A direct search optimization method that models the objective and constraint
functions by linear interpolation". This report has the number DAMTP 1992/NA5,
University of Cambridge, and it has been published in the proceedings of the
conference on Numerical Analysis and Optimization that was held in Oaxaca,
Mexico in January, 1992, which is the book "Advances in Optimization and
Numerical Analysis" (eds. Susana Gomez and Jean-Pierre Hennart), Kluwer
Academic Publishers (1994).

     The instructions for using the Fortran code are given in the comments of
SUBROUTINE COBYLA, which is the interface between the user and the main
calculation that is done by SUBROUTINE COBYLB. There is a need for a linear
programming problem to be solved subject to a Euclidean norm trust region
constraint. Therefore SUBROUTINE TRSTLP is provided too, but you may have some
software that you prefer to use instead. These 3 subroutines are separated by
lines of hyphens below. Further, there follows the main program, the CALCFC
subroutine and the output that are appropriate to the numerical examples that
are discussed in the last section of DAMTP 1992/NA5. Please note, however,
that some cosmetic restructuring of the software has caused the given output
to differ slightly from Table 1 of the report.

     There are no restrictions on the use of the software, nor do I offer any
guarantees of success. Indeed, at the time of writing this note I had applied
it only to test problems that have up to 10 variables.

                        Mike Powell (May 7th, 1992).

-------------------------------------------------------------------------------
      SUBROUTINE COBYLA (N,M,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W,IACT)
      DIMENSION X(*),W(*),IACT(*)
C
C     This subroutine minimizes an objective function F(X) subject to M
C     inequality constraints on X, where X is a vector of variables that has
C     N components. The algorithm employs linear approximations to the
C     objective and constraint functions, the approximations being formed by
C     linear interpolation at N+1 points in the space of the variables.
C     We regard these interpolation points as vertices of a simplex. The
C     parameter RHO controls the size of the simplex and it is reduced
C     automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
C     to achieve a good vector of variables for the current size, and then
C     RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and
C     RHOEND should be set to reasonable initial changes to and the required   
C     accuracy in the variables respectively, but this accuracy should be
C     viewed as a subject for experimentation because it is not guaranteed.
C     The subroutine has an advantage over many of its competitors, however,
C     which is that it treats each constraint individually when calculating
C     a change to the variables, instead of lumping the constraints together
C     into a single penalty function. The name of the subroutine is derived
C     from the phrase Constrained Optimization BY Linear Approximations.
C
C     The user must set the values of N, M, RHOBEG and RHOEND, and must
C     provide an initial vector of variables in X. Further, the value of
C     IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
C     printing during the calculation. Specifically, there is no output if
C     IPRINT=0 and there is output only at the end of the calculation if
C     IPRINT=1. Otherwise each new value of RHO and SIGMA is printed.
C     Further, the vector of variables and some function information are
C     given either when RHO is reduced or when each new value of F(X) is
C     computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
C     is a penalty parameter, it being assumed that a change to X is an
C     improvement if it reduces the merit function
C                F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)),
C     where C1,C2,...,CM denote the constraint functions that should become
C     nonnegative eventually, at least to the precision of RHOEND. In the
C     printed output the displayed term that is multiplied by SIGMA is
C     called MAXCV, which stands for 'MAXimum Constraint Violation'. The
C     argument MAXFUN is an integer variable that must be set by the user to a
C     limit on the number of calls of CALCFC, the purpose of this routine being
C     given below. The value of MAXFUN will be altered to the number of calls
C     of CALCFC that are made. The arguments W and IACT provide real and
C     integer arrays that are used as working space. Their lengths must be at
C     least N*(3*N+2*M+11)+4*M+6 and M+1 respectively.
C
C     In order to define the objective and constraint functions, we require
C     a subroutine that has the name and arguments
C                SUBROUTINE CALCFC (N,M,X,F,CON)
C                DIMENSION X(*),CON(*)  .
C     The values of N and M are fixed and have been defined already, while
C     X is now the current vector of variables. The subroutine should return
C     the objective and constraint functions at X in F and CON(1),CON(2),
C     ...,CON(M). Note that we are trying to adjust X so that F(X) is as
C     small as possible subject to the constraint functions being nonnegative.
C
C     Partition the working space array W to provide the storage that is needed
C     for the main calculation.
C
      MPP=M+2
      ICON=1
      ISIM=ICON+MPP
      ISIMI=ISIM+N*N+N
      IDATM=ISIMI+N*N
      IA=IDATM+N*MPP+MPP
      IVSIG=IA+M*N+N
      IVETA=IVSIG+N
      ISIGB=IVETA+N
      IDX=ISIGB+N
      IWORK=IDX+N
      CALL COBYLB (N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W(ICON),
     1  W(ISIM),W(ISIMI),W(IDATM),W(IA),W(IVSIG),W(IVETA),W(ISIGB),
     2  W(IDX),W(IWORK),IACT)
      RETURN
      END
C------------------------------------------------------------------------------
      SUBROUTINE COBYLB (N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,
     1  CON,SIM,SIMI,DATMAT,A,VSIG,VETA,SIGBAR,DX,W,IACT)
      DIMENSION X(*),CON(*),SIM(N,*),SIMI(N,*),DATMAT(MPP,*),
     1  A(N,*),VSIG(*),VETA(*),SIGBAR(*),DX(*),W(*),IACT(*)
C
C     Set the initial values of some parameters. The last column of SIM holds
C     the optimal vertex of the current simplex, and the preceding N columns
C     hold the displacements from the optimal vertex to the other vertices.
C     Further, SIMI holds the inverse of the matrix that is contained in the
C     first N columns of SIM.
C
      IPTEM=MIN0(N,5)
      IPTEMP=IPTEM+1
      NP=N+1
      MP=M+1
      ALPHA=0.25
      BETA=2.1
      GAMMA=0.5
      DELTA=1.1
      RHO=RHOBEG
      PARMU=0.0
      IF (IPRINT .GE. 2) PRINT 10, RHO
   10 FORMAT (/3X,'The initial value of RHO is',1PE13.6,2X,
     1  'and PARMU is set to zero.')
      NFVALS=0
      TEMP=1.0/RHO
      DO 30 I=1,N
      SIM(I,NP)=X(I)
      DO 20 J=1,N
      SIM(I,J)=0.0
   20 SIMI(I,J)=0.0
      SIM(I,I)=RHO
   30 SIMI(I,I)=TEMP
      JDROP=NP
      IBRNCH=0
C
C     Make the next call of the user-supplied subroutine CALCFC. These
C     instructions are also used for calling CALCFC during the iterations of
C     the algorithm.
C
   40 IF (NFVALS .GE. MAXFUN .AND. NFVALS .GT. 0) THEN
          IF (IPRINT .GE. 1) PRINT 50
   50     FORMAT (/3X,'Return from subroutine COBYLA because the ',
     1      'MAXFUN limit has been reached.')
          GOTO 600
      END IF
      NFVALS=NFVALS+1
      CALL CALCFC (N,M,X,F,CON)
      RESMAX=0.0
      IF (M .GT. 0) THEN
          DO 60 K=1,M
   60     RESMAX=AMAX1(RESMAX,-CON(K))
      END IF
      IF (NFVALS .EQ. IPRINT-1 .OR. IPRINT .EQ. 3) THEN
          PRINT 70, NFVALS,F,RESMAX,(X(I),I=1,IPTEM)
   70     FORMAT (/3X,'NFVALS =',I5,3X,'F =',1PE13.6,4X,'MAXCV =',
     1      1PE13.6/3X,'X =',1PE13.6,1P4E15.6)
          IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
   80     FORMAT (1PE19.6,1P4E15.6)
      END IF
      CON(MP)=F
      CON(MPP)=RESMAX
      IF (IBRNCH .EQ. 1) GOTO 440
C
C     Set the recently calculated function values in a column of DATMAT. This
C     array has a column for each vertex of the current simplex, the entries of
C     each column being the values of the constraint functions (if any)
C     followed by the objective function and the greatest constraint violation
C     at the vertex.
C
      DO 90 K=1,MPP
   90 DATMAT(K,JDROP)=CON(K)
      IF (NFVALS .GT. NP) GOTO 130
C
C     Exchange the new vertex of the initial simplex with the optimal vertex if
C     necessary. Then, if the initial simplex is not complete, pick its next
C     vertex and calculate the function values there.
C
      IF (JDROP .LE. N) THEN
          IF (DATMAT(MP,NP) .LE. F) THEN
              X(JDROP)=SIM(JDROP,NP)
          ELSE
              SIM(JDROP,NP)=X(JDROP)
              DO 100 K=1,MPP
              DATMAT(K,JDROP)=DATMAT(K,NP)
  100         DATMAT(K,NP)=CON(K)
              DO 120 K=1,JDROP
              SIM(JDROP,K)=-RHO
              TEMP=0.0
              DO 110 I=K,JDROP
  110         TEMP=TEMP-SIMI(I,K)
  120         SIMI(JDROP,K)=TEMP
          END IF
      END IF
      IF (NFVALS .LE. N) THEN
          JDROP=NFVALS
          X(JDROP)=X(JDROP)+RHO
          GOTO 40
      END IF
  130 IBRNCH=1
C
C     Identify the optimal vertex of the current simplex.
C
  140 PHIMIN=DATMAT(MP,NP)+PARMU*DATMAT(MPP,NP)
      NBEST=NP
      DO 150 J=1,N
      TEMP=DATMAT(MP,J)+PARMU*DATMAT(MPP,J)
      IF (TEMP .LT. PHIMIN) THEN
          NBEST=J
          PHIMIN=TEMP
      ELSE IF (TEMP .EQ. PHIMIN .AND. PARMU .EQ. 0.0) THEN
          IF (DATMAT(MPP,J) .LT. DATMAT(MPP,NBEST)) NBEST=J
      END IF
  150 CONTINUE
C
C     Switch the best vertex into pole position if it is not there already,
C     and also update SIM, SIMI and DATMAT.
C
      IF (NBEST .LE. N) THEN
          DO 160 I=1,MPP
          TEMP=DATMAT(I,NP)
          DATMAT(I,NP)=DATMAT(I,NBEST)
  160     DATMAT(I,NBEST)=TEMP
          DO 180 I=1,N
          TEMP=SIM(I,NBEST)
          SIM(I,NBEST)=0.0
          SIM(I,NP)=SIM(I,NP)+TEMP
          TEMPA=0.0
          DO 170 K=1,N
          SIM(I,K)=SIM(I,K)-TEMP
  170     TEMPA=TEMPA-SIMI(K,I)
  180     SIMI(NBEST,I)=TEMPA
      END IF
C
C     Make an error return if SIGI is a poor approximation to the inverse of
C     the leading N by N submatrix of SIG.
C
      ERROR=0.0
      DO 200 I=1,N
      DO 200 J=1,N
      TEMP=0.0
      IF (I .EQ. J) TEMP=TEMP-1.0
      DO 190 K=1,N
  190 TEMP=TEMP+SIMI(I,K)*SIM(K,J)
  200 ERROR=AMAX1(ERROR,ABS(TEMP))
      IF (ERROR .GT. 0.1) THEN
          IF (IPRINT .GE. 1) PRINT 210
  210     FORMAT (/3X,'Return from subroutine COBYLA because ',
     1      'rounding errors are becoming damaging.')
          GOTO 600
      END IF
C
C     Calculate the coefficients of the linear approximations to the objective
C     and constraint functions, placing minus the objective function gradient
C     after the constraint gradients in the array A. The vector W is used for
C     working space.
C
      DO 240 K=1,MP
      CON(K)=-DATMAT(K,NP)
      DO 220 J=1,N
  220 W(J)=DATMAT(K,J)+CON(K)
      DO 240 I=1,N
      TEMP=0.0
      DO 230 J=1,N
  230 TEMP=TEMP+W(J)*SIMI(J,I)
      IF (K .EQ. MP) TEMP=-TEMP
  240 A(I,K)=TEMP
C
C     Calculate the values of sigma and eta, and set IFLAG=0 if the current
C     simplex is not acceptable.
C
      IFLAG=1
      PARSIG=ALPHA*RHO
      PARETA=BETA*RHO
      DO 260 J=1,N
      WSIG=0.0
      WETA=0.0
      DO 250 I=1,N
      WSIG=WSIG+SIMI(J,I)**2
  250 WETA=WETA+SIM(I,J)**2
      VSIG(J)=1.0/SQRT(WSIG)
      VETA(J)=SQRT(WETA)
      IF (VSIG(J) .LT. PARSIG .OR. VETA(J) .GT. PARETA) IFLAG=0
  260 CONTINUE
C
C     If a new vertex is needed to improve acceptability, then decide which
C     vertex to drop from the simplex.
C
      IF (IBRNCH .EQ. 1 .OR. IFLAG .EQ. 1) GOTO 370
      JDROP=0
      TEMP=PARETA
      DO 270 J=1,N
      IF (VETA(J) .GT. TEMP) THEN
          JDROP=J
          TEMP=VETA(J)
      END IF
  270 CONTINUE
      IF (JDROP .EQ. 0) THEN
          DO 280 J=1,N
          IF (VSIG(J) .LT. TEMP) THEN
              JDROP=J
              TEMP=VSIG(J)
          END IF
  280     CONTINUE
      END IF
C
C     Calculate the step to the new vertex and its sign.
C
      TEMP=GAMMA*RHO*VSIG(JDROP)
      DO 290 I=1,N
  290 DX(I)=TEMP*SIMI(JDROP,I)
      CVMAXP=0.0
      CVMAXM=0.0
      DO 310 K=1,MP
      SUM=0.0
      DO 300 I=1,N
  300 SUM=SUM+A(I,K)*DX(I)
      IF (K .LT. MP) THEN
          TEMP=DATMAT(K,NP)
          CVMAXP=AMAX1(CVMAXP,-SUM-TEMP)
          CVMAXM=AMAX1(CVMAXM,SUM-TEMP)
      END IF
  310 CONTINUE
      DXSIGN=1.0
      IF (PARMU*(CVMAXP-CVMAXM) .GT. SUM+SUM) DXSIGN=-1.0
C
C     Update the elements of SIM and SIMI, and set the next X.
C
      TEMP=0.0
      DO 320 I=1,N
      DX(I)=DXSIGN*DX(I)
      SIM(I,JDROP)=DX(I)
  320 TEMP=TEMP+SIMI(JDROP,I)*DX(I)
      DO 330 I=1,N
  330 SIMI(JDROP,I)=SIMI(JDROP,I)/TEMP
      DO 360 J=1,N
      IF (J .NE. JDROP) THEN
          TEMP=0.0
          DO 340 I=1,N
  340     TEMP=TEMP+SIMI(J,I)*DX(I)
          DO 350 I=1,N
  350     SIMI(J,I)=SIMI(J,I)-TEMP*SIMI(JDROP,I)
      END IF
  360 X(J)=SIM(J,NP)+DX(J)
      GOTO 40
C
C     Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO.
C
  370 IZ=1
      IZDOTA=IZ+N*N
      IVMC=IZDOTA+N
      ISDIRN=IVMC+MP
      IDXNEW=ISDIRN+N
      IVMD=IDXNEW+N
      CALL TRSTLP (N,M,A,CON,RHO,DX,IFULL,IACT,W(IZ),W(IZDOTA),
     1  W(IVMC),W(ISDIRN),W(IDXNEW),W(IVMD))
      IF (IFULL .EQ. 0) THEN
          TEMP=0.0
          DO 380 I=1,N
  380     TEMP=TEMP+DX(I)**2
          IF (TEMP .LT. 0.25*RHO*RHO) THEN
              IBRNCH=1
              GOTO 550
          END IF
      END IF
C
C     Predict the change to F and the new maximum constraint violation if the
C     variables are altered from x(0) to x(0)+DX.
C
      RESNEW=0.0
      CON(MP)=0.0
      DO 400 K=1,MP
      SUM=CON(K)
      DO 390 I=1,N
  390 SUM=SUM-A(I,K)*DX(I)
      IF (K .LT. MP) RESNEW=AMAX1(RESNEW,SUM)
  400 CONTINUE
C
C     Increase PARMU if necessary and branch back if this change alters the
C     optimal vertex. Otherwise PREREM and PREREC will be set to the predicted
C     reductions in the merit function and the maximum constraint violation
C     respectively.
C
      BARMU=0.0
      PREREC=DATMAT(MPP,NP)-RESNEW
      IF (PREREC .GT. 0.0) BARMU=SUM/PREREC
      IF (PARMU .LT. 1.5*BARMU) THEN
          PARMU=2.0*BARMU
          IF (IPRINT .GE. 2) PRINT 410, PARMU
  410     FORMAT (/3X,'Increase in PARMU to',1PE13.6)
          PHI=DATMAT(MP,NP)+PARMU*DATMAT(MPP,NP)
          DO 420 J=1,N
          TEMP=DATMAT(MP,J)+PARMU*DATMAT(MPP,J)
          IF (TEMP .LT. PHI) GOTO 140
          IF (TEMP .EQ. PHI .AND. PARMU .EQ. 0.0) THEN
              IF (DATMAT(MPP,J) .LT. DATMAT(MPP,NP)) GOTO 140
          END IF
  420     CONTINUE
      END IF
      PREREM=PARMU*PREREC-SUM
C
C     Calculate the constraint and objective functions at x(*). Then find the
C     actual reduction in the merit function.
C
      DO 430 I=1,N
  430 X(I)=SIM(I,NP)+DX(I)
      IBRNCH=1
      GOTO 40
  440 VMOLD=DATMAT(MP,NP)+PARMU*DATMAT(MPP,NP)
      VMNEW=F+PARMU*RESMAX
      TRURED=VMOLD-VMNEW
      IF (PARMU .EQ. 0.0 .AND. F .EQ. DATMAT(MP,NP)) THEN
          PREREM=PREREC
          TRURED=DATMAT(MPP,NP)-RESMAX
      END IF
C
C     Begin the operations that decide whether x(*) should replace one of the
C     vertices of the current simplex, the change being mandatory if TRURED is
C     positive. Firstly, JDROP is set to the index of the vertex that is to be
C     replaced.
C
      RATIO=0.0
      IF (TRURED .LE. 0.0) RATIO=1.0
      JDROP=0
      DO 460 J=1,N
      TEMP=0.0
      DO 450 I=1,N
  450 TEMP=TEMP+SIMI(J,I)*DX(I)
      TEMP=ABS(TEMP)
      IF (TEMP .GT. RATIO) THEN
          JDROP=J
          RATIO=TEMP
      END IF
  460 SIGBAR(J)=TEMP*VSIG(J)
C
C     Calculate the value of ell.
C
      EDGMAX=DELTA*RHO
      L=0
      DO 480 J=1,N
      IF (SIGBAR(J) .GE. PARSIG .OR. SIGBAR(J) .GE. VSIG(J)) THEN
          TEMP=VETA(J)
          IF (TRURED .GT. 0.0) THEN
              TEMP=0.0
              DO 470 I=1,N
  470         TEMP=TEMP+(DX(I)-SIM(I,J))**2
              TEMP=SQRT(TEMP)
          END IF
          IF (TEMP .GT. EDGMAX) THEN
              L=J
              EDGMAX=TEMP
          END IF
      END IF
  480 CONTINUE
      IF (L .GT. 0) JDROP=L
      IF (JDROP .EQ. 0) GOTO 550
C
C     Revise the simplex by updating the elements of SIM, SIMI and DATMAT.
C
      TEMP=0.0
      DO 490 I=1,N
      SIM(I,JDROP)=DX(I)
  490 TEMP=TEMP+SIMI(JDROP,I)*DX(I)
      DO 500 I=1,N
  500 SIMI(JDROP,I)=SIMI(JDROP,I)/TEMP
      DO 530 J=1,N
      IF (J .NE. JDROP) THEN
          TEMP=0.0
          DO 510 I=1,N
  510     TEMP=TEMP+SIMI(J,I)*DX(I)
          DO 520 I=1,N
  520     SIMI(J,I)=SIMI(J,I)-TEMP*SIMI(JDROP,I)
      END IF
  530 CONTINUE
      DO 540 K=1,MPP
  540 DATMAT(K,JDROP)=CON(K)
C
C     Branch back for further iterations with the current RHO.
C
      IF (TRURED .GT. 0.0 .AND. TRURED .GE. 0.1*PREREM) GOTO 140
  550 IF (IFLAG .EQ. 0) THEN
          IBRNCH=0
          GOTO 140
      END IF
C
C     Otherwise reduce RHO if it is not at its least value and reset PARMU.
C
      IF (RHO .GT. RHOEND) THEN
          RHO=0.5*RHO
          IF (RHO .LE. 1.5*RHOEND) RHO=RHOEND
          IF (PARMU .GT. 0.0) THEN
              DENOM=0.0
              DO 570 K=1,MP
              CMIN=DATMAT(K,NP)
              CMAX=CMIN
              DO 560 I=1,N
              CMIN=AMIN1(CMIN,DATMAT(K,I))
  560         CMAX=AMAX1(CMAX,DATMAT(K,I))
              IF (K .LE. M .AND. CMIN .LT. 0.5*CMAX) THEN
                  TEMP=AMAX1(CMAX,0.0)-CMIN
                  IF (DENOM .LE. 0.0) THEN
                      DENOM=TEMP
                  ELSE
                      DENOM=AMIN1(DENOM,TEMP)
                  END IF
              END IF
  570         CONTINUE
              IF (DENOM .EQ. 0.0) THEN
                  PARMU=0.0
              ELSE IF (CMAX-CMIN .LT. PARMU*DENOM) THEN
                  PARMU=(CMAX-CMIN)/DENOM
              END IF
          END IF
          IF (IPRINT .GE. 2) PRINT 580, RHO,PARMU
  580     FORMAT (/3X,'Reduction in RHO to',1PE13.6,'  and PARMU =',
     1      1PE13.6)
          IF (IPRINT .EQ. 2) THEN
              PRINT 70, NFVALS,DATMAT(MP,NP),DATMAT(MPP,NP),
     1          (SIM(I,NP),I=1,IPTEM)
              IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
          END IF
          GOTO 140
      END IF
C
C     Return the best calculated values of the variables.
C
      IF (IPRINT .GE. 1) PRINT 590
  590 FORMAT (/3X,'Normal return from subroutine COBYLA')
      IF (IFULL .EQ. 1) GOTO 620
  600 DO 610 I=1,N
  610 X(I)=SIM(I,NP)
      F=DATMAT(MP,NP)
      RESMAX=DATMAT(MPP,NP)
  620 IF (IPRINT .GE. 1) THEN
          PRINT 70, NFVALS,F,RESMAX,(X(I),I=1,IPTEM)
          IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
      END IF
      MAXFUN=NFVALS
      RETURN
      END
C------------------------------------------------------------------------------
      SUBROUTINE TRSTLP (N,M,A,B,RHO,DX,IFULL,IACT,Z,ZDOTA,VMULTC,
     1  SDIRN,DXNEW,VMULTD) 
      DIMENSION A(N,*),B(*),DX(*),IACT(*),Z(N,*),ZDOTA(*),
     1  VMULTC(*),SDIRN(*),DXNEW(*),VMULTD(*)
C
C     This subroutine calculates an N-component vector DX by applying the
C     following two stages. In the first stage, DX is set to the shortest
C     vector that minimizes the greatest violation of the constraints
C       A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K=2,3,...,M,
C     subject to the Euclidean length of DX being at most RHO. If its length is
C     strictly less than RHO, then we use the resultant freedom in DX to
C     minimize the objective function
C              -A(1,M+1)*DX(1)-A(2,M+1)*DX(2)-...-A(N,M+1)*DX(N)
C     subject to no increase in any greatest constraint violation. This
C     notation allows the gradient of the objective function to be regarded as
C     the gradient of a constraint. Therefore the two stages are distinguished
C     by MCON .EQ. M and MCON .GT. M respectively. It is possible that a
C     degeneracy may prevent DX from attaining the target length RHO. Then the
C     value IFULL=0 would be set, but usually IFULL=1 on return.
C
C     In general NACT is the number of constraints in the active set and
C     IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT
C     contains a permutation of the remaining constraint indices. Further, Z is
C     an orthogonal matrix whose first NACT columns can be regarded as the
C     result of Gram-Schmidt applied to the active constraint gradients. For
C     J=1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th
C     column of Z with the gradient of the J-th active constraint. DX is the
C     current vector of variables and here the residuals of the active
C     constraints should be zero. Further, the active constraints have
C     nonnegative Lagrange multipliers that are held at the beginning of
C     VMULTC. The remainder of this vector holds the residuals of the inactive
C     constraints at DX, the ordering of the components of VMULTC being in
C     agreement with the permutation of the indices of the constraints that is
C     in IACT. All these residuals are nonnegative, which is achieved by the
C     shift RESMAX that makes the least residual zero.
C
C     Initialize Z and some other variables. The value of RESMAX will be
C     appropriate to DX=0, while ICON will be the index of a most violated
C     constraint if RESMAX is positive. Usually during the first stage the
C     vector SDIRN gives a search direction that reduces all the active
C     constraint violations by one simultaneously.
C
      IFULL=1
      MCON=M
      NACT=0
      RESMAX=0.0
      DO 20 I=1,N
      DO 10 J=1,N
   10 Z(I,J)=0.0
      Z(I,I)=1.0
   20 DX(I)=0.0
      IF (M .GE. 1) THEN
          DO 30 K=1,M
          IF (B(K) .GT. RESMAX) THEN
              RESMAX=B(K)
              ICON=K
          END IF
   30     CONTINUE
          DO 40 K=1,M
          IACT(K)=K
   40     VMULTC(K)=RESMAX-B(K)
      END IF
      IF (RESMAX .EQ. 0.0) GOTO 480
      DO 50 I=1,N
   50 SDIRN(I)=0.0
C
C     End the current stage of the calculation if 3 consecutive iterations
C     have either failed to reduce the best calculated value of the objective
C     function or to increase the number of active constraints since the best
C     value was calculated. This strategy prevents cycling, but there is a
C     remote possibility that it will cause premature termination.
C
   60 OPTOLD=0.0
      ICOUNT=0
   70 IF (MCON .EQ. M) THEN
          OPTNEW=RESMAX
      ELSE
          OPTNEW=0.0
          DO 80 I=1,N
   80     OPTNEW=OPTNEW-DX(I)*A(I,MCON)
      END IF
      IF (ICOUNT .EQ. 0 .OR. OPTNEW .LT. OPTOLD) THEN
          OPTOLD=OPTNEW
          NACTX=NACT
          ICOUNT=3
      ELSE IF (NACT .GT. NACTX) THEN
          NACTX=NACT
          ICOUNT=3
      ELSE
          ICOUNT=ICOUNT-1
          IF (ICOUNT .EQ. 0) GOTO 490
      END IF
C
C     If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to
C     the active set. Apply Givens rotations so that the last N-NACT-1 columns
C     of Z are orthogonal to the gradient of the new constraint, a scalar
C     product being set to zero if its nonzero value could be due to computer
C     rounding errors. The array DXNEW is used for working space.
C
      IF (ICON .LE. NACT) GOTO 260
      KK=IACT(ICON)
      DO 90 I=1,N
   90 DXNEW(I)=A(I,KK)
      TOT=0.0
      K=N
  100 IF (K .GT. NACT) THEN
          SP=0.0
          SPABS=0.0
          DO 110 I=1,N
          TEMP=Z(I,K)*DXNEW(I)
          SP=SP+TEMP
  110     SPABS=SPABS+ABS(TEMP)
          ACCA=SPABS+0.1*ABS(SP)
          ACCB=SPABS+0.2*ABS(SP)
          IF (SPABS .GE. ACCA .OR. ACCA .GE. ACCB) SP=0.0
          IF (TOT .EQ. 0.0) THEN
              TOT=SP
          ELSE
              KP=K+1
              TEMP=SQRT(SP*SP+TOT*TOT)
              ALPHA=SP/TEMP
              BETA=TOT/TEMP
              TOT=TEMP
              DO 120 I=1,N
              TEMP=ALPHA*Z(I,K)+BETA*Z(I,KP)
              Z(I,KP)=ALPHA*Z(I,KP)-BETA*Z(I,K)
  120         Z(I,K)=TEMP
          END IF
          K=K-1
          GOTO 100
      END IF
C
C     Add the new constraint if this can be done without a deletion from the
C     active set.
C
      IF (TOT .NE. 0.0) THEN
          NACT=NACT+1
          ZDOTA(NACT)=TOT
          VMULTC(ICON)=VMULTC(NACT)
          VMULTC(NACT)=0.0
          GOTO 210
      END IF
C
C     The next instruction is reached if a deletion has to be made from the
C     active set in order to make room for the new active constraint, because
C     the new constraint gradient is a linear combination of the gradients of
C     the old active constraints. Set the elements of VMULTD to the multipliers
C     of the linear combination. Further, set IOUT to the index of the
C     constraint to be deleted, but branch if no suitable index can be found.
C
      RATIO=-1.0
      K=NACT
  130 ZDOTV=0.0
      ZDVABS=0.0
      DO 140 I=1,N
      TEMP=Z(I,K)*DXNEW(I)
      ZDOTV=ZDOTV+TEMP
  140 ZDVABS=ZDVABS+ABS(TEMP)
      ACCA=ZDVABS+0.1*ABS(ZDOTV)
      ACCB=ZDVABS+0.2*ABS(ZDOTV)
      IF (ZDVABS .LT. ACCA .AND. ACCA .LT. ACCB) THEN
          TEMP=ZDOTV/ZDOTA(K)
          IF (TEMP .GT. 0.0 .AND. IACT(K) .LE. M) THEN
              TEMPA=VMULTC(K)/TEMP
              IF (RATIO .LT. 0.0 .OR. TEMPA .LT. RATIO) THEN
                  RATIO=TEMPA
                  IOUT=K
              END IF
           END IF
          IF (K .GE. 2) THEN
              KW=IACT(K)
              DO 150 I=1,N
  150         DXNEW(I)=DXNEW(I)-TEMP*A(I,KW)
          END IF
          VMULTD(K)=TEMP
      ELSE
          VMULTD(K)=0.0
      END IF
      K=K-1
      IF (K .GT. 0) GOTO 130
      IF (RATIO .LT. 0.0) GOTO 490
C
C     Revise the Lagrange multipliers and reorder the active constraints so
C     that the one to be replaced is at the end of the list. Also calculate the
C     new value of ZDOTA(NACT) and branch if it is not acceptable.
C
      DO 160 K=1,NACT
  160 VMULTC(K)=AMAX1(0.0,VMULTC(K)-RATIO*VMULTD(K))
      IF (ICON .LT. NACT) THEN
          ISAVE=IACT(ICON)
          VSAVE=VMULTC(ICON)
          K=ICON
  170     KP=K+1
          KW=IACT(KP)
          SP=0.0
          DO 180 I=1,N
  180     SP=SP+Z(I,K)*A(I,KW)
          TEMP=SQRT(SP*SP+ZDOTA(KP)**2)
          ALPHA=ZDOTA(KP)/TEMP
          BETA=SP/TEMP
          ZDOTA(KP)=ALPHA*ZDOTA(K)
          ZDOTA(K)=TEMP
          DO 190 I=1,N
          TEMP=ALPHA*Z(I,KP)+BETA*Z(I,K)
          Z(I,KP)=ALPHA*Z(I,K)-BETA*Z(I,KP)
  190     Z(I,K)=TEMP
          IACT(K)=KW
          VMULTC(K)=VMULTC(KP)
          K=KP
          IF (K .LT. NACT) GOTO 170
          IACT(K)=ISAVE
          VMULTC(K)=VSAVE
      END IF
      TEMP=0.0
      DO 200 I=1,N
  200 TEMP=TEMP+Z(I,NACT)*A(I,KK)
      IF (TEMP .EQ. 0.0) GOTO 490
      ZDOTA(NACT)=TEMP
      VMULTC(ICON)=0.0
      VMULTC(NACT)=RATIO
C
C     Update IACT and ensure that the objective function continues to be
C     treated as the last active constraint when MCON>M.
C
  210 IACT(ICON)=IACT(NACT)
      IACT(NACT)=KK
      IF (MCON .GT. M .AND. KK .NE. MCON) THEN
          K=NACT-1
          SP=0.0
          DO 220 I=1,N
  220     SP=SP+Z(I,K)*A(I,KK)
          TEMP=SQRT(SP*SP+ZDOTA(NACT)**2)
          ALPHA=ZDOTA(NACT)/TEMP
          BETA=SP/TEMP
          ZDOTA(NACT)=ALPHA*ZDOTA(K)
          ZDOTA(K)=TEMP
          DO 230 I=1,N
          TEMP=ALPHA*Z(I,NACT)+BETA*Z(I,K)
          Z(I,NACT)=ALPHA*Z(I,K)-BETA*Z(I,NACT)
  230     Z(I,K)=TEMP
          IACT(NACT)=IACT(K)
          IACT(K)=KK
          TEMP=VMULTC(K)
          VMULTC(K)=VMULTC(NACT)
          VMULTC(NACT)=TEMP
      END IF
C
C     If stage one is in progress, then set SDIRN to the direction of the next
C     change to the current vector of variables.
C
      IF (MCON .GT. M) GOTO 320
      KK=IACT(NACT)
      TEMP=0.0
      DO 240 I=1,N
  240 TEMP=TEMP+SDIRN(I)*A(I,KK)
      TEMP=TEMP-1.0
      TEMP=TEMP/ZDOTA(NACT)
      DO 250 I=1,N
  250 SDIRN(I)=SDIRN(I)-TEMP*Z(I,NACT)
      GOTO 340
C
C     Delete the constraint that has the index IACT(ICON) from the active set.
C
  260 IF (ICON .LT. NACT) THEN
          ISAVE=IACT(ICON)
          VSAVE=VMULTC(ICON)
          K=ICON
  270     KP=K+1
          KK=IACT(KP)
          SP=0.0
          DO 280 I=1,N
  280     SP=SP+Z(I,K)*A(I,KK)
          TEMP=SQRT(SP*SP+ZDOTA(KP)**2)
          ALPHA=ZDOTA(KP)/TEMP
          BETA=SP/TEMP
          ZDOTA(KP)=ALPHA*ZDOTA(K)
          ZDOTA(K)=TEMP
          DO 290 I=1,N
          TEMP=ALPHA*Z(I,KP)+BETA*Z(I,K)
          Z(I,KP)=ALPHA*Z(I,K)-BETA*Z(I,KP)
  290     Z(I,K)=TEMP
          IACT(K)=KK
          VMULTC(K)=VMULTC(KP)
          K=KP
          IF (K .LT. NACT) GOTO 270
          IACT(K)=ISAVE
          VMULTC(K)=VSAVE
      END IF
      NACT=NACT-1
C
C     If stage one is in progress, then set SDIRN to the direction of the next
C     change to the current vector of variables.
C
      IF (MCON .GT. M) GOTO 320
      TEMP=0.0
      DO 300 I=1,N
  300 TEMP=TEMP+SDIRN(I)*Z(I,NACT+1)
      DO 310 I=1,N
  310 SDIRN(I)=SDIRN(I)-TEMP*Z(I,NACT+1)
      GO TO 340
C
C     Pick the next search direction of stage two.
C
  320 TEMP=1.0/ZDOTA(NACT)
      DO 330 I=1,N
  330 SDIRN(I)=TEMP*Z(I,NACT)
C
C     Calculate the step to the boundary of the trust region or take the step
C     that reduces RESMAX to zero. The two statements below that include the
C     factor 1.0E-6 prevent some harmless underflows that occurred in a test
C     calculation. Further, we skip the step if it could be zero within a
C     reasonable tolerance for computer rounding errors.
C
  340 DD=RHO*RHO
      SD=0.0
      SS=0.0
      DO 350 I=1,N
      IF (ABS(DX(I)) .GE. 1.0E-6*RHO) DD=DD-DX(I)**2
      SD=SD+DX(I)*SDIRN(I)
  350 SS=SS+SDIRN(I)**2
      IF (DD .LE. 0.0) GOTO 490
      TEMP=SQRT(SS*DD)
      IF (ABS(SD) .GE. 1.0E-6*TEMP) TEMP=SQRT(SS*DD+SD*SD)
      STPFUL=DD/(TEMP+SD)
      STEP=STPFUL
      IF (MCON .EQ. M) THEN
          ACCA=STEP+0.1*RESMAX
          ACCB=STEP+0.2*RESMAX
          IF (STEP .GE. ACCA .OR. ACCA .GE. ACCB) GOTO 480
          STEP=AMIN1(STEP,RESMAX)
      END IF
C
C     Set DXNEW to the new variables if STEP is the steplength, and reduce
C     RESMAX to the corresponding maximum residual if stage one is being done.
C     Because DXNEW will be changed during the calculation of some Lagrange
C     multipliers, it will be restored to the following value later.
C
      DO 360 I=1,N
  360 DXNEW(I)=DX(I)+STEP*SDIRN(I)
      IF (MCON .EQ. M) THEN
          RESOLD=RESMAX
          RESMAX=0.0
          DO 380 K=1,NACT
          KK=IACT(K)
          TEMP=B(KK)
          DO 370 I=1,N
  370     TEMP=TEMP-A(I,KK)*DXNEW(I)
          RESMAX=AMAX1(RESMAX,TEMP)
  380     CONTINUE
      END IF
C
C     Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A
C     device is included to force VMULTD(K)=0.0 if deviations from this value
C     can be attributed to computer rounding errors. First calculate the new
C     Lagrange multipliers.
C
      K=NACT
  390 ZDOTW=0.0
      ZDWABS=0.0
      DO 400 I=1,N
      TEMP=Z(I,K)*DXNEW(I)
      ZDOTW=ZDOTW+TEMP
  400 ZDWABS=ZDWABS+ABS(TEMP)
      ACCA=ZDWABS+0.1*ABS(ZDOTW)
      ACCB=ZDWABS+0.2*ABS(ZDOTW)
      IF (ZDWABS .GE. ACCA .OR. ACCA .GE. ACCB) ZDOTW=0.0
      VMULTD(K)=ZDOTW/ZDOTA(K)
      IF (K .GE. 2) THEN
          KK=IACT(K)
          DO 410 I=1,N
  410     DXNEW(I)=DXNEW(I)-VMULTD(K)*A(I,KK)
          K=K-1
          GOTO 390
      END IF
      IF (MCON .GT. M) VMULTD(NACT)=AMAX1(0.0,VMULTD(NACT))
C
C     Complete VMULTC by finding the new constraint residuals.
C
      DO 420 I=1,N
  420 DXNEW(I)=DX(I)+STEP*SDIRN(I)
      IF (MCON .GT. NACT) THEN
          KL=NACT+1
          DO 440 K=KL,MCON
          KK=IACT(K)
          SUM=RESMAX-B(KK)
          SUMABS=RESMAX+ABS(B(KK))
          DO 430 I=1,N
          TEMP=A(I,KK)*DXNEW(I)
          SUM=SUM+TEMP
  430     SUMABS=SUMABS+ABS(TEMP)
          ACCA=SUMABS+0.1*ABS(SUM)
          ACCB=SUMABS+0.2*ABS(SUM)
          IF (SUMABS .GE. ACCA .OR. ACCA .GE. ACCB) SUM=0.0
  440     VMULTD(K)=SUM
      END IF
C
C     Calculate the fraction of the step from DX to DXNEW that will be taken.
C
      RATIO=1.0
      ICON=0
      DO 450 K=1,MCON
      IF (VMULTD(K) .LT. 0.0) THEN
          TEMP=VMULTC(K)/(VMULTC(K)-VMULTD(K))
          IF (TEMP .LT. RATIO) THEN
              RATIO=TEMP
              ICON=K
          END IF
      END IF
  450 CONTINUE
C
C     Update DX, VMULTC and RESMAX.
C
      TEMP=1.0-RATIO
      DO 460 I=1,N
  460 DX(I)=TEMP*DX(I)+RATIO*DXNEW(I)
      DO 470 K=1,MCON
  470 VMULTC(K)=AMAX1(0.0,TEMP*VMULTC(K)+RATIO*VMULTD(K))
      IF (MCON .EQ. M) RESMAX=RESOLD+RATIO*(RESMAX-RESOLD)
C
C     If the full step is not acceptable then begin another iteration.
C     Otherwise switch to stage two or end the calculation.
C
      IF (ICON .GT. 0) GOTO 70
      IF (STEP .EQ. STPFUL) GOTO 500
  480 MCON=M+1
      ICON=MCON
      IACT(MCON)=MCON
      VMULTC(MCON)=0.0
      GOTO 60
C
C     We employ any freedom that may be available to reduce the objective
C     function before returning a DX whose length is less than RHO.
C
  490 IF (MCON .EQ. M) GOTO 480
      IFULL=0
  500 RETURN
      END
C------------------------------------------------------------------------------
C     Main program of test problems in Report DAMTP 1992/NA5.
C------------------------------------------------------------------------------
      COMMON NPROB
      DIMENSION X(10),XOPT(10),W(3000),IACT(51)
      DO 180 NPROB=1,10
      IF (NPROB .EQ. 1) THEN
C
C     Minimization of a simple quadratic function of two variables.
C
          PRINT 10
   10     FORMAT (/7X,'Output from test problem 1 (Simple quadratic)')
          N=2
          M=0
          XOPT(1)=-1.0
          XOPT(2)=0.0
      ELSE IF (NPROB .EQ. 2) THEN
C
C     Easy two dimensional minimization in unit circle.
C
          PRINT 20
   20     FORMAT (/7X,'Output from test problem 2 (2D unit circle ',
     1      'calculation)')
          N=2
          M=1
          XOPT(1)=SQRT(0.5)
          XOPT(2)=-XOPT(1)
      ELSE IF (NPROB .EQ. 3) THEN
C
C     Easy three dimensional minimization in ellipsoid.
C
          PRINT 30
   30     FORMAT (/7X,'Output from test problem 3 (3D ellipsoid ',
     1      'calculation)')
          N=3
          M=1
          XOPT(1)=1.0/SQRT(3.0)
          XOPT(2)=1.0/SQRT(6.0)
          XOPT(3)=-1.0/3.0
      ELSE IF (NPROB .EQ. 4) THEN
C
C     Weak version of Rosenbrock's problem.
C
          PRINT 40
   40     FORMAT (/7X,'Output from test problem 4 (Weak Rosenbrock)')
          N=2
          M=0
          XOPT(1)=-1.0
          XOPT(2)=1.0
      ELSE IF (NPROB .EQ. 5) THEN
C
C     Intermediate version of Rosenbrock's problem.
C
          PRINT 50
   50     FORMAT (/7X,'Output from test problem 5 (Intermediate ',
     1      'Rosenbrock)')
          N=2
          M=0
          XOPT(1)=-1.0
          XOPT(2)=1.0
      ELSE IF (NPROB .EQ. 6) THEN
C
C     This problem is taken from Fletcher's book Practical Methods of
C     Optimization and has the equation number (9.1.15).
C
          PRINT 60
   60     FORMAT (/7X,'Output from test problem 6 (Equation ',
     1      '(9.1.15) in Fletcher)')
          N=2
          M=2
          XOPT(1)=SQRT(0.5)
          XOPT(2)=XOPT(1)
      ELSE IF (NPROB .EQ. 7) THEN
C
C     This problem is taken from Fletcher's book Practical Methods of
C     Optimization and has the equation number (14.4.2).
C
          PRINT 70
   70     FORMAT (/7X,'Output from test problem 7 (Equation ',
     1      '(14.4.2) in Fletcher)')
          N=3
          M=3
          XOPT(1)=0.0
          XOPT(2)=-3.0
          XOPT(3)=-3.0
      ELSE IF (NPROB .EQ. 8) THEN
C
C     This problem is taken from page 66 of Hock and Schittkowski's book Test
C     Examples for Nonlinear Programming Codes. It is their test problem Number
C     43, and has the name Rosen-Suzuki.
C
          PRINT 80
   80     FORMAT (/7X,'Output from test problem 8 (Rosen-Suzuki)')
          N=4
          M=3
          XOPT(1)=0.0
          XOPT(2)=1.0
          XOPT(3)=2.0
          XOPT(4)=-1.0
      ELSE IF (NPROB .EQ. 9) THEN
C
C     This problem is taken from page 111 of Hock and Schittkowski's
C     book Test Examples for Nonlinear Programming Codes. It is their
C     test problem Number 100.
C
          PRINT 90
   90     FORMAT (/7X,'Output from test problem 9 (Hock and ',
     1      'Schittkowski 100)')
          N=7
          M=4
          XOPT(1)=2.330499
          XOPT(2)=1.951372
          XOPT(3)=-0.4775414
          XOPT(4)=4.365726
          XOPT(5)=-0.624487
          XOPT(6)=1.038131
          XOPT(7)=1.594227
      ELSE IF (NPROB .EQ. 10) THEN
C
C     This problem is taken from page 415 of Luenberger's book Applied
C     Nonlinear Programming. It is to maximize the area of a hexagon of
C     unit diameter.
C
          PRINT 100
  100     FORMAT (/7X,'Output from test problem 10 (Hexagon area)')
          N=9
          M=14
      END IF
      DO 160 ICASE=1,2
      DO 120 I=1,N
  120 X(I)=1.0
      RHOBEG=0.5
      RHOEND=0.001
      IF (ICASE .EQ. 2) RHOEND=0.0001
      IPRINT=1
      MAXFUN=2000
      CALL COBYLA (N,M,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W,IACT)
      IF (NPROB .EQ. 10) THEN
          TEMPA=X(1)+X(3)+X(5)+X(7)
          TEMPB=X(2)+X(4)+X(6)+X(8)
          TEMPC=0.5/SQRT(TEMPA*TEMPA+TEMPB*TEMPB)
          TEMPD=TEMPC*SQRT(3.0)
          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 130 I=1,4
  130     XOPT(I+4)=XOPT(I)
      END IF
      TEMP=0.0
      DO 140 I=1,N
  140 TEMP=TEMP+(X(I)-XOPT(I))**2
      PRINT 150, SQRT(TEMP)
  150 FORMAT (/5X,'Least squares error in variables =',1PE16.6)
  160 CONTINUE
      PRINT 170
  170 FORMAT (2X,'----------------------------------------------',
     1  '--------------------')
  180 CONTINUE
      STOP
      END
C------------------------------------------------------------------------------
      SUBROUTINE CALCFC (N,M,X,F,CON)
      COMMON NPROB
      DIMENSION X(*),CON(*)
      IF (NPROB .EQ. 1) THEN
C
C     Test problem 1 (Simple quadratic)
C     
          F=10.0*(X(1)+1.0)**2+X(2)**2
      ELSE IF (NPROB .EQ. 2) THEN
C
C    Test problem 2 (2D unit circle calculation)
C
          F=X(1)*X(2)
          CON(1)=1.0-X(1)**2-X(2)**2
      ELSE IF (NPROB .EQ. 3) THEN
C
C     Test problem 3 (3D ellipsoid calculation)
C
          F=X(1)*X(2)*X(3)
          CON(1)=1.0-X(1)**2-2.0*X(2)**2-3.0*X(3)**2
      ELSE IF (NPROB .EQ. 4) THEN
C
C     Test problem 4 (Weak Rosenbrock)
C
          F=(X(1)**2-X(2))**2+(1.0+X(1))**2
      ELSE IF (NPROB .EQ. 5) THEN
C
C     Test problem 5 (Intermediate Rosenbrock)
C
          F=10.0*(X(1)**2-X(2))**2+(1.0+X(1))**2
      ELSE IF (NPROB .EQ. 6) THEN
C
C     Test problem 6 (Equation (9.1.15) in Fletcher's book)
C
          F=-X(1)-X(2)
          CON(1)=X(2)-X(1)**2
          CON(2)=1.0-X(1)**2-X(2)**2
      ELSE IF (NPROB .EQ. 7) THEN
C
C     Test problem 7 (Equation (14.4.2) in Fletcher's book)
C
          F=X(3)
          CON(1)=5.0*X(1)-X(2)+X(3)
          CON(2)=X(3)-X(1)**2-X(2)**2-4.0*X(2)
          CON(3)=X(3)-5.0*X(1)-X(2)
      ELSE IF (NPROB .EQ. 8) THEN
C
C     Test problem 8 (Rosen-Suzuki)
C
          F=X(1)**2+X(2)**2+2.0*X(3)**2+X(4)**2-5.0*X(1)-5.0*X(2)
     1      -21.0*X(3)+7.0*X(4)
          CON(1)=8.0-X(1)**2-X(2)**2-X(3)**2-X(4)**2-X(1)+X(2)
     1      -X(3)+X(4)
          CON(2)=10.0-X(1)**2-2.0*X(2)**2-X(3)**2-2.0*X(4)**2+X(1)+X(4)
          CON(3)=5.0-2.0*X(1)**2-X(2)**2-X(3)**2-2.0*X(1)+X(2)+X(4)
      ELSE IF (NPROB .EQ. 9) THEN
C
C     Test problem 9 (Hock and Schittkowski 100)
C
          F=(X(1)-10.0)**2+5.0*(X(2)-12.0)**2+X(3)**4+3.0*(X(4)-11.0)**2
     1      +10.0*X(5)**6+7.0*X(6)**2+X(7)**4-4.0*X(6)*X(7)-10.0*X(6)
     2      -8.0*X(7)
          CON(1)=127.0-2.0*X(1)**2-3.0*X(2)**4-X(3)-4.0*X(4)**2-5.0*X(5)
          CON(2)=282.0-7.0*X(1)-3.0*X(2)-10.0*X(3)**2-X(4)+X(5)
          CON(3)=196.0-23.0*X(1)-X(2)**2-6.0*X(6)**2+8.0*X(7)
          CON(4)=-4.0*X(1)**2-X(2)**2+3.0*X(1)*X(2)-2.0*X(3)**2-5.0*X(6)
     1      +11.0*X(7)
      ELSE IF (NPROB .EQ. 10) THEN
C
C     Test problem 10 (Hexagon area)
C
          F=-0.5*(X(1)*X(4)-X(2)*X(3)+X(3)*X(9)-X(5)*X(9)+X(5)*X(8)
     1      -X(6)*X(7))
          CON(1)=1.0-X(3)**2-X(4)**2
          CON(2)=1.0-X(9)**2
          CON(3)=1.0-X(5)**2-X(6)**2
          CON(4)=1.0-X(1)**2-(X(2)-X(9))**2
          CON(5)=1.0-(X(1)-X(5))**2-(X(2)-X(6))**2
          CON(6)=1.0-(X(1)-X(7))**2-(X(2)-X(8))**2
          CON(7)=1.0-(X(3)-X(5))**2-(X(4)-X(6))**2
          CON(8)=1.0-(X(3)-X(7))**2-(X(4)-X(8))**2
          CON(9)=1.0-X(7)**2-(X(8)-X(9))**2
          CON(10)=X(1)*X(4)-X(2)*X(3)
          CON(11)=X(3)*X(9)
          CON(12)=-X(5)*X(9)
          CON(13)=X(5)*X(8)-X(6)*X(7)
          CON(14)=X(9)
      END IF
      RETURN
      END
-------------------------------------------------------------------------------

      Output from test problem 1 (Simple quadratic)

  Normal return from subroutine COBYLA

  NFVALS =  37    F = 1.809750E-05    MAXCV = 0.000000E+00
  X =-1.000879E+00   3.220609E-03

    Least squares error in variables =    3.338389E-03

  Normal return from subroutine COBYLA

  NFVALS =  65    F = 1.153291E-07    MAXCV = 0.000000E+00
  X =-9.999341E-01   2.682342E-04

    Least squares error in variables =    2.762020E-04
 ------------------------------------------------------------------

      Output from test problem 2 (2D unit circle calculation)

  Normal return from subroutine COBYLA

  NFVALS =  37    F =-4.999994E-01    MAXCV = 2.026558E-06
  X = 7.062163E-01  -7.079976E-01

    Least squares error in variables =    1.259601E-03

  Normal return from subroutine COBYLA

  NFVALS =  44    F =-5.000000E-01    MAXCV = 5.960464E-08
  X = 7.071500E-01  -7.070636E-01

    Least squares error in variables =    6.107080E-05
 ------------------------------------------------------------------

      Output from test problem 3 (3D ellipsoid calculation)

  Normal return from subroutine COBYLA

  NFVALS =  52    F =-7.856688E-02    MAXCV = 6.288290E-06
  X = 5.780203E-01   4.069204E-01  -3.340311E-01

    Least squares error in variables =    1.642899E-03

  Normal return from subroutine COBYLA

  NFVALS =  65    F =-7.856742E-02    MAXCV = 8.940697E-08
  X = 5.773363E-01   4.082997E-01  -3.332995E-01

    Least squares error in variables =    6.312904E-05
 ------------------------------------------------------------------

      Output from test problem 4 (Weak Rosenbrock)

  Normal return from subroutine COBYLA

  NFVALS = 100    F = 3.125441E-05    MAXCV = 0.000000E+00
  X =-9.958720E-01   9.879909E-01

    Least squares error in variables =    1.269875E-02

  Normal return from subroutine COBYLA

  NFVALS = 173    F = 6.409362E-07    MAXCV = 0.000000E+00
  X =-9.994782E-01   9.983495E-01

    Least squares error in variables =    1.730967E-03
 ------------------------------------------------------------------

      Output from test problem 5 (Intermediate Rosenbrock)

  Normal return from subroutine COBYLA

  NFVALS = 347    F = 4.008353E-03    MAXCV = 0.000000E+00
  X =-9.366965E-01   8.777190E-01

    Least squares error in variables =    1.376952E-01

  Normal return from subroutine COBYLA

  NFVALS = 698    F = 9.516375E-05    MAXCV = 0.000000E+00
  X =-9.904447E-01   9.803594E-01

    Least squares error in variables =    2.184159E-02
 ------------------------------------------------------------------

      Output from test problem 6 (Equation (9.1.15) in Fletcher)

  Normal return from subroutine COBYLA

  NFVALS =  30    F =-1.414216E+00    MAXCV = 2.950430E-06
  X = 7.071948E-01   7.070208E-01

    Least squares error in variables =    1.230355E-04

  Normal return from subroutine COBYLA

  NFVALS =  40    F =-1.414214E+00    MAXCV = 0.000000E+00
  X = 7.071791E-01   7.070344E-01

    Least squares error in variables =    1.023325E-04
 ------------------------------------------------------------------

      Output from test problem 7 (Equation (14.4.2) in Fletcher)

  Normal return from subroutine COBYLA

  NFVALS =  28    F =-2.999881E+00    MAXCV = 0.000000E+00
  X = 1.362514E-08  -2.999881E+00  -2.999881E+00

    Least squares error in variables =    1.689246E-04

  Normal return from subroutine COBYLA

  NFVALS =  32    F =-3.000046E+00    MAXCV = 4.673004E-05
  X = 1.207445E-08  -3.000000E+00  -3.000046E+00

    Least squares error in variables =    4.649224E-05
 ------------------------------------------------------------------

      Output from test problem 8 (Rosen-Suzuki)

  Normal return from subroutine COBYLA

  NFVALS =  66    F =-4.400000E+01    MAXCV = 3.099442E-06
  X =-6.020486E-04   9.995968E-01   2.000546E+00  -9.994259E-01

    Least squares error in variables =    1.073541E-03

  Normal return from subroutine COBYLA

  NFVALS =  79    F =-4.400000E+01    MAXCV = 1.251698E-06
  X =-2.246869E-04   9.996516E-01   2.000260E+00  -9.997578E-01

    Least squares error in variables =    5.460466E-04
 ------------------------------------------------------------------

      Output from test problem 9 (Hock and Schittkowski 100)

  Normal return from subroutine COBYLA

  NFVALS = 237    F = 6.806309E+02    MAXCV = 0.000000E+00
  X = 2.332463E+00   1.951341E+00  -4.587620E-01   4.364742E+00  -6.244753E-01
      1.038812E+00   1.593632E+00

    Least squares error in variables =    1.892897E-02

  Normal return from subroutine COBYLA

  NFVALS = 248    F = 6.806310E+02    MAXCV = 1.907349E-05
  X = 2.332446E+00   1.951307E+00  -4.577461E-01   4.364753E+00  -6.241184E-01
      1.039491E+00   1.593760E+00

    Least squares error in variables =    1.996995E-02
 ------------------------------------------------------------------

      Output from test problem 10 (Hexagon area)

  Normal return from subroutine COBYLA

  NFVALS = 150    F =-8.660254E-01    MAXCV = 1.192093E-06
  X = 6.605685E-01   7.507660E-01  -3.188329E-01   9.478114E-01   6.614124E-01
      7.500232E-01  -3.198982E-01   9.474520E-01  -6.671554E-11

    Least squares error in variables =    1.124314E-03

  Normal return from subroutine COBYLA

  NFVALS = 173    F =-8.660254E-01    MAXCV = 3.352761E-07
  X = 6.606672E-01   7.506790E-01  -3.195507E-01   9.475691E-01   6.608437E-01
      7.505235E-01  -3.197733E-01   9.474941E-01  -3.822812E-11

    Least squares error in variables =    2.350494E-04
 ------------------------------------------------------------------