SUBROUTINE QPCOLR( NOCURV, POSDEF, RENEWR, UNITQ, QPHESS,
     *                   N, NCOLR, NCTOTL, NFREE, NQ, NROWH, NCOLH,
     *                   NROWRT, NCOLRT, NHESS, KFREE,
     *                   CSLAST, SNLAST, DRMAX, EMAX, HSIZE, RDLAST,
     *                   HESS, RT, SCALE, ZY, HZ, WRK )
C
C     IMPLICIT           REAL*8(A-H,O-Z)
      INTEGER            N, NCOLR, NCTOTL, NFREE, NQ, NROWH, NCOLH,
     *                   NROWRT, NCOLRT, NHESS
      INTEGER            KFREE(N)
      DOUBLE PRECISION   CSLAST, SNLAST, DRMAX, EMAX, HSIZE, RDLAST
      DOUBLE PRECISION   HESS(NROWH,NCOLH), RT(NROWRT,NCOLRT), HZ(N),
     *                   SCALE(NCTOTL), ZY(NQ,NQ)
      DOUBLE PRECISION   WRK(N)
      LOGICAL            NOCURV, POSDEF, RENEWR, UNITQ
      EXTERNAL           QPHESS
C
      INTEGER            NOUT, MSG, ISTART
      DOUBLE PRECISION   WMACH
      COMMON    /SOLMCH/ WMACH(15)
      COMMON    /SOL1CM/ NOUT, MSG, ISTART
C
      LOGICAL            SCLDQP
      COMMON    /SOL2LP/ SCLDQP
C
C  *********************************************************************
C  QPCOLR  IS USED TO COMPUTE ELEMENTS OF THE  (NCOLR)-TH  COLUMN OF  R,
C  THE CHOLESKY FACTOR OF THE PROJECTED HESSIAN.  IF  RENEWR  IS  TRUE
C  ON ENTRY,  THE COMPLETE COLUMN IS TO BE COMPUTED.  OTHERWISE, ONLY
C  THE LAST DIAGONAL ELEMENT IS REQUIRED.
C  IF THE RESULTING PROJECTED HESSIAN IS SINGULAR OR INDEFINITE, ITS
C  LAST DIAGONAL ELEMENT IS INCREASED BY AN AMOUNT  EMAX  THAT ENSURES
C  POSITIVE DEFINITENESS.  THIS DIAGONAL MODIFICATION WILL ALTER THE
C  SCALE OF THE QP SEARCH VECTOR  P, BUT NOT ITS DIRECTION.
C
C  ON EXIT,  QPCOLR WILL HAVE STORED THE  NCOLR  ELEMENTS OF THE NEW
C  COLUMN OF  R  IN THE ARRAY  RT,  AND SET THE VARIABLES  NOCURV,
C  POSDEF,  RENEWR,  DRMAX,  EMAX  AND  HSIZE.
C
C  SYSTEMS OPTIMIZATION LABORATORY, STANFORD UNIVERSITY.
C  ORIGINAL VERSION MARCH 1982.  REV. APRIL 1982.
C  *********************************************************************
C
      INTEGER            J, JTHCOL, K, NCOLR1
      DOUBLE PRECISION   EPSMCH, RDSMAX, RDSMIN, RDSQ, RNORM, S, ZTHZ
      DOUBLE PRECISION   ZERO, ONE, TWO, TEN
      DOUBLE PRECISION   DSQRT, V2NORM
      DOUBLE PRECISION   DABS, DMAX1
      DATA               ZERO  , ONE   , TWO   , TEN
     *                  /0.0D+0, 1.0D+0, 2.0D+0, 10.0D+0/
C
      EPSMCH = WMACH(3)
C
      IF (RENEWR) GO TO 300
C
C  ---------------------------------------------------------------------
C  ONLY THE LAST ELEMENT OF THE NEW COLUMN OF  R  NEED BE COMPUTED.
C  THIS SITUATION CAN ONLY OCCUR WHEN A CONSTRAINT IS ADDED TO THE
C  WORKING SET WITH  ZTHZ  NOT POSITIVE DEFINITE.
C  ---------------------------------------------------------------------
C  THE LAST DIAGONAL ELEMENT OF  R  IS THAT OF  ZTHZ  PLUS A DIAGONAL
C  MODIFICATION.  THE SQUARE OF THE TRUE DIAGONAL IS RECOVERED FROM THE
C  ROTATIONS USED TO UPDATE  R  WHEN THE CONSTRAINT WAS ADDED TO THE
C  WORKING SET.
C
      RDLAST = RT(NCOLR,NCOLR)
      S      = DABS( SNLAST )
      RDSQ   = ( (CSLAST - S)*RDLAST )*( (CSLAST + S)*RDLAST )
      GO TO 600
C
C  ---------------------------------------------------------------------
C  THE PROJECTED HESSIAN IS EXPANDED BY A ROW AND COLUMN.  COMPUTE THE
C  FIRST  (NCOLR - 1)  ELEMENTS OF THE NEW COLUMN OF THE CHOLESKY FACTOR
C  R.  ALSO, COMPUTE  RDSQ,  THE SQUARE OF THE LAST DIAGONAL ELEMENT.
C  ---------------------------------------------------------------------
  300 CALL ZEROVC( N, WRK, N, 1 )
      IF (UNITQ) GO TO 320
C
C  EXPAND THE NEW COLUMN OF Z IN TO AN N-VECTOR.
C
      DO 310 K = 1, NFREE
         J      = KFREE(K)
         WRK(J) = ZY(K,NCOLR)
  310 CONTINUE
      IF (SCLDQP) CALL DSCALE( N, SCALE, N, 1, WRK, N, 1 )
      JTHCOL = 0
      GO TO 330
C
C  ONLY BOUNDS ARE IN THE WORKING SET (NFREE  IS EQUAL TO  NCOLZ).  THE
C  (NCOLR)-TH  COLUMN OF  Z  IS JUST A COLUMN OF THE IDENTITY MATRIX.
C
  320 JTHCOL      = KFREE(NCOLR)
      WRK(JTHCOL) = ONE
C
C  COMPUTE THE HESSIAN TIMES THE LAST COLUMN OF Z.
C
  330 CALL QPHESS( N, NROWH, NCOLH, JTHCOL, HESS, WRK, HZ )
      NHESS  = NHESS + 1
C
      IF (UNITQ  .AND.  SCLDQP) CALL SSCALE( N, SCALE(JTHCOL), HZ, N,1 )
      IF (              SCLDQP) CALL DSCALE( N, SCALE, N, 1,   HZ, N,1 )
C
C  COMPUTE THE  (NCOLR)-TH  COLUMN OF  Z(T)H Z.
C
      CALL ZYPROD( 4, N, NFREE, NCOLR, NFREE, NQ, UNITQ,
     *             KFREE, KFREE, HZ, ZY, WRK )
C
      CALL COPYVC( NCOLR, HZ, NCOLR, 1, RT(1,NCOLR), NCOLR, 1 )
C
C  COMPUTE THE FIRST  (NCOLR - 1)  ELEMENTS OF THE LAST COLUMN OF  R.
C
      NCOLR1 = NCOLR - 1
      ZTHZ   = RT(NCOLR,NCOLR)
      RDSQ   = ZTHZ
      IF (NCOLR1 .EQ. 0) GO TO 370
      CALL RSOLVE( 2, NROWRT, NCOLR1, RT, RT(1,NCOLR) )
      RNORM  = V2NORM( NCOLR1, RT(1,NCOLR), NCOLR1, 1 )
C
C  COMPUTE THE SQUARE OF THE LAST DIAGONAL ELEMENT OF  R.
C
      RDSQ   = ZTHZ - RNORM*RNORM
C
C  UPDATE THE ESTIMATE OF THE NORM OF THE HESSIAN.
C
  370 HSIZE  = DMAX1( HSIZE, DABS( ZTHZ ) )
C
C  ---------------------------------------------------------------------
C  COMPUTE  RDLAST, THE LAST DIAGONAL OF  R.  THE VARIABLES POSDEF  AND
C  NOCURV  ARE SET HERE.  THEY ARE USED TO INDICATE IF THE NEW PROJECTED
C  HESSIAN IS POSITIVE DEFINITE OR SINGULAR.  IF  POSDEF  IS SET TO
C  FALSE,  RDLAST  WILL BE THAT OF  ZTHZ  PLUS A DIAGONAL MODIFICATION.
C  IF THE REQUIRED DIAGONAL MODIFICATION IS LARGE,  RENEWR  WILL BE SET
C  TO BE  TRUE,  INDICATING THAT THE LAST ROW AND COLUMN OF  R  MUST BE
C  RECOMPUTED WHEN A CONSTRAINT IS ADDED TO THE WORKING SET DURING THE
C  NEXT ITERATION.
C  ---------------------------------------------------------------------
  600 NOCURV = .FALSE.
      RENEWR = .FALSE.
      EMAX   = ZERO
C
C  RDSMIN  IS THE SQUARE OF THE SMALLEST ALLOWABLE DIAGONAL ELEMENT
C  FOR A POSITIVE-DEFINITE CHOLESKY FACTOR.  NOTE THAT THE TEST FOR A
C  SINGULAR MATRIX IS SCALE DEPENDENT.
C
      IF (NCOLR .EQ. 1) RDSMIN =  EPSMCH*HSIZE
      IF (NCOLR .GT. 1) RDSMIN = (EPSMCH*DRMAX) * DRMAX
      POSDEF = RDSQ .GT. RDSMIN
      IF (POSDEF) GO TO 900
C
      IF (RDSQ .LT. ( - RDSMIN )) GO TO 610
C
C  ---------------------------------------------------------------------
C  THE PROJECTED HESSIAN IS SINGULAR.
C  ---------------------------------------------------------------------
C  THE QUADRATIC HAS NO CURVATURE ALONG AT LEAST ONE DIRECTION.  THE
C  PERTURBATION  EMAX  IS CHOSEN TO MAKE THE NEW EIGENVALUE OF  ZTHZ
C  SMALL AND POSITIVE.
C
      EMAX   = RDSMIN - RDSQ
      RDSQ   = RDSMIN
      NOCURV = .TRUE.
      GO TO 900
C
C  ---------------------------------------------------------------------
C  THE PROJECTED HESSIAN IS INDEFINITE.  THERE ARE TWO CASES.
C  ---------------------------------------------------------------------
C  CASE 1.  THE MODULUS OF THE NEW LAST DIAGONAL OF  R  IS NOT TOO
C  LARGE.  THE MODULUS OF  RDSQ  IS USED FOR THE SQUARE ROOT.
C
  610 RDSMAX  = TEN * HSIZE
      IF (RDSQ .LT. ( - RDSMAX )) GO TO 620
      EMAX    = - TWO*RDSQ
      GO TO 900
C
C  CASE 2.  THE MODULUS OF THE LAST DIAGONAL OF  R  IS JUDGED TO BE TOO
C  LARGE (SOME LOSS OF PRECISION MAY HAVE OCCURRED).  SET  RENEWR  SO
C  THAT THE LAST COLUMN IS RECOMPUTED LATER.
C
  620 EMAX   = RDSMAX - RDSQ
      RENEWR = .TRUE.
      RDSQ   = RDSMAX
C
C  COMPUTE THE LAST DIAGONAL ELEMENT.
C
  900 RDLAST = DSQRT( DABS( RDSQ ) )
      RT(NCOLR,NCOLR) = RDLAST
C
      IF (MSG .GE. 80  .AND.  (.NOT. POSDEF))
     *WRITE (NOUT, 9000) POSDEF, NOCURV, EMAX, RDLAST
C
      RETURN
C
 9000 FORMAT(/ 54H //QPCOLR//  POSDEF NOCURV          EMAX        RDLAST
     *       / 13H //QPCOLR//  , L6, 1X, L6, 2(1PE14.4) )
C
C  END OF QPCOLR
      END