outlets/lincoa.f

      SUBROUTINE LINCOA (N,NPT,M,A,IA,B,X,RHOBEG,RHOEND,IPRINT,
     1  MAXFUN,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(IA,*),B(*),X(*),W(*)
C
C     This subroutine seeks the least value of a function of many variables,
C       subject to general linear inequality constraints, by a trust region
C       method that forms quadratic models by interpolation. Usually there
C       is much freedom in each new model after satisfying the interpolation
C       conditions, which is taken up by minimizing the Frobenius norm of
C       the change to the second derivative matrix of the model. One new
C       function value is calculated on each iteration, usually at a point
C       where the current model predicts a reduction in the least value so
C       far of the objective function subject to the linear constraints.
C       Alternatively, a new vector of variables may be chosen to replace
C       an interpolation point that may be too far away for reliability, and
C       then the new point does not have to satisfy the linear constraints.
C       The arguments of the subroutine are as follows.
C
C     N must be set to the number of variables and must be at least two.
C     NPT must be set to the number of interpolation conditions, which is
C       required to be in the interval [N+2,(N+1)(N+2)/2]. Typical choices
C       of the author are NPT=N+6 and NPT=2*N+1. Larger values tend to be
C       highly inefficent when the number of variables is substantial, due
C       to the amount of work and extra difficulty of adjusting more points.
C     M must be set to the number of linear inequality constraints.
C     A is a matrix whose columns are the constraint gradients, which are
C       required to be nonzero.
C     IA is the first dimension of the array A, which must be at least N.
C     B is the vector of right hand sides of the constraints, the J-th
C       constraint being that the scalar product of A(.,J) with X(.) is at
C       most B(J). The initial vector X(.) is made feasible by increasing
C       the value of B(J) if necessary.
C     X is the vector of variables. Initial values of X(1),X(2),...,X(N)
C       must be supplied. If they do not satisfy the constraints, then B
C       is increased as mentioned above. X contains on return the variables
C       that have given the least calculated F subject to the constraints.
C     RHOBEG and RHOEND must be set to the initial and final values of a
C       trust region radius, so both must be positive with RHOEND<=RHOBEG.
C       Typically, RHOBEG should be about one tenth of the greatest expected
C       change to a variable, and RHOEND should indicate the accuracy that
C       is required in the final values of the variables.
C     The value of IPRINT should be set to 0, 1, 2 or 3, which controls the
C       amount of printing. Specifically, there is no output if IPRINT=0 and
C       there is output only at the return if IPRINT=1. Otherwise, the best
C       feasible vector of variables so far and the corresponding value of
C       the objective function are printed whenever RHO is reduced, where
C       RHO is the current lower bound on the trust region radius. Further,
C       each new value of F with its variables are output if IPRINT=3.
C     MAXFUN must be set to an upper bound on the number of calls of CALFUN,
C       its value being at least NPT+1.
C     W is an array used for working space. Its length must be at least
C       M*(2+N) + NPT*(4+N+NPT) + N*(9+3*N) + MAX [ M+3*N, 2*M+N, 2*NPT ].
C       On return, W(1) is set to the final value of F, and W(2) is set to
C       the total number of function evaluations plus 0.5.
C
C     SUBROUTINE CALFUN (N,X,F) has to be provided by the user. It must set
C       F to the value of the objective function for the variables X(1),
C       X(2),...,X(N). The value of the argument F is positive when CALFUN
C       is called if and only if the current X satisfies the constraints
C       to working accuracy.
C
C     Check that N, NPT and MAXFUN are acceptable.
C
      ZERO=0.0D0
      SMALLX=1.0D-6*RHOEND
      NP=N+1
      NPTM=NPT-NP
      IF (N .LE. 1) THEN
          PRINT 10
   10     FORMAT (/4X,'Return from LINCOA because N is less than 2.')
          GOTO 80
      END IF
      IF (NPT .LT. N+2 .OR. NPT .GT. ((N+2)*NP)/2) THEN
          PRINT 20
   20     FORMAT (/4X,'Return from LINCOA because NPT is not in',
     1      ' the required interval.')
          GOTO 80
      END IF
      IF (MAXFUN .LE. NPT) THEN
          PRINT 30
   30     FORMAT (/4X,'Return from LINCOA because MAXFUN is less',
     1      ' than NPT+1.')
          GOTO 80
      END IF
C
C     Normalize the constraints, and copy the resultant constraint matrix
C       and right hand sides into working space, after increasing the right
C       hand sides if necessary so that the starting point is feasible.
C
      IAMAT=MAX0(M+3*N,2*M+N,2*NPT)+1
      IB=IAMAT+M*N
      IFLAG=0
      IF (M .GT. 0) THEN
          IW=IAMAT-1
          DO 60 J=1,M
          SUM=ZERO
          TEMP=ZERO
          DO 40 I=1,N
          SUM=SUM+A(I,J)*X(I)
   40     TEMP=TEMP+A(I,J)**2
          IF (TEMP .EQ. ZERO) THEN
              PRINT 50
   50         FORMAT (/4X,'Return from LINCOA because the gradient of',
     1          ' a constraint is zero.')
              GOTO 80
          END IF
          TEMP=DSQRT(TEMP)
          IF (SUM-B(J) .GT. SMALLX*TEMP) IFLAG=1
          W(IB+J-1)=DMAX1(B(J),SUM)/TEMP
          DO 60 I=1,N
          IW=IW+1
   60     W(IW)=A(I,J)/TEMP
      END IF
      IF (IFLAG .EQ. 1) THEN
          IF (IPRINT .GT. 0) PRINT 70
   70     FORMAT (/4X,'LINCOA has made the initial X feasible by',
     1      ' increasing part(s) of B.')
      END IF
C
C     Partition the working space array, so that different parts of it can be
C     treated separately by the subroutine that performs the main calculation.
C
      NDIM=NPT+N
      IXB=IB+M
      IXP=IXB+N
      IFV=IXP+N*NPT
      IXS=IFV+NPT
      IXO=IXS+N
      IGO=IXO+N
      IHQ=IGO+N
      IPQ=IHQ+(N*NP)/2
      IBMAT=IPQ+NPT
      IZMAT=IBMAT+NDIM*N
      ISTP=IZMAT+NPT*NPTM
      ISP=ISTP+N
      IXN=ISP+NPT+NPT
      IAC=IXN+N
      IRC=IAC+N
      IQF=IRC+M
      IRF=IQF+N*N
      IPQW=IRF+(N*NP)/2
C
C     The above settings provide a partition of W for subroutine LINCOB.
C
      CALL LINCOB (N,NPT,M,W(IAMAT),W(IB),X,RHOBEG,RHOEND,IPRINT,
     1  MAXFUN,W(IXB),W(IXP),W(IFV),W(IXS),W(IXO),W(IGO),W(IHQ),
     2  W(IPQ),W(IBMAT),W(IZMAT),NDIM,W(ISTP),W(ISP),W(IXN),W(IAC),
     3  W(IRC),W(IQF),W(IRF),W(IPQW),W)
   80 RETURN
      END

      SUBROUTINE GETACT (N,M,AMAT,B,NACT,IACT,QFAC,RFAC,SNORM,
     1  RESNEW,RESACT,G,DW,VLAM,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION AMAT(N,*),B(*),IACT(*),QFAC(N,*),RFAC(*),
     1  RESNEW(*),RESACT(*),G(*),DW(*),VLAM(*),W(*)
C
C     N, M, AMAT, B, NACT, IACT, QFAC and RFAC are the same as the terms
C       with these names in SUBROUTINE LINCOB. The current values must be
C       set on entry. NACT, IACT, QFAC and RFAC are kept up to date when
C       GETACT changes the current active set.
C     SNORM, RESNEW, RESACT, G and DW are the same as the terms with these
C       names in SUBROUTINE TRSTEP. The elements of RESNEW and RESACT are
C       also kept up to date.
C     VLAM and W are used for working space, the vector VLAM being reserved
C       for the Lagrange multipliers of the calculation. Their lengths must
C       be at least N.
C     The main purpose of GETACT is to pick the current active set. It is
C       defined by the property that the projection of -G into the space
C       orthogonal to the active constraint normals is as large as possible,
C       subject to this projected steepest descent direction moving no closer
C       to the boundary of every constraint whose current residual is at most
C       0.2*SNORM. On return, the settings in NACT, IACT, QFAC and RFAC are
C       all appropriate to this choice of active set.
C     Occasionally this projected direction is zero, and then the final value
C       of W(1) is set to zero. Otherwise, the direction itself is returned
C       in DW, and W(1) is set to the square of the length of the direction.
C
C     Set some constants and a temporary VLAM.
C
      ONE=1.0D0
      TINY=1.0D-60
      ZERO=0.0D0
      TDEL=0.2D0*SNORM
      DDSAV=ZERO
      DO 10 I=1,N
      DDSAV=DDSAV+G(I)**2
   10 VLAM(I)=ZERO
      DDSAV=DDSAV+DDSAV
C
C     Set the initial QFAC to the identity matrix in the case NACT=0.
C
      IF (NACT .EQ. 0) THEN
          DO 30 I=1,N
          DO 20 J=1,N
   20     QFAC(I,J)=ZERO
   30     QFAC(I,I)=ONE
          GOTO 100
      END IF
C
C     Remove any constraints from the initial active set whose residuals
C       exceed TDEL.
C
      IFLAG=1
      IC=NACT
   40 IF (RESACT(IC) .GT. TDEL) GOTO 800
   50 IC=IC-1
      IF (IC .GT. 0) GOTO 40
C
C     Remove any constraints from the initial active set whose Lagrange
C       multipliers are nonnegative, and set the surviving multipliers.
C
      IFLAG=2
   60 IF (NACT .EQ. 0) GOTO 100
      IC=NACT
   70 TEMP=ZERO
      DO 80 I=1,N
   80 TEMP=TEMP+QFAC(I,IC)*G(I)
      IDIAG=(IC*IC+IC)/2
      IF (IC .LT. NACT) THEN
          JW=IDIAG+IC
          DO 90 J=IC+1,NACT
          TEMP=TEMP-RFAC(JW)*VLAM(J)
   90     JW=JW+J
      END IF
      IF (TEMP .GE. ZERO) GOTO 800
      VLAM(IC)=TEMP/RFAC(IDIAG)
      IC=IC-1
      IF (IC .GT. 0) GOTO 70
C
C     Set the new search direction D. Terminate if the 2-norm of D is zero
C       or does not decrease, or if NACT=N holds. The situation NACT=N
C       occurs for sufficiently large SNORM if the origin is in the convex
C       hull of the constraint gradients.
C
  100 IF (NACT .EQ. N) GOTO 290
      DO 110 J=NACT+1,N
      W(J)=ZERO
      DO 110 I=1,N
  110 W(J)=W(J)+QFAC(I,J)*G(I)
      DD=ZERO
      DO 130 I=1,N
      DW(I)=ZERO
      DO 120 J=NACT+1,N
  120 DW(I)=DW(I)-W(J)*QFAC(I,J)
  130 DD=DD+DW(I)**2
      IF (DD .GE. DDSAV) GOTO 290
      IF (DD .EQ. ZERO) GOTO 300
      DDSAV=DD
      DNORM=DSQRT(DD)
C
C     Pick the next integer L or terminate, a positive value of L being
C       the index of the most violated constraint. The purpose of CTOL
C       below is to estimate whether a positive value of VIOLMX may be
C       due to computer rounding errors.
C
      L=0
      IF (M .GT. 0) THEN
          TEST=DNORM/SNORM
          VIOLMX=ZERO
          DO 150 J=1,M
          IF (RESNEW(J) .GT. ZERO .AND. RESNEW(J) .LE. TDEL) THEN
              SUM=ZERO
              DO 140 I=1,N
  140         SUM=SUM+AMAT(I,J)*DW(I)
              IF (SUM .GT. TEST*RESNEW(J)) THEN
                  IF (SUM .GT. VIOLMX) THEN
                      L=J
                      VIOLMX=SUM
                  END IF
              END IF
          END IF
  150     CONTINUE
          CTOL=ZERO
          TEMP=0.01D0*DNORM
          IF (VIOLMX .GT. ZERO .AND. VIOLMX .LT. TEMP) THEN
              IF (NACT .GT. 0) THEN
                  DO 170 K=1,NACT
                  J=IACT(K)
                  SUM=ZERO
                  DO 160 I=1,N
  160             SUM=SUM+DW(I)*AMAT(I,J)
  170             CTOL=DMAX1(CTOL,DABS(SUM))
              END IF
          END IF
      END IF
      W(1)=ONE
      IF (L .EQ. 0) GOTO 300
      IF (VIOLMX .LE. 10.0D0*CTOL) GOTO 300
C
C     Apply Givens rotations to the last (N-NACT) columns of QFAC so that
C       the first (NACT+1) columns of QFAC are the ones required for the
C       addition of the L-th constraint, and add the appropriate column
C       to RFAC.
C
      NACTP=NACT+1
      IDIAG=(NACTP*NACTP-NACTP)/2
      RDIAG=ZERO
      DO 200 J=N,1,-1
      SPROD=ZERO
      DO 180 I=1,N
  180 SPROD=SPROD+QFAC(I,J)*AMAT(I,L)
      IF (J .LE. NACT) THEN
          RFAC(IDIAG+J)=SPROD
      ELSE
          IF (DABS(RDIAG) .LE. 1.0D-20*DABS(SPROD)) THEN
              RDIAG=SPROD
          ELSE
              TEMP=DSQRT(SPROD*SPROD+RDIAG*RDIAG)
              COSV=SPROD/TEMP
              SINV=RDIAG/TEMP
              RDIAG=TEMP
              DO 190 I=1,N
              TEMP=COSV*QFAC(I,J)+SINV*QFAC(I,J+1)
              QFAC(I,J+1)=-SINV*QFAC(I,J)+COSV*QFAC(I,J+1)
  190         QFAC(I,J)=TEMP
          END IF
      END IF
  200 CONTINUE
      IF (RDIAG .LT. ZERO) THEN
          DO 210 I=1,N
  210     QFAC(I,NACTP)=-QFAC(I,NACTP)
      END IF
      RFAC(IDIAG+NACTP)=DABS(RDIAG)
      NACT=NACTP
      IACT(NACT)=L
      RESACT(NACT)=RESNEW(L)
      VLAM(NACT)=ZERO
      RESNEW(L)=ZERO
C
C     Set the components of the vector VMU in W.
C
  220 W(NACT)=ONE/RFAC((NACT*NACT+NACT)/2)**2
      IF (NACT .GT. 1) THEN
          DO 240 I=NACT-1,1,-1
          IDIAG=(I*I+I)/2
          JW=IDIAG+I
          SUM=ZERO
          DO 230 J=I+1,NACT
          SUM=SUM-RFAC(JW)*W(J)
  230     JW=JW+J
  240     W(I)=SUM/RFAC(IDIAG)
      END IF
C
C     Calculate the multiple of VMU to subtract from VLAM, and update VLAM.
C
      VMULT=VIOLMX
      IC=0
      J=1
  250 IF (J .LT. NACT) THEN
          IF (VLAM(J) .GE. VMULT*W(J)) THEN
              IC=J
              VMULT=VLAM(J)/W(J)
          END IF
          J=J+1
          GOTO 250
      END IF
      DO 260 J=1,NACT
  260 VLAM(J)=VLAM(J)-VMULT*W(J)
      IF (IC .GT. 0) VLAM(IC)=ZERO
      VIOLMX=DMAX1(VIOLMX-VMULT,ZERO)
      IF (IC .EQ. 0) VIOLMX=ZERO
C
C     Reduce the active set if necessary, so that all components of the
C       new VLAM are negative, with resetting of the residuals of the
C       constraints that become inactive.
C
      IFLAG=3
      IC=NACT
  270 IF (VLAM(IC) .LT. ZERO) GOTO 280
      RESNEW(IACT(IC))=DMAX1(RESACT(IC),TINY)
      GOTO 800
  280 IC=IC-1
      IF (IC .GT. 0) GOTO 270
C
C     Calculate the next VMU if VIOLMX is positive. Return if NACT=N holds,
C       as then the active constraints imply D=0. Otherwise, go to label
C       100, to calculate the new D and to test for termination.
C
      IF (VIOLMX .GT. ZERO) GOTO 220
      IF (NACT .LT. N) GOTO 100
  290 DD=ZERO
  300 W(1)=DD
      RETURN
C
C     These instructions rearrange the active constraints so that the new
C       value of IACT(NACT) is the old value of IACT(IC). A sequence of
C       Givens rotations is applied to the current QFAC and RFAC. Then NACT
C       is reduced by one.
C
  800 RESNEW(IACT(IC))=DMAX1(RESACT(IC),TINY)
      JC=IC
  810 IF (JC .LT. NACT) THEN
          JCP=JC+1
          IDIAG=JC*JCP/2
          JW=IDIAG+JCP
          TEMP=DSQRT(RFAC(JW-1)**2+RFAC(JW)**2)
          CVAL=RFAC(JW)/TEMP
          SVAL=RFAC(JW-1)/TEMP
          RFAC(JW-1)=SVAL*RFAC(IDIAG)
          RFAC(JW)=CVAL*RFAC(IDIAG)
          RFAC(IDIAG)=TEMP
          IF (JCP .LT. NACT) THEN
              DO 820 J=JCP+1,NACT
              TEMP=SVAL*RFAC(JW+JC)+CVAL*RFAC(JW+JCP)
              RFAC(JW+JCP)=CVAL*RFAC(JW+JC)-SVAL*RFAC(JW+JCP)
              RFAC(JW+JC)=TEMP
  820         JW=JW+J
          END IF
          JDIAG=IDIAG-JC
          DO 830 I=1,N
          IF (I .LT. JC) THEN
              TEMP=RFAC(IDIAG+I)
              RFAC(IDIAG+I)=RFAC(JDIAG+I)
              RFAC(JDIAG+I)=TEMP
          END IF
          TEMP=SVAL*QFAC(I,JC)+CVAL*QFAC(I,JCP)
          QFAC(I,JCP)=CVAL*QFAC(I,JC)-SVAL*QFAC(I,JCP)
  830     QFAC(I,JC)=TEMP
          IACT(JC)=IACT(JCP)
          RESACT(JC)=RESACT(JCP)
          VLAM(JC)=VLAM(JCP)
          JC=JCP
          GOTO 810
      END IF
      NACT=NACT-1
      GOTO (50,60,280),IFLAG
      END

      SUBROUTINE LINCOB (N,NPT,M,AMAT,B,X,RHOBEG,RHOEND,IPRINT,
     1  MAXFUN,XBASE,XPT,FVAL,XSAV,XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,
     2  STEP,SP,XNEW,IACT,RESCON,QFAC,RFAC,PQW,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION AMAT(N,*),B(*),X(*),XBASE(*),XPT(NPT,*),FVAL(*),
     1  XSAV(*),XOPT(*),GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),
     2  ZMAT(NPT,*),STEP(*),SP(*),XNEW(*),IACT(*),RESCON(*),
     3  QFAC(N,*),RFAC(*),PQW(*),W(*)
C
C     The arguments N, NPT, M, X, RHOBEG, RHOEND, IPRINT and MAXFUN are
C       identical to the corresponding arguments in SUBROUTINE LINCOA.
C     AMAT is a matrix whose columns are the constraint gradients, scaled
C       so that they have unit length.
C     B contains on entry the right hand sides of the constraints, scaled
C       as above, but later B is modified for variables relative to XBASE.
C     XBASE holds a shift of origin that should reduce the contributions
C       from rounding errors to values of the model and Lagrange functions.
C     XPT contains the interpolation point coordinates relative to XBASE.
C     FVAL holds the values of F at the interpolation points.
C     XSAV holds the best feasible vector of variables so far, without any
C       shift of origin.
C     XOPT is set to XSAV-XBASE, which is the displacement from XBASE of
C       the feasible vector of variables that provides the least calculated
C       F so far, this vector being the current trust region centre.
C     GOPT holds the gradient of the quadratic model at XSAV = XBASE+XOPT.
C     HQ holds the explicit second derivatives of the quadratic model.
C     PQ contains the parameters of the implicit second derivatives of the
C       quadratic model.
C     BMAT holds the last N columns of the big inverse matrix H.
C     ZMAT holds the factorization of the leading NPT by NPT submatrix
C       of H, this factorization being ZMAT times Diag(DZ) times ZMAT^T,
C       where the elements of DZ are plus or minus one, as specified by IDZ.
C     NDIM is the first dimension of BMAT and has the value NPT+N.
C     STEP is employed for trial steps from XOPT. It is also used for working
C       space when XBASE is shifted and in PRELIM.
C     SP is reserved for the scalar products XOPT^T XPT(K,.), K=1,2,...,NPT,
C       followed by STEP^T XPT(K,.), K=1,2,...,NPT.
C     XNEW is the displacement from XBASE of the vector of variables for
C       the current calculation of F, except that SUBROUTINE TRSTEP uses it
C       for working space.
C     IACT is an integer array for the indices of the active constraints.
C     RESCON holds useful information about the constraint residuals. Every
C       nonnegative RESCON(J) is the residual of the J-th constraint at the
C       current trust region centre. Otherwise, if RESCON(J) is negative, the
C       J-th constraint holds as a strict inequality at the trust region
C       centre, its residual being at least |RESCON(J)|; further, the value
C       of |RESCON(J)| is at least the current trust region radius DELTA.
C     QFAC is the orthogonal part of the QR factorization of the matrix of
C       active constraint gradients, these gradients being ordered in
C       accordance with IACT. When NACT is less than N, columns are added
C       to QFAC to complete an N by N orthogonal matrix, which is important
C       for keeping calculated steps sufficiently close to the boundaries
C       of the active constraints.
C     RFAC is the upper triangular part of this QR factorization, beginning
C       with the first diagonal element, followed by the two elements in the
C       upper triangular part of the second column and so on.
C     PQW is used for working space, mainly for storing second derivative
C       coefficients of quadratic functions. Its length is NPT+N.
C     The array W is also used for working space. The required number of
C       elements, namely MAX[M+3*N,2*M+N,2*NPT], is set in LINCOA.
C
C     Set some constants.
C
      HALF=0.5D0
      ONE=1.0D0
      TENTH=0.1D0
      ZERO=0.0D0
      NP=N+1
      NH=(N*NP)/2
      NPTM=NPT-NP
C
C     Set the elements of XBASE, XPT, FVAL, XSAV, XOPT, GOPT, HQ, PQ, BMAT,
C       ZMAT and SP for the first iteration. An important feature is that,
C       if the interpolation point XPT(K,.) is not feasible, where K is any
C       integer from [1,NPT], then a change is made to XPT(K,.) if necessary
C       so that the constraint violation is at least 0.2*RHOBEG. Also KOPT
C       is set so that XPT(KOPT,.) is the initial trust region centre.
C
      CALL PRELIM (N,NPT,M,AMAT,B,X,RHOBEG,IPRINT,XBASE,XPT,FVAL,
     1  XSAV,XOPT,GOPT,KOPT,HQ,PQ,BMAT,ZMAT,IDZ,NDIM,SP,RESCON,
     2  STEP,PQW,W)
C
C     Begin the iterative procedure.
C
      NF=NPT
      FOPT=FVAL(KOPT)
      RHO=RHOBEG
      DELTA=RHO
      IFEAS=0
      NACT=0
      ITEST=3
   10 KNEW=0
      NVALA=0
      NVALB=0
C
C     Shift XBASE if XOPT may be too far from XBASE. First make the changes
C       to BMAT that do not depend on ZMAT.
C
   20 FSAVE=FOPT
      XOPTSQ=ZERO
      DO 30 I=1,N
   30 XOPTSQ=XOPTSQ+XOPT(I)**2
      IF (XOPTSQ .GE. 1.0D4*DELTA*DELTA) THEN
          QOPTSQ=0.25D0*XOPTSQ
          DO 50 K=1,NPT
          SUM=ZERO
          DO 40 I=1,N
   40     SUM=SUM+XPT(K,I)*XOPT(I)
          SUM=SUM-HALF*XOPTSQ
          W(NPT+K)=SUM
          SP(K)=ZERO
          DO 50 I=1,N
          XPT(K,I)=XPT(K,I)-HALF*XOPT(I)
          STEP(I)=BMAT(K,I)
          W(I)=SUM*XPT(K,I)+QOPTSQ*XOPT(I)
          IP=NPT+I
          DO 50 J=1,I
   50     BMAT(IP,J)=BMAT(IP,J)+STEP(I)*W(J)+W(I)*STEP(J)
C
C     Then the revisions of BMAT that depend on ZMAT are calculated.
C
          DO 90 K=1,NPTM
          SUMZ=ZERO
          DO 60 I=1,NPT
          SUMZ=SUMZ+ZMAT(I,K)
   60     W(I)=W(NPT+I)*ZMAT(I,K)
          DO 80 J=1,N
          SUM=QOPTSQ*SUMZ*XOPT(J)
          DO 70 I=1,NPT
   70     SUM=SUM+W(I)*XPT(I,J)
          STEP(J)=SUM
          IF (K .LT. IDZ) SUM=-SUM
          DO 80 I=1,NPT
   80     BMAT(I,J)=BMAT(I,J)+SUM*ZMAT(I,K)
          DO 90 I=1,N
          IP=I+NPT
          TEMP=STEP(I)
          IF (K .LT. IDZ) TEMP=-TEMP
          DO 90 J=1,I
   90     BMAT(IP,J)=BMAT(IP,J)+TEMP*STEP(J)
C
C     Update the right hand sides of the constraints.
C
          IF (M .GT. 0) THEN
              DO 110 J=1,M
              TEMP=ZERO
              DO 100 I=1,N
  100         TEMP=TEMP+AMAT(I,J)*XOPT(I)
  110         B(J)=B(J)-TEMP
          END IF
C
C     The following instructions complete the shift of XBASE, including the
C       changes to the parameters of the quadratic model.
C
          IH=0
          DO 130 J=1,N
          W(J)=ZERO
          DO 120 K=1,NPT
          W(J)=W(J)+PQ(K)*XPT(K,J)
  120     XPT(K,J)=XPT(K,J)-HALF*XOPT(J)
          DO 130 I=1,J
          IH=IH+1
          HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J)
  130     BMAT(NPT+I,J)=BMAT(NPT+J,I)
          DO 140 J=1,N
          XBASE(J)=XBASE(J)+XOPT(J)
          XOPT(J)=ZERO
  140     XPT(KOPT,J)=ZERO
      END IF
C
C     In the case KNEW=0, generate the next trust region step by calling
C       TRSTEP, where SNORM is the current trust region radius initially.
C       The final value of SNORM is the length of the calculated step,
C       except that SNORM is zero on return if the projected gradient is
C       unsuitable for starting the conjugate gradient iterations.
C
      DELSAV=DELTA
      KSAVE=KNEW
      IF (KNEW .EQ. 0) THEN
          SNORM=DELTA
          DO 150 I=1,N
  150     XNEW(I)=GOPT(I)
          CALL TRSTEP (N,NPT,M,AMAT,B,XPT,HQ,PQ,NACT,IACT,RESCON,
     1      QFAC,RFAC,SNORM,STEP,XNEW,W,W(M+1),PQW,PQW(NP),W(M+NP))
C
C     A trust region step is applied whenever its length, namely SNORM, is at
C       least HALF*DELTA. It is also applied if its length is at least 0.1999
C       times DELTA and if a line search of TRSTEP has caused a change to the
C       active set. Otherwise there is a branch below to label 530 or 560.
C
          TEMP=HALF*DELTA
          IF (XNEW(1) .GE. HALF) TEMP=0.1999D0*DELTA
          IF (SNORM .LE. TEMP) THEN
              DELTA=HALF*DELTA
              IF (DELTA .LE. 1.4D0*RHO) DELTA=RHO
              NVALA=NVALA+1
              NVALB=NVALB+1
              TEMP=SNORM/RHO
              IF (DELSAV .GT. RHO) TEMP=ONE
              IF (TEMP .GE. HALF) NVALA=ZERO
              IF (TEMP .GE. TENTH) NVALB=ZERO
              IF (DELSAV .GT. RHO) GOTO 530
              IF (NVALA .LT. 5 .AND. NVALB .LT. 3) GOTO 530
              IF (SNORM .GT. ZERO) KSAVE=-1
              GOTO 560
          END IF
          NVALA=ZERO
          NVALB=ZERO
C
C     Alternatively, KNEW is positive. Then the model step is calculated
C       within a trust region of radius DEL, after setting the gradient at
C       XBASE and the second derivative parameters of the KNEW-th Lagrange
C       function in W(1) to W(N) and in PQW(1) to PQW(NPT), respectively.
C
      ELSE
          DEL=DMAX1(TENTH*DELTA,RHO)
          DO 160 I=1,N
  160     W(I)=BMAT(KNEW,I)
          DO 170 K=1,NPT
  170     PQW(K)=ZERO
          DO 180 J=1,NPTM
          TEMP=ZMAT(KNEW,J)
          IF (J .LT. IDZ) TEMP=-TEMP
          DO 180 K=1,NPT
  180     PQW(K)=PQW(K)+TEMP*ZMAT(K,J)
          CALL QMSTEP (N,NPT,M,AMAT,B,XPT,XOPT,NACT,IACT,RESCON,
     1      QFAC,KOPT,KNEW,DEL,STEP,W,PQW,W(NP),W(NP+M),IFEAS)
      END IF
C
C     Set VQUAD to the change to the quadratic model when the move STEP is
C       made from XOPT. If STEP is a trust region step, then VQUAD should be
C       negative. If it is nonnegative due to rounding errors in this case,
C       there is a branch to label 530 to try to improve the model.
C
      VQUAD=ZERO
      IH=0
      DO 190 J=1,N
      VQUAD=VQUAD+STEP(J)*GOPT(J)
      DO 190 I=1,J
      IH=IH+1
      TEMP=STEP(I)*STEP(J)
      IF (I .EQ. J) TEMP=HALF*TEMP
  190 VQUAD=VQUAD+TEMP*HQ(IH)
      DO 210 K=1,NPT
      TEMP=ZERO
      DO 200 J=1,N
      TEMP=TEMP+XPT(K,J)*STEP(J)
  200 SP(NPT+K)=TEMP
  210 VQUAD=VQUAD+HALF*PQ(K)*TEMP*TEMP
      IF (KSAVE .EQ. 0 .AND. VQUAD .GE. ZERO) GOTO 530
C
C     Calculate the next value of the objective function. The difference
C       between the actual new value of F and the value predicted by the
C       model is recorded in DIFF.
C
  220 NF=NF+1
      IF (NF .GT. MAXFUN) THEN
          NF=NF-1
          IF (IPRINT .GT. 0) PRINT 230
  230     FORMAT (/4X,'Return from LINCOA because CALFUN has been',
     1      ' called MAXFUN times.')
          GOTO 600
      END IF
      XDIFF=ZERO
      DO 240 I=1,N
      XNEW(I)=XOPT(I)+STEP(I)
      X(I)=XBASE(I)+XNEW(I)
  240 XDIFF=XDIFF+(X(I)-XSAV(I))**2
      XDIFF=DSQRT(XDIFF)
      IF (KSAVE .EQ. -1) XDIFF=RHO
      IF (XDIFF .LE. TENTH*RHO .OR. XDIFF .GE. DELTA+DELTA) THEN
          IFEAS=0
          IF (IPRINT .GT. 0) PRINT 250
  250     FORMAT (/4X,'Return from LINCOA because rounding errors',
     1      ' prevent reasonable changes to X.')
          GOTO 600
      END IF
      IF (KSAVE .LE. 0) IFEAS=1
      F=DFLOAT(IFEAS)
      CALL CALFUN (N,X,F)
      IF (IPRINT .EQ. 3) THEN
          PRINT 260, NF,F,(X(I),I=1,N)
  260     FORMAT (/4X,'Function number',I6,'    F =',1PD18.10,
     1      '    The corresponding X is:'/(2X,5D15.6))
      END IF
      IF (KSAVE .EQ. -1) GOTO 600
      DIFF=F-FOPT-VQUAD
C
C     If X is feasible, then set DFFALT to the difference between the new
C       value of F and the value predicted by the alternative model.
C
      IF (IFEAS .EQ. 1 .AND. ITEST .LT. 3) THEN
          DO 270 K=1,NPT
          PQW(K)=ZERO
  270     W(K)=FVAL(K)-FVAL(KOPT)
          DO 290 J=1,NPTM
          SUM=ZERO
          DO 280 I=1,NPT
  280     SUM=SUM+W(I)*ZMAT(I,J)
          IF (J .LT. IDZ) SUM=-SUM
          DO 290 K=1,NPT
  290     PQW(K)=PQW(K)+SUM*ZMAT(K,J)
          VQALT=ZERO
          DO 310 K=1,NPT
          SUM=ZERO
          DO 300 J=1,N
  300     SUM=SUM+BMAT(K,J)*STEP(J)
          VQALT=VQALT+SUM*W(K)
  310     VQALT=VQALT+PQW(K)*SP(NPT+K)*(HALF*SP(NPT+K)+SP(K))
          DFFALT=F-FOPT-VQALT
      END IF
      IF (ITEST .EQ. 3) THEN
          DFFALT=DIFF
          ITEST=0
      END IF
C
C     Pick the next value of DELTA after a trust region step.
C
      IF (KSAVE .EQ. 0) THEN
          RATIO=(F-FOPT)/VQUAD
          IF (RATIO .LE. TENTH) THEN
              DELTA=HALF*DELTA
          ELSE IF (RATIO .LE. 0.7D0) THEN
              DELTA=DMAX1(HALF*DELTA,SNORM)
          ELSE 
              TEMP=DSQRT(2.0D0)*DELTA
              DELTA=DMAX1(HALF*DELTA,SNORM+SNORM)
              DELTA=DMIN1(DELTA,TEMP)
          END IF
          IF (DELTA .LE. 1.4D0*RHO) DELTA=RHO
      END IF
C
C     Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point
C       can be moved. If STEP is a trust region step, then KNEW is zero at
C       present, but a positive value is picked by subroutine UPDATE.
C
      CALL UPDATE (N,NPT,XPT,BMAT,ZMAT,IDZ,NDIM,SP,STEP,KOPT,
     1  KNEW,PQW,W)
      IF (KNEW .EQ. 0) THEN
          IF (IPRINT .GT. 0) PRINT 320
  320     FORMAT (/4X,'Return from LINCOA because the denominator'
     1      ' of the updating formula is zero.')
          GOTO 600
      END IF
C
C     If ITEST is increased to 3, then the next quadratic model is the
C       one whose second derivative matrix is least subject to the new
C       interpolation conditions. Otherwise the new model is constructed
C       by the symmetric Broyden method in the usual way.
C
      IF (IFEAS .EQ. 1) THEN
          ITEST=ITEST+1
          IF (DABS(DFFALT) .GE. TENTH*DABS(DIFF)) ITEST=0
      END IF
C
C     Update the second derivatives of the model by the symmetric Broyden
C       method, using PQW for the second derivative parameters of the new
C       KNEW-th Lagrange function. The contribution from the old parameter
C       PQ(KNEW) is included in the second derivative matrix HQ. W is used
C       later for the gradient of the new KNEW-th Lagrange function.       
C
      IF (ITEST .LT. 3) THEN
          DO 330 K=1,NPT
  330     PQW(K)=ZERO
          DO 350 J=1,NPTM
          TEMP=ZMAT(KNEW,J)
          IF (TEMP .NE. ZERO) THEN
              IF (J .LT. IDZ) TEMP=-TEMP
              DO 340 K=1,NPT
  340         PQW(K)=PQW(K)+TEMP*ZMAT(K,J)
          END IF
  350     CONTINUE
          IH=0
          DO 360 I=1,N
          W(I)=BMAT(KNEW,I)
          TEMP=PQ(KNEW)*XPT(KNEW,I)
          DO 360 J=1,I
          IH=IH+1
  360     HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J)
          PQ(KNEW)=ZERO
          DO 370 K=1,NPT
  370     PQ(K)=PQ(K)+DIFF*PQW(K)
      END IF
C
C     Include the new interpolation point with the corresponding updates of
C       SP. Also make the changes of the symmetric Broyden method to GOPT at
C       the old XOPT if ITEST is less than 3.
C
      FVAL(KNEW)=F
      SP(KNEW)=SP(KOPT)+SP(NPT+KOPT)
      SSQ=ZERO
      DO 380 I=1,N
      XPT(KNEW,I)=XNEW(I)
  380 SSQ=SSQ+STEP(I)**2
      SP(NPT+KNEW)=SP(NPT+KOPT)+SSQ
      IF (ITEST .LT. 3) THEN
          DO 390 K=1,NPT
          TEMP=PQW(K)*SP(K)
          DO 390 I=1,N
  390     W(I)=W(I)+TEMP*XPT(K,I)
          DO 400 I=1,N
  400     GOPT(I)=GOPT(I)+DIFF*W(I)
      END IF
C
C     Update FOPT, XSAV, XOPT, KOPT, RESCON and SP if the new F is the
C       least calculated value so far with a feasible vector of variables.
C
      IF (F .LT. FOPT .AND. IFEAS .EQ. 1) THEN
          FOPT=F
          DO 410 J=1,N
          XSAV(J)=X(J)
  410     XOPT(J)=XNEW(J)
          KOPT=KNEW
          SNORM=DSQRT(SSQ)
          DO 430 J=1,M
          IF (RESCON(J) .GE. DELTA+SNORM) THEN
              RESCON(J)=SNORM-RESCON(J)
          ELSE
              RESCON(J)=RESCON(J)+SNORM
              IF (RESCON(J)+DELTA .GT. ZERO) THEN
                  TEMP=B(J)
                  DO 420 I=1,N
  420             TEMP=TEMP-XOPT(I)*AMAT(I,J)
                  TEMP=DMAX1(TEMP,ZERO)
                  IF (TEMP .GE. DELTA) TEMP=-TEMP
                  RESCON(J)=TEMP
              END IF
          END IF
  430     CONTINUE
          DO 440 K=1,NPT
  440     SP(K)=SP(K)+SP(NPT+K)
C
C     Also revise GOPT when symmetric Broyden updating is applied.
C
          IF (ITEST .LT. 3) THEN
              IH=0
              DO 450 J=1,N
              DO 450 I=1,J
              IH=IH+1
              IF (I .LT. J) GOPT(J)=GOPT(J)+HQ(IH)*STEP(I)
  450         GOPT(I)=GOPT(I)+HQ(IH)*STEP(J)
              DO 460 K=1,NPT
              TEMP=PQ(K)*SP(NPT+K)
              DO 460 I=1,N
  460         GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
          END IF
      END IF
C
C     Replace the current model by the least Frobenius norm interpolant if
C       this interpolant gives substantial reductions in the predictions
C       of values of F at feasible points.
C
      IF (ITEST .EQ. 3) THEN
          DO 470 K=1,NPT
          PQ(K)=ZERO
  470     W(K)=FVAL(K)-FVAL(KOPT)
          DO 490 J=1,NPTM
          SUM=ZERO
          DO 480 I=1,NPT
  480     SUM=SUM+W(I)*ZMAT(I,J)
          IF (J .LT. IDZ) SUM=-SUM
          DO 490 K=1,NPT
  490     PQ(K)=PQ(K)+SUM*ZMAT(K,J)
          DO 500 J=1,N
          GOPT(J)=ZERO
          DO 500 I=1,NPT
  500     GOPT(J)=GOPT(J)+W(I)*BMAT(I,J)
          DO 510 K=1,NPT
          TEMP=PQ(K)*SP(K)
          DO 510 I=1,N
  510     GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
          DO 520 IH=1,NH
  520     HQ(IH)=ZERO
      END IF
C
C     If a trust region step has provided a sufficient decrease in F, then
C       branch for another trust region calculation. Every iteration that
C       takes a model step is followed by an attempt to take a trust region
C       step.
C
      KNEW=0
      IF (KSAVE .GT. 0) GOTO 20
      IF (RATIO .GE. TENTH) GOTO 20
C
C     Alternatively, find out if the interpolation points are close enough
C       to the best point so far.
C
  530 DISTSQ=DMAX1(DELTA*DELTA,4.0D0*RHO*RHO)
      DO 550 K=1,NPT
      SUM=ZERO
      DO 540 J=1,N
  540 SUM=SUM+(XPT(K,J)-XOPT(J))**2
      IF (SUM .GT. DISTSQ) THEN
          KNEW=K
          DISTSQ=SUM
      END IF
  550 CONTINUE
C
C     If KNEW is positive, then branch back for the next iteration, which
C       will generate a "model step". Otherwise, if the current iteration
C       has reduced F, or if DELTA was above its lower bound when the last
C       trust region step was calculated, then try a "trust region" step
C       instead.
C
      IF (KNEW .GT. 0) GOTO 20
      KNEW=0
      IF (FOPT .LT. FSAVE) GOTO 20
      IF (DELSAV .GT. RHO) GOTO 20
C
C     The calculations with the current value of RHO are complete.
C       Pick the next value of RHO.
C
  560 IF (RHO .GT. RHOEND) THEN
          DELTA=HALF*RHO
          IF (RHO .GT. 250.0D0*RHOEND) THEN
              RHO=TENTH*RHO
          ELSE IF (RHO .LE. 16.0D0*RHOEND) THEN
              RHO=RHOEND
          ELSE
              RHO=DSQRT(RHO*RHOEND)
          END IF 
          DELTA=DMAX1(DELTA,RHO)
          IF (IPRINT .GE. 2) THEN
              IF (IPRINT .GE. 3) PRINT 570
  570         FORMAT (5X)
              PRINT 580, RHO,NF
  580         FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of',
     1          ' function values =',I6)
              PRINT 590, FOPT,(XBASE(I)+XOPT(I),I=1,N)
  590         FORMAT (4X,'Least value of F =',1PD23.15,9X,
     1          'The corresponding X is:'/(2X,5D15.6))
          END IF
          GOTO 10
      END IF
C
C     Return from the calculation, after branching to label 220 for another
C       Newton-Raphson step if it has not been tried before.
C
      IF (KSAVE .EQ. -1) GOTO 220
  600 IF (FOPT .LE. F .OR. IFEAS .EQ. 0) THEN
          DO 610 I=1,N
  610     X(I)=XSAV(I)
          F=FOPT
      END IF
      IF (IPRINT .GE. 1) THEN
          PRINT 620, NF
  620     FORMAT (/4X,'At the return from LINCOA',5X,
     1      'Number of function values =',I6)
          PRINT 590, F,(X(I),I=1,N)
      END IF
      W(1)=F
      W(2)=DFLOAT(NF)+HALF
      RETURN
      END

      SUBROUTINE PRELIM (N,NPT,M,AMAT,B,X,RHOBEG,IPRINT,XBASE,
     1  XPT,FVAL,XSAV,XOPT,GOPT,KOPT,HQ,PQ,BMAT,ZMAT,IDZ,NDIM,
     2  SP,RESCON,STEP,PQW,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION AMAT(N,*),B(*),X(*),XBASE(*),XPT(NPT,*),FVAL(*),
     1  XSAV(*),XOPT(*),GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),
     2  SP(*),RESCON(*),STEP(*),PQW(*),W(*)
C
C     The arguments N, NPT, M, AMAT, B, X, RHOBEG, IPRINT, XBASE, XPT, FVAL,
C       XSAV, XOPT, GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SP and RESCON are the
C       same as the corresponding arguments in SUBROUTINE LINCOB.
C     KOPT is set to the integer such that XPT(KOPT,.) is the initial trust
C       region centre.
C     IDZ is going to be set to one, so that every element of Diag(DZ) is
C       one in the product ZMAT times Diag(DZ) times ZMAT^T, which is the
C       factorization of the leading NPT by NPT submatrix of H.
C     STEP, PQW and W are used for working space, the arrays STEP and PQW
C       being taken from LINCOB. The length of W must be at least N+NPT.
C
C     SUBROUTINE PRELIM provides the elements of XBASE, XPT, BMAT and ZMAT
C       for the first iteration, an important feature being that, if any of
C       of the columns of XPT is an infeasible point, then the largest of
C       the constraint violations there is at least 0.2*RHOBEG. It also sets
C       the initial elements of FVAL, XOPT, GOPT, HQ, PQ, SP and RESCON.
C
C     Set some constants.
C
      HALF=0.5D0
      ONE=1.0D0
      ZERO=0.0D0
      NPTM=NPT-N-1
      RHOSQ=RHOBEG*RHOBEG
      RECIP=ONE/RHOSQ
      RECIQ=DSQRT(HALF)/RHOSQ
      TEST=0.2D0*RHOBEG
      IDZ=1
      KBASE=1
C
C     Set the initial elements of XPT, BMAT, SP and ZMAT to zero. 
C
      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 K=1,NPT
      SP(K)=ZERO
      DO 30 J=1,NPT-N-1
   30 ZMAT(K,J)=ZERO
C
C     Set the nonzero coordinates of XPT(K,.), K=1,2,...,min[2*N+1,NPT],
C       but they may be altered later to make a constraint violation
C       sufficiently large. The initial nonzero elements of BMAT and of
C       the first min[N,NPT-N-1] columns of ZMAT are set also.
C
      DO 40 J=1,N
      XPT(J+1,J)=RHOBEG
      IF (J .LT. NPT-N) THEN
          JP=N+J+1
          XPT(JP,J)=-RHOBEG
          BMAT(J+1,J)=HALF/RHOBEG
          BMAT(JP,J)=-HALF/RHOBEG
          ZMAT(1,J)=-RECIQ-RECIQ
          ZMAT(J+1,J)=RECIQ
          ZMAT(JP,J)=RECIQ
      ELSE
          BMAT(1,J)=-ONE/RHOBEG
          BMAT(J+1,J)=ONE/RHOBEG
          BMAT(NPT+J,J)=-HALF*RHOSQ
      END IF
   40 CONTINUE
C
C     Set the remaining initial nonzero elements of XPT and ZMAT when the
C       number of interpolation points exceeds 2*N+1.
C
      IF (NPT .GT. 2*N+1) THEN
          DO 50 K=N+1,NPT-N-1
          ITEMP=(K-1)/N
          IPT=K-ITEMP*N
          JPT=IPT+ITEMP
          IF (JPT .GT. N) JPT=JPT-N
          XPT(N+K+1,IPT)=RHOBEG
          XPT(N+K+1,JPT)=RHOBEG
          ZMAT(1,K)=RECIP
          ZMAT(IPT+1,K)=-RECIP
          ZMAT(JPT+1,K)=-RECIP
   50     ZMAT(N+K+1,K)=RECIP
      END IF
C
C     Update the constraint right hand sides to allow for the shift XBASE.
C
      IF (M .GT. 0) THEN
          DO 70 J=1,M
          TEMP=ZERO
          DO 60 I=1,N
   60     TEMP=TEMP+AMAT(I,J)*XBASE(I)
   70     B(J)=B(J)-TEMP
      END IF
C
C     Go through the initial points, shifting every infeasible point if
C       necessary so that its constraint violation is at least 0.2*RHOBEG.
C
      DO 150 NF=1,NPT
      FEAS=ONE
      BIGV=ZERO
      J=0
   80 J=J+1
      IF (J .LE. M .AND. NF .GE. 2) THEN
          RESID=-B(J)
          DO 90 I=1,N
   90     RESID=RESID+XPT(NF,I)*AMAT(I,J)
          IF (RESID .LE. BIGV) GOTO 80
          BIGV=RESID
          JSAV=J
          IF (RESID .LE. TEST) THEN
              FEAS=-ONE
              GOTO 80
          END IF
          FEAS=ZERO
      END IF
      IF (FEAS .LT. ZERO) THEN
          DO 100 I=1,N
  100     STEP(I)=XPT(NF,I)+(TEST-BIGV)*AMAT(I,JSAV)
          DO 110 K=1,NPT
          SP(NPT+K)=ZERO
          DO 110 J=1,N
  110     SP(NPT+K)=SP(NPT+K)+XPT(K,J)*STEP(J)
          CALL UPDATE (N,NPT,XPT,BMAT,ZMAT,IDZ,NDIM,SP,STEP,
     1      KBASE,NF,PQW,W)
          DO 120 I=1,N
  120     XPT(NF,I)=STEP(I)
      END IF
C
C     Calculate the objective function at the current interpolation point,
C       and set KOPT to the index of the first trust region centre.
C
      DO 130 J=1,N
  130 X(J)=XBASE(J)+XPT(NF,J)
      F=FEAS
      CALL CALFUN (N,X,F)
      IF (IPRINT .EQ. 3) THEN
          PRINT 140, NF,F,(X(I),I=1,N)
  140     FORMAT (/4X,'Function number',I6,'    F =',1PD18.10,
     1      '    The corresponding X is:'/(2X,5D15.6))
      END IF
      IF (NF .EQ. 1) THEN
          KOPT=1
      ELSE IF (F .LT. FVAL(KOPT) .AND. FEAS .GT. ZERO) THEN
          KOPT=NF
      END IF
  150 FVAL(NF)=F
C
C     Set PQ for the first quadratic model.
C
      DO 160 J=1,NPTM
      W(J)=ZERO
      DO 160 K=1,NPT
  160 W(J)=W(J)+ZMAT(K,J)*FVAL(K)
      DO 170 K=1,NPT
      PQ(K)=ZERO
      DO 170 J=1,NPTM
  170 PQ(K)=PQ(K)+ZMAT(K,J)*W(J)
C
C     Set XOPT, SP, GOPT and HQ for the first quadratic model.
C
      DO 180 J=1,N
      XOPT(J)=XPT(KOPT,J)
      XSAV(J)=XBASE(J)+XOPT(J)
  180 GOPT(J)=ZERO
      DO 200 K=1,NPT
      SP(K)=ZERO
      DO 190 J=1,N
  190 SP(K)=SP(K)+XPT(K,J)*XOPT(J)
      TEMP=PQ(K)*SP(K)
      DO 200 J=1,N
  200 GOPT(J)=GOPT(J)+FVAL(K)*BMAT(K,J)+TEMP*XPT(K,J)
      DO 210 I=1,(N*N+N)/2
  210 HQ(I)=ZERO
C
C     Set the initial elements of RESCON.
C
      DO 230 J=1,M
      TEMP=B(J)
      DO 220 I=1,N
  220 TEMP=TEMP-XOPT(I)*AMAT(I,J)
      TEMP=DMAX1(TEMP,ZERO)
      IF (TEMP .GE. RHOBEG) TEMP=-TEMP
  230 RESCON(J)=TEMP  
      RETURN
      END

      SUBROUTINE QMSTEP (N,NPT,M,AMAT,B,XPT,XOPT,NACT,IACT,
     1  RESCON,QFAC,KOPT,KNEW,DEL,STEP,GL,PQW,RSTAT,W,IFEAS)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION AMAT(N,*),B(*),XPT(NPT,*),XOPT(*),IACT(*),
     1  RESCON(*),QFAC(N,*),STEP(*),GL(*),PQW(*),RSTAT(*),W(*)
C
C     N, NPT, M, AMAT, B, XPT, XOPT, NACT, IACT, RESCON, QFAC, KOPT are the
C       same as the terms with these names in SUBROUTINE LINCOB.
C     KNEW is the index of the interpolation point that is going to be moved.
C     DEL is the current restriction on the length of STEP, which is never
C       greater than the current trust region radius DELTA.
C     STEP will be set to the required step from XOPT to the new point.
C     GL must be set on entry to the gradient of LFUNC at XBASE, where LFUNC
C       is the KNEW-th Lagrange function. It is used also for some other
C       gradients of LFUNC.
C     PQW provides the second derivative parameters of LFUNC.
C     RSTAT and W are used for working space. Their lengths must be at least
C       M and N, respectively. RSTAT(J) is set to -1.0, 0.0, or 1.0 if the
C       J-th constraint is irrelevant, active, or both inactive and relevant,
C       respectively.
C     IFEAS will be set to 0 or 1 if XOPT+STEP is infeasible or feasible.
C
C     STEP is chosen to provide a relatively large value of the modulus of
C       LFUNC(XOPT+STEP), subject to ||STEP|| .LE. DEL. A projected STEP is
C       calculated too, within the trust region, that does not alter the
C       residuals of the active constraints. The projected step is preferred
C       if its value of | LFUNC(XOPT+STEP) | is at least one fifth of the
C       original one, but the greatest violation of a linear constraint must
C       be at least 0.2*DEL, in order to keep the interpolation points apart.
C       The remedy when the maximum constraint violation is too small is to
C       restore the original step, which is perturbed if necessary so that
C       its maximum constraint violation becomes 0.2*DEL.
C
C     Set some constants.
C
      HALF=0.5D0
      ONE=1.0D0
      TENTH=0.1D0
      ZERO=0.0D0
      TEST=0.2D0*DEL
C
C     Replace GL by the gradient of LFUNC at the trust region centre, and
C       set the elements of RSTAT.
C
      DO 20 K=1,NPT
      TEMP=ZERO
      DO 10 J=1,N
   10 TEMP=TEMP+XPT(K,J)*XOPT(J)
      TEMP=PQW(K)*TEMP
      DO 20 I=1,N
   20 GL(I)=GL(I)+TEMP*XPT(K,I)
      IF (M .GT. 0) THEN
          DO 30 J=1,M
          RSTAT(J)=ONE
   30     IF (DABS(RESCON(J)) .GE. DEL) RSTAT(J)=-ONE
          DO 40 K=1,NACT
   40     RSTAT(IACT(K))=ZERO
      END IF
C
C     Find the greatest modulus of LFUNC on a line through XOPT and
C       another interpolation point within the trust region.
C
      IFLAG=0
      VBIG=ZERO
      DO 60 K=1,NPT
      IF (K .EQ. KOPT) GOTO 60
      SS=ZERO
      SP=ZERO
      DO 50 I=1,N
      TEMP=XPT(K,I)-XOPT(I)
      SS=SS+TEMP*TEMP
   50 SP=SP+GL(I)*TEMP
      STP=-DEL/DSQRT(SS)
      IF (K .EQ. KNEW) THEN
          IF (SP*(SP-ONE) .LT. ZERO) STP=-STP
          VLAG=DABS(STP*SP)+STP*STP*DABS(SP-ONE)
      ELSE
          VLAG=DABS(STP*(ONE-STP)*SP)
      END IF
      IF (VLAG .GT. VBIG) THEN
          KSAV=K
          STPSAV=STP
          VBIG=VLAG
      END IF
   60 CONTINUE
C
C     Set STEP to the move that gives the greatest modulus calculated above.
C       This move may be replaced by a steepest ascent step from XOPT.
C
      GG=ZERO
      DO 70 I=1,N
      GG=GG+GL(I)**2
   70 STEP(I)=STPSAV*(XPT(KSAV,I)-XOPT(I))
      VGRAD=DEL*DSQRT(GG)
      IF (VGRAD .LE. TENTH*VBIG) GOTO 220
C
C     Make the replacement if it provides a larger value of VBIG.
C
      GHG=ZERO
      DO 90 K=1,NPT
      TEMP=ZERO
      DO 80 J=1,N
   80 TEMP=TEMP+XPT(K,J)*GL(J)
   90 GHG=GHG+PQW(K)*TEMP*TEMP
      VNEW=VGRAD+DABS(HALF*DEL*DEL*GHG/GG)
      IF (VNEW .GT. VBIG) THEN
          VBIG=VNEW
          STP=DEL/DSQRT(GG)
          IF (GHG .LT. ZERO) STP=-STP
          DO 100 I=1,N
  100     STEP(I)=STP*GL(I)
      END IF
      IF (NACT .EQ. 0 .OR. NACT .EQ. N) GOTO 220
C
C     Overwrite GL by its projection. Then set VNEW to the greatest
C       value of |LFUNC| on the projected gradient from XOPT subject to
C       the trust region bound. If VNEW is sufficiently large, then STEP
C       may be changed to a move along the projected gradient.
C
      DO 110 K=NACT+1,N
      W(K)=ZERO
      DO 110 I=1,N
  110 W(K)=W(K)+GL(I)*QFAC(I,K)
      GG=ZERO
      DO 130 I=1,N
      GL(I)=ZERO
      DO 120 K=NACT+1,N
  120 GL(I)=GL(I)+QFAC(I,K)*W(K)
  130 GG=GG+GL(I)**2
      VGRAD=DEL*DSQRT(GG)
      IF (VGRAD .LE. TENTH*VBIG) GOTO 220
      GHG=ZERO
      DO 150 K=1,NPT
      TEMP=ZERO
      DO 140 J=1,N
  140 TEMP=TEMP+XPT(K,J)*GL(J)
  150 GHG=GHG+PQW(K)*TEMP*TEMP
      VNEW=VGRAD+DABS(HALF*DEL*DEL*GHG/GG)
C
C     Set W to the possible move along the projected gradient.
C
      STP=DEL/DSQRT(GG)
      IF (GHG .LT. ZERO) STP=-STP
      WW=ZERO
      DO 160 I=1,N
      W(I)=STP*GL(I)
  160 WW=WW+W(I)**2
C
C     Set STEP to W if W gives a sufficiently large value of the modulus
C       of the Lagrange function, and if W either preserves feasibility
C       or gives a constraint violation of at least 0.2*DEL. The purpose
C       of CTOL below is to provide a check on feasibility that includes
C       a tolerance for contributions from computer rounding errors.
C
      IF (VNEW/VBIG .GE. 0.2D0) THEN
          IFEAS=1
          BIGV=ZERO
          J=0
  170     J=J+1
          IF (J .LE. M) THEN
              IF (RSTAT(J) .EQ. ONE) THEN
                  TEMP=-RESCON(J)
                  DO 180 I=1,N
  180             TEMP=TEMP+W(I)*AMAT(I,J)
                  BIGV=DMAX1(BIGV,TEMP)
              END IF
              IF (BIGV .LT. TEST) GOTO 170
              IFEAS=0
          END IF
          CTOL=ZERO
          TEMP=0.01D0*DSQRT(WW)
          IF (BIGV .GT. ZERO .AND. BIGV .LT. TEMP) THEN
              DO 200 K=1,NACT
              J=IACT(K)
              SUM=ZERO
              DO 190 I=1,N
  190         SUM=SUM+W(I)*AMAT(I,J)
  200         CTOL=DMAX1(CTOL,DABS(SUM))
          END IF
          IF (BIGV .LE. 10.0D0*CTOL .OR. BIGV .GE. TEST) THEN
              DO 210 I=1,N
  210         STEP(I)=W(I)
              GOTO 260
          END IF
      END IF
C
C     Calculate the greatest constraint violation at XOPT+STEP with STEP at
C       its original value. Modify STEP if this violation is unacceptable.
C
  220 IFEAS=1
      BIGV=ZERO
      RESMAX=ZERO
      J=0
  230 J=J+1
      IF (J .LE. M) THEN
          IF (RSTAT(J) .LT. ZERO) GOTO 230
          TEMP=-RESCON(J)
          DO 240 I=1,N
  240     TEMP=TEMP+STEP(I)*AMAT(I,J)
          RESMAX=DMAX1(RESMAX,TEMP)
          IF (TEMP .LT. TEST) THEN
              IF (TEMP .LE. BIGV) GOTO 230
              BIGV=TEMP
              JSAV=J
              IFEAS=-1
              GOTO 230
          END IF
          IFEAS=0
      END IF
      IF (IFEAS .EQ. -1) THEN
          DO 250 I=1,N
  250     STEP(I)=STEP(I)+(TEST-BIGV)*AMAT(I,JSAV)
          IFEAS=0
      END IF
C
C     Return the calculated STEP and the value of IFEAS.
C
  260 RETURN
      END

      SUBROUTINE TRSTEP (N,NPT,M,AMAT,B,XPT,HQ,PQ,NACT,IACT,RESCON,
     1  QFAC,RFAC,SNORM,STEP,G,RESNEW,RESACT,D,DW,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION AMAT(N,*),B(*),XPT(NPT,*),HQ(*),PQ(*),IACT(*),
     1  RESCON(*),QFAC(N,*),RFAC(*),STEP(*),G(*),RESNEW(*),RESACT(*),
     2  D(*),DW(*),W(*)
C
C     N, NPT, M, AMAT, B, XPT, HQ, PQ, NACT, IACT, RESCON, QFAC and RFAC
C       are the same as the terms with these names in LINCOB. If RESCON(J)
C       is negative, then |RESCON(J)| must be no less than the trust region
C       radius, so that the J-th constraint can be ignored.
C     SNORM is set to the trust region radius DELTA initially. On the
C       return, however, it is the length of the calculated STEP, which is
C       set to zero if the constraints do not allow a long enough step.
C     STEP is the total calculated step so far from the trust region centre,
C       its final value being given by the sequence of CG iterations, which
C       terminate if the trust region boundary is reached.
C     G must be set on entry to the gradient of the quadratic model at the
C       trust region centre. It is used as working space, however, and is
C       always the gradient of the model at the current STEP, except that
C       on return the value of G(1) is set to ONE instead of to ZERO if
C       and only if GETACT is called more than once.
C     RESNEW, RESACT, D, DW and W are used for working space. A negative
C       value of RESNEW(J) indicates that the J-th constraint does not
C       restrict the CG steps of the current trust region calculation, a
C       zero value of RESNEW(J) indicates that the J-th constraint is active,
C       and otherwise RESNEW(J) is set to the greater of TINY and the actual
C       residual of the J-th constraint for the current STEP. RESACT holds
C       the residuals of the active constraints, which may be positive.
C       D is the search direction of each line search. DW is either another
C       search direction or the change in gradient along D. The length of W
C       must be at least MAX[M,2*N].
C
C     Set some numbers for the conjugate gradient iterations.
C
      HALF=0.5D0
      ONE=1.0D0
      TINY=1.0D-60
      ZERO=0.0D0
      CTEST=0.01D0
      SNSQ=SNORM*SNORM
C
C     Set the initial elements of RESNEW, RESACT and STEP.
C
      IF (M .GT. 0) THEN
          DO 10 J=1,M
          RESNEW(J)=RESCON(J)
          IF (RESCON(J) .GE. SNORM) THEN
              RESNEW(J)=-ONE
          ELSE IF (RESCON(J) .GE. ZERO) THEN
              RESNEW(J)=DMAX1(RESNEW(J),TINY)
          END IF
   10     CONTINUE
          IF (NACT .GT. 0) THEN
              DO 20 K=1,NACT
              RESACT(K)=RESCON(IACT(K))
   20         RESNEW(IACT(K))=ZERO
          END IF
      END IF
      DO 30 I=1,N
   30 STEP(I)=ZERO
      SS=ZERO
      REDUCT=ZERO
      NCALL=0
C
C     GETACT picks the active set for the current STEP. It also sets DW to
C       the vector closest to -G that is orthogonal to the normals of the
C       active constraints. DW is scaled to have length 0.2*SNORM, as then
C       a move of DW from STEP is allowed by the linear constraints.
C
   40 NCALL=NCALL+1
      CALL GETACT (N,M,AMAT,B,NACT,IACT,QFAC,RFAC,SNORM,RESNEW,
     1  RESACT,G,DW,W,W(N+1))
      IF (W(N+1) .EQ. ZERO) GOTO 320
      SCALE=0.2D0*SNORM/DSQRT(W(N+1))
      DO 50 I=1,N
   50 DW(I)=SCALE*DW(I)
C
C     If the modulus of the residual of an active constraint is substantial,
C       then set D to the shortest move from STEP to the boundaries of the
C       active constraints.
C
      RESMAX=ZERO
      IF (NACT .GT. 0) THEN
          DO 60 K=1,NACT
   60     RESMAX=DMAX1(RESMAX,RESACT(K))
      END IF
      GAMMA=ZERO
      IF (RESMAX .GT. 1.0D-4*SNORM) THEN
          IR=0
          DO 80 K=1,NACT
          TEMP=RESACT(K)
          IF (K .GE. 2) THEN
              DO 70 I=1,K-1
              IR=IR+1
   70         TEMP=TEMP-RFAC(IR)*W(I)
          END IF
          IR=IR+1
   80     W(K)=TEMP/RFAC(IR)
          DO 90 I=1,N
          D(I)=ZERO
          DO 90 K=1,NACT
   90     D(I)=D(I)+W(K)*QFAC(I,K)
C
C     The vector D that has just been calculated is also the shortest move
C       from STEP+DW to the boundaries of the active constraints. Set GAMMA
C       to the greatest steplength of this move that satisfies the trust
C       region bound.
C
          RHS=SNSQ
          DS=ZERO
          DD=ZERO
          DO 100 I=1,N
          SUM=STEP(I)+DW(I)
          RHS=RHS-SUM*SUM
          DS=DS+D(I)*SUM
  100     DD=DD+D(I)**2
          IF (RHS .GT. ZERO) THEN
              TEMP=DSQRT(DS*DS+DD*RHS)
              IF (DS .LE. ZERO) THEN
                  GAMMA=(TEMP-DS)/DD
              ELSE
                  GAMMA=RHS/(TEMP+DS)
              END IF
          END IF
C
C     Reduce the steplength GAMMA if necessary so that the move along D
C       also satisfies the linear constraints.
C
          J=0
  110     IF (GAMMA .GT. ZERO) THEN
              J=J+1
              IF (RESNEW(J) .GT. ZERO) THEN
                  AD=ZERO
                  ADW=ZERO
                  DO 120 I=1,N
                  AD=AD+AMAT(I,J)*D(I)
  120             ADW=ADW+AMAT(I,J)*DW(I)
                  IF (AD .GT. ZERO) THEN
                      TEMP=DMAX1((RESNEW(J)-ADW)/AD,ZERO)
                      GAMMA=DMIN1(GAMMA,TEMP)
                  END IF
              END IF
              IF (J .LT. M) GOTO 110
          END IF
          GAMMA=DMIN1(GAMMA,ONE)
      END IF
C
C     Set the next direction for seeking a reduction in the model function
C       subject to the trust region bound and the linear constraints.
C
      IF (GAMMA .LE. ZERO) THEN
          DO 130 I=1,N
  130     D(I)=DW(I)
          ICOUNT=NACT
      ELSE
          DO 140 I=1,N
  140     D(I)=DW(I)+GAMMA*D(I)
          ICOUNT=NACT-1
      END IF
      ALPBD=ONE
C
C     Set ALPHA to the steplength from STEP along D to the trust region
C       boundary. Return if the first derivative term of this step is
C       sufficiently small or if no further progress is possible.
C
  150 ICOUNT=ICOUNT+1
      RHS=SNSQ-SS
      IF (RHS .LE. ZERO) GOTO 320
      DG=ZERO
      DS=ZERO
      DD=ZERO
      DO 160 I=1,N
      DG=DG+D(I)*G(I)
      DS=DS+D(I)*STEP(I)
  160 DD=DD+D(I)**2
      IF (DG .GE. ZERO) GOTO 320
      TEMP=DSQRT(RHS*DD+DS*DS)
      IF (DS .LE. ZERO) THEN
          ALPHA=(TEMP-DS)/DD
      ELSE
          ALPHA=RHS/(TEMP+DS)
      END IF
      IF (-ALPHA*DG .LE. CTEST*REDUCT) GOTO 320
C
C     Set DW to the change in gradient along D.
C
      IH=0
      DO 170 J=1,N
      DW(J)=ZERO
      DO 170 I=1,J
      IH=IH+1
      IF (I .LT. J) DW(J)=DW(J)+HQ(IH)*D(I)
  170 DW(I)=DW(I)+HQ(IH)*D(J)
      DO 190 K=1,NPT
      TEMP=ZERO
      DO 180 J=1,N
  180 TEMP=TEMP+XPT(K,J)*D(J)
      TEMP=PQ(K)*TEMP
      DO 190 I=1,N
  190 DW(I)=DW(I)+TEMP*XPT(K,I)
C
C     Set DGD to the curvature of the model along D. Then reduce ALPHA if
C       necessary to the value that minimizes the model.
C
      DGD=ZERO
      DO 200 I=1,N
  200 DGD=DGD+D(I)*DW(I)
      ALPHT=ALPHA
      IF (DG+ALPHA*DGD .GT. ZERO) THEN
          ALPHA=-DG/DGD
      END IF
C
C     Make a further reduction in ALPHA if necessary to preserve feasibility,
C       and put some scalar products of D with constraint gradients in W.
C
      ALPHM=ALPHA
      JSAV=0
      IF (M .GT. 0) THEN
          DO 220 J=1,M
          AD=ZERO
          IF (RESNEW(J) .GT. ZERO) THEN
              DO 210 I=1,N
  210         AD=AD+AMAT(I,J)*D(I)
              IF (ALPHA*AD .GT. RESNEW(J)) THEN
                  ALPHA=RESNEW(J)/AD
                  JSAV=J
              END IF
          END IF
  220     W(J)=AD
      END IF
      ALPHA=DMAX1(ALPHA,ALPBD)
      ALPHA=DMIN1(ALPHA,ALPHM)
      IF (ICOUNT .EQ. NACT) ALPHA=DMIN1(ALPHA,ONE)
C
C     Update STEP, G, RESNEW, RESACT and REDUCT.
C
      SS=ZERO
      DO 230 I=1,N
      STEP(I)=STEP(I)+ALPHA*D(I)
      SS=SS+STEP(I)**2
  230 G(I)=G(I)+ALPHA*DW(I)
      IF (M .GT. 0) THEN
          DO 240 J=1,M
          IF (RESNEW(J) .GT. ZERO) THEN
              RESNEW(J)=DMAX1(RESNEW(J)-ALPHA*W(J),TINY)
          END IF
  240     CONTINUE
      END IF
      IF (ICOUNT .EQ. NACT .AND. NACT .GT. 0) THEN
          DO 250 K=1,NACT
  250     RESACT(K)=(ONE-GAMMA)*RESACT(K)
      END IF
      REDUCT=REDUCT-ALPHA*(DG+HALF*ALPHA*DGD)
C
C     Test for termination. Branch to label 40 if there is a new active
C       constraint and if the distance from STEP to the trust region
C       boundary is at least 0.2*SNORM.
C
      IF (ALPHA .EQ. ALPHT) GOTO 320
      TEMP=-ALPHM*(DG+HALF*ALPHM*DGD)
      IF (TEMP .LE. CTEST*REDUCT) GOTO 320
      IF (JSAV .GT. 0) THEN
          IF (SS .LE. 0.64D0*SNSQ) GOTO 40
          GOTO 320
      END IF
      IF (ICOUNT .EQ. N) GOTO 320
C
C     Calculate the next search direction, which is conjugate to the
C       previous one except in the case ICOUNT=NACT.
C
      IF (NACT .GT. 0) THEN
          DO 260 J=NACT+1,N
          W(J)=ZERO
          DO 260 I=1,N
  260     W(J)=W(J)+G(I)*QFAC(I,J)
          DO 280 I=1,N
          TEMP=ZERO
          DO 270 J=NACT+1,N
  270     TEMP=TEMP+QFAC(I,J)*W(J)
  280     W(N+I)=TEMP
      ELSE
          DO 290 I=1,N
  290     W(N+I)=G(I)
      END IF
      IF (ICOUNT .EQ. NACT) THEN
          BETA=ZERO
      ELSE
          WGD=ZERO
          DO 300 I=1,N
  300     WGD=WGD+W(N+I)*DW(I)
          BETA=WGD/DGD
      END IF
      DO 310 I=1,N
  310 D(I)=-W(N+I)+BETA*D(I)
      ALPBD=ZERO
      GOTO 150
C
C     Return from the subroutine.
C
  320 SNORM=ZERO
      IF (REDUCT .GT. ZERO) SNORM=DSQRT(SS)
      G(1)=ZERO
      IF (NCALL .GT. 1) G(1)=ONE
      RETURN
      END

      SUBROUTINE UPDATE (N,NPT,XPT,BMAT,ZMAT,IDZ,NDIM,SP,STEP,
     1  KOPT,KNEW,VLAG,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),SP(*),STEP(*),
     2  VLAG(*),W(*)
C
C     The arguments N, NPT, XPT, BMAT, ZMAT, IDZ, NDIM ,SP and STEP are
C       identical to the corresponding arguments in SUBROUTINE LINCOB.
C     KOPT is such that XPT(KOPT,.) is the current trust region centre.
C     KNEW on exit is usually positive, and then it is the index of an
C       interpolation point to be moved to the position XPT(KOPT,.)+STEP(.).
C       It is set on entry either to its final value or to 0. In the latter
C       case, the final value of KNEW is chosen to maximize the denominator
C       of the matrix updating formula times a weighting factor.
C     VLAG and W are used for working space, the first NPT+N elements of
C       both of these vectors being required.
C
C     The arrays BMAT and ZMAT with IDZ are updated, the new matrices being
C       the ones that are suitable after the shift of the KNEW-th point to
C       the new position XPT(KOPT,.)+STEP(.). A return with KNEW set to zero
C       occurs if the calculation fails due to a zero denominator in the
C       updating formula, which should never happen.
C
C     Set some constants.
C
      HALF=0.5D0
      ONE=1.0D0
      ZERO=0.0D0
      NPTM=NPT-N-1
C
C     Calculate VLAG and BETA for the current choice of STEP. The first NPT
C       elements of VLAG are set to the values of the Lagrange functions at
C       XPT(KOPT,.)+STEP(.). The first NPT components of W_check are held
C       in W, where W_check is defined in a paper on the updating method.
C
      DO 20 K=1,NPT
      W(K)=SP(NPT+K)*(HALF*SP(NPT+K)+SP(K))
      SUM=ZERO
      DO 10 J=1,N
   10 SUM=SUM+BMAT(K,J)*STEP(J)
   20 VLAG(K)=SUM
      BETA=ZERO
      DO 40 K=1,NPTM
      SUM=ZERO
      DO 30 I=1,NPT
   30 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 40 I=1,NPT
   40 VLAG(I)=VLAG(I)+SUM*ZMAT(I,K)
      BSUM=ZERO
      DX=ZERO
      SSQ=ZERO
      DO 70 J=1,N
      SUM=ZERO
      DO 50 I=1,NPT
   50 SUM=SUM+W(I)*BMAT(I,J)
      BSUM=BSUM+SUM*STEP(J)
      JP=NPT+J
      DO 60 K=1,N
   60 SUM=SUM+BMAT(JP,K)*STEP(K)
      VLAG(JP)=SUM
      BSUM=BSUM+SUM*STEP(J)
      DX=DX+STEP(J)*XPT(KOPT,J)
   70 SSQ=SSQ+STEP(J)**2
      BETA=DX*DX+SSQ*(SP(KOPT)+DX+DX+HALF*SSQ)+BETA-BSUM
      VLAG(KOPT)=VLAG(KOPT)+ONE
C
C     If KNEW is zero initially, then pick the index of the interpolation
C       point to be deleted, by maximizing the absolute value of the
C       denominator of the updating formula times a weighting factor.
C       
C
      IF (KNEW .EQ. 0) THEN
          DENMAX=ZERO
          DO 100 K=1,NPT
          HDIAG=ZERO
          DO 80 J=1,NPTM
          TEMP=ONE
          IF (J .LT. IDZ) TEMP=-ONE
   80     HDIAG=HDIAG+TEMP*ZMAT(K,J)**2
          DENABS=DABS(BETA*HDIAG+VLAG(K)**2)
          DISTSQ=ZERO
          DO 90 J=1,N
   90     DISTSQ=DISTSQ+(XPT(K,J)-XPT(KOPT,J))**2
          TEMP=DENABS*DISTSQ*DISTSQ
          IF (TEMP .GT. DENMAX) THEN
              DENMAX=TEMP
              KNEW=K
          END IF
  100     CONTINUE
      END IF
C
C     Apply the rotations that put zeros in the KNEW-th row of ZMAT.
C
      JL=1
      IF (NPTM .GE. 2) THEN
          DO 120 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 110 I=1,NPT
              TEMP=TEMPA*ZMAT(I,JL)+TEMPB*ZMAT(I,J)
              ZMAT(I,J)=TEMPA*ZMAT(I,J)-TEMPB*ZMAT(I,JL)
  110         ZMAT(I,JL)=TEMP
              ZMAT(KNEW,J)=ZERO
          END IF
  120     CONTINUE
      END IF
C
C     Put the first NPT components of the KNEW-th column of the Z Z^T matrix
C       into W, and calculate the parameters of the updating formula.
C
      TEMPA=ZMAT(KNEW,1)
      IF (IDZ .GE. 2) TEMPA=-TEMPA
      IF (JL .GT. 1) TEMPB=ZMAT(KNEW,JL)
      DO 130 I=1,NPT
      W(I)=TEMPA*ZMAT(I,1)
      IF (JL .GT. 1) W(I)=W(I)+TEMPB*ZMAT(I,JL)
  130 CONTINUE
      ALPHA=W(KNEW)
      TAU=VLAG(KNEW)
      TAUSQ=TAU*TAU
      DENOM=ALPHA*BETA+TAUSQ
      VLAG(KNEW)=VLAG(KNEW)-ONE
      IF (DENOM .EQ. ZERO) THEN
          KNEW=0
          GOTO 180
      END IF
      SQRTDN=DSQRT(DABS(DENOM))
C
C     Complete the updating of ZMAT when there is only one nonzero element
C       in the KNEW-th row of the new matrix ZMAT. IFLAG is set to one when
C       the value of IDZ is going to be reduced.
C
      IFLAG=0
      IF (JL .EQ. 1) THEN
          TEMPA=TAU/SQRTDN
          TEMPB=ZMAT(KNEW,1)/SQRTDN
          DO 140 I=1,NPT
  140     ZMAT(I,1)=TEMPA*ZMAT(I,1)-TEMPB*VLAG(I)
          IF (DENOM .LT. ZERO) THEN
              IF (IDZ .EQ. 1) THEN
                  IDZ=2
              ELSE
                  IFLAG=1
              END IF
          END IF
      ELSE
C
C     Complete the updating of ZMAT in the alternative case.
C
          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*SQRTDN
          DO 150 I=1,NPT
          ZMAT(I,JA)=SCALA*(TAU*ZMAT(I,JA)-TEMP*VLAG(I))
  150     ZMAT(I,JB)=SCALB*(ZMAT(I,JB)-TEMPA*W(I)-TEMPB*VLAG(I))
          IF (DENOM .LE. ZERO) THEN
              IF (BETA .LT. ZERO) THEN
                  IDZ=IDZ+1
              ELSE
                  IFLAG=1
              END IF
          END IF
      END IF
C
C     Reduce IDZ when the diagonal part of the ZMAT times Diag(DZ) times
C       ZMAT^T factorization gains another positive element. Then exchange
C       the first and IDZ-th columns of ZMAT.
C
      IF (IFLAG .EQ. 1) THEN
          IDZ=IDZ-1
          DO 160 I=1,NPT
          TEMP=ZMAT(I,1)
          ZMAT(I,1)=ZMAT(I,IDZ)
  160     ZMAT(I,IDZ)=TEMP
      END IF
C
C     Finally, update the matrix BMAT.
C
      DO 170 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 170 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)
  170 CONTINUE
  180 RETURN
      END