C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     FILE CMSUBS66 FORTRAN
C
C     ADDCON   ALLOC    BDPERT   BNDALF   CHKDAT   DELCON
C     FINDP    GETLAM   PRTSOL   RSOLVE   TQADD    TSOLVE   ZYPROD
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE ADDCON( MODFYG, MODFYR, ORTHOG, UNITQ, INFORM,
     *                   IFIX, IADD, JADD, NACTIV, NCOLR, NCOLZ, NFREE,
     *                   N, NQ, NROWA, NROWRT, NCOLRT, KFREE,
     *                   CONDMX, CSLAST, SNLAST,
     *                   A, QTG, RT, ZY, WRK1, WRK2 )
C
C     IMPLICIT           REAL*8(A-H,O-Z)
      LOGICAL            MODFYG, MODFYR, ORTHOG, UNITQ
      INTEGER            INFORM, IFIX, IADD, JADD, NACTIV, NCOLR,
     *                   NCOLZ, NFREE, N, NQ, NROWA, NROWRT, NCOLRT
      INTEGER            KFREE(N)
      DOUBLE PRECISION   CONDMX, CSLAST, SNLAST
      DOUBLE PRECISION   A(NROWA,N), QTG(N), RT(NROWRT,NCOLRT),
     *                   ZY(NQ,NQ), WRK1(N), WRK2(N)
C
      INTEGER            NOUT, MSG, ISTART
      DOUBLE PRECISION   ASIZE, DTMAX, DTMIN
      DOUBLE PRECISION   WMACH
      COMMON    /SOLMCH/ WMACH(15)
      COMMON    /SOL1CM/ NOUT, MSG, ISTART
      COMMON    /SOL5CM/ ASIZE, DTMAX, DTMIN
C
C  *********************************************************************
C  ADDCON  UPDATES THE FACTORIZATION OF THE MATRIX OF
C  CONSTRAINTS IN THE WORKING SET,  A(FREE) * (Z Y) = (0 T).
C  IF THE LOGICAL VARIABLE  MODFYR  IS TRUE, THE CHOLESKY FACTORIZATION
C  OF THE PROJECTED HESSIAN, R(T)*R, IS UPDATED ALSO.
C
C  THERE ARE THREE SEPARATE CASES TO CONSIDER (ALTHOUGH EACH CASE
C  SHARES CODE WITH ANOTHER)...
C
C  (1) A FREE VARIABLE BECOMES FIXED ON ONE OF ITS BOUNDS WHEN THERE
C      ARE ALREADY SOME GENERAL CONSTRAINTS IN THE WORKING SET.
C
C  (2) A FREE VARIABLE BECOMES FIXED ON ONE OF ITS BOUNDS WHEN THERE
C      ARE ONLY BOUND CONSTRAINTS IN THE WORKING SET.
C
C  (3) A GENERAL CONSTRAINT (CORRESPONDING TO ROW  IADD  OF  A) IS
C      ADDED TO THE WORKING SET.
C
C      IN CASES (1) AND (2), WE ASSUME THAT  KFREE(IFIX) = JADD.
C  IN ALL CASES,  JADD  IS THE INDEX OF THE CONSTRAINT BEING ADDED.
C
C  IF THERE ARE NO GENERAL CONSTRAINTS IN THE WORKING SET,  THE
C  MATRIX  Q = (Z Y)  IS THE IDENTITY AND WILL NOT BE TOUCHED.
C
C  IF  MODFYR  IS TRUE AND  NCOLZ IS GREATER THAN ONE ON ENTRY,
C  CSLAST AND SNLAST CONTAIN THE LAST OF THE SEQUENCE OF GIVENS
C  ROTATIONS USED TO REDUCE THE INTERMEDIATE UPPER-HESSENBERG MATRIX
C  TO UPPER-TRIANGULAR FORM.  THESE ELEMENTS ARE NEEDED BY QPCORE.
C
C  IF  MODFYG  IS TRUE ON ENTRY, THE COLUMN TRANSFORMATIONS ARE
C  APPLIED TO THE VECTOR  Q(T)GRAD,  STORED IN  QTG.
C
C  SYSTEMS OPTIMIZATION LABORATORY, STANFORD UNIVERSITY.
C  VERSION OF JANUARY 1982.  REV. OCT. 1982.  MARCH 1983.
C  MARCH 1983.  HOUSEHOLDER REFLECTION USED TO ADD GENERAL CONSTRAINT.
C  APRIL 1983.  ELIMINATIONS ADDED AS AN OPTION.
C  *********************************************************************
C
      INTEGER            I, INCT, ISWAP, ITRANS, J, K, KP1, LDIAG, LENQ,
     *                   LENRT, LROWA, LROWR, NACT1, NCOLZ1, NELM,NFREE1
      DOUBLE PRECISION   BETA, COND, CONDBD, CS, D, DELTA, DTNEW,
     *                   EPSMCH, ONE, POINT9, SN, TDTMAX, TDTMIN,
     *                   ZERO
      DOUBLE PRECISION   DOT, QUOTNT, V2NORM
      DOUBLE PRECISION   DMAX1, DMIN1
      DATA               ZERO  , POINT9, ONE
     *                  /0.0D+0, 0.9D+0, 1.0D+0/
C
      EPSMCH = WMACH(3)
C
C  IF THE CONDITION ESTIMATOR OF THE UPDATED FACTORS IS GREATER THAN
C  CONDBD,  A WARNING MESSAGE IS PRINTED.
C
      CONDBD = EPSMCH**(- POINT9)
      LENQ   = NQ   *(NQ - 1) + 1
      LROWA  = NROWA*(N  - 1) + 1
      NCOLZ1 = NCOLZ - 1
      IF (JADD .GT. N) GO TO 200
C
C  ---------------------------------------------------------------------
C  A SIMPLE BOUND HAS ENTERED THE WORKING SET.  IADD  IS NOT USED.
C  ---------------------------------------------------------------------
      IF (MSG .GE. 80)
     *WRITE (NOUT, 1010) NACTIV, NCOLZ, NFREE, IFIX, JADD, UNITQ
C
C  SET  WRK1 = APPROPRIATE ROW OF  Q.
C  REORDER THE ELEMENTS OF  KFREE (THIS REQUIRES REORDERING THE
C  CORRESPONDING ROWS OF  Q).
C
      NFREE1 = NFREE - 1
      NACT1  = NACTIV
      IF (UNITQ) GO TO 120
C
C  Q  IS STORED EXPLICITLY.  INTERCHANGE COMPONENTS  IFIX  AND  NFREE
C  OF  KFREE  AND SWAP THE CORRESPONDING ROWS OF  Q.
C
      CALL COPYVC( NFREE, ZY(IFIX,1), LENQ, NQ, WRK1, N, 1 )
      IF (IFIX .EQ. NFREE) GO TO 400
      KFREE(IFIX) = KFREE(NFREE)
      CALL COPYVC( NFREE, ZY(NFREE,1), LENQ, NQ, ZY(IFIX,1), LENQ, NQ )
      GO TO 400
C
C  Q  IS NOT STORED, BUT  KFREE  DEFINES AN ORDERING OF THE COLUMNS
C  OF THE IDENTITY MATRIX THAT IMPLICITLY DEFINE  Z.
C  REORDER  KFREE  SO THAT VARIABLES  IFIX+1,...,NFREE  ARE MOVED ONE
C  POSITION TO THE LEFT.
C
  120 CALL ZEROVC( NFREE, WRK1, N, 1 )
      WRK1(IFIX) = ONE
      IF (IFIX .EQ. NFREE) GO TO 400
      DO 130 I = IFIX, NFREE1
         KFREE(I) = KFREE(I+1)
  130 CONTINUE
      GO TO 400
C
C  ---------------------------------------------------------------------
C  A GENERAL CONSTRAINT HAS ENTERED THE WORKING SET.  IFIX IS NOT USED.
C  ---------------------------------------------------------------------
  200 IF (MSG .GE. 80)
     *WRITE (NOUT, 1020) NACTIV, NCOLZ, NFREE, IADD, JADD, UNITQ
C
      NACT1  = NACTIV + 1
C
C  TRANSFORM THE INCOMING ROW OF  A  BY  Q(T).
C
      CALL COPYVC( N, A(IADD,1), LROWA, NROWA, WRK1, N, 1 )
      CALL ZYPROD( 8, N, NACTIV, NCOLZ, NFREE, NQ, UNITQ,
     *             KFREE, KFREE, WRK1, ZY, WRK2 )
C
      IF (.NOT. UNITQ) GO TO 250
C
C  THIS IS THE FIRST GENERAL CONSTRAINT TO BE ADDED  --  SET  Q = I.
C
      DO 220 J = 1, NFREE
         CALL ZEROVC( NFREE, ZY(1,J), NFREE, 1 )
         ZY(J,J) = ONE
  220 CONTINUE
      UNITQ  = .FALSE.
C
C  CHECK THAT THE INCOMING ROW IS NOT DEPENDENT UPON THOSE
C  ALREADY IN THE WORKING SET.
C
  250 DTNEW  = V2NORM( NCOLZ, WRK1, NCOLZ, 1 )
      IF (NACT1 .GT. 1) GO TO 300
C
C  THIS IS THE ONLY GENERAL CONSTRAINT IN THE WORKING SET.
C
      COND   = QUOTNT( ASIZE, DTNEW )
      IF (COND .GE. CONDMX) GO TO 910
      IF (COND .GE. CONDBD) WRITE (NOUT, 2000) JADD
      DTMAX  = DTNEW
      DTMIN  = DTNEW
      GO TO 400
C
C  THERE ARE ALREADY SOME GENERAL CONSTRAINTS IN THE WORKING SET.
C  UPDATE THE ESTIMATE OF THE CONDITION NUMBER.
C
  300 TDTMAX = DMAX1( DTNEW, DTMAX )
      TDTMIN = DMIN1( DTNEW, DTMIN )
      COND   = QUOTNT( TDTMAX, TDTMIN )
      IF (COND .GE. CONDMX) GO TO 910
      IF (COND .GE. CONDBD) WRITE (NOUT, 2000) JADD
      DTMAX  = TDTMAX
      DTMIN  = TDTMIN
C
C  ---------------------------------------------------------------------
C  USE ONE OR MORE COLUMN TRANSFORMATIONS TO REDUCE THE FIRST  NCOLZ1
C  ELEMENTS OF  WRK1  TO ZERO.  THIS AFFECTS  ZY,  EXCEPT IF (UNITQ).
C  THE TRANSFORMATIONS MAY ALSO BE APPLIED TO  QTG  AND  R.
C  ---------------------------------------------------------------------
C
  400 IF (NCOLZ1 .EQ. 0    ) GO TO 600
      IF (MODFYR .OR. UNITQ) GO TO 500
C
C  ---------------------------------------------------------------------
C  THERE IS NO  R.  USE A SINGLE ELIMINATION OR HOUSEHOLDER MATRIX.
C  ---------------------------------------------------------------------
      IF (ORTHOG) GO TO 440
C
C  *********************************************************************
C  ELIMINATION.
C  WE USE   ELM( ..., ZERO, ZERO )   TO PERFORM AN INTERCHANGE.
C  *********************************************************************
      CALL ETAGEN( NCOLZ1, WRK1(NCOLZ), WRK1, NCOLZ1, 1, ISWAP, ITRANS )
      IF (ISWAP .GT. 0)
     *   CALL ELM   ( ORTHOG, NFREE, ZY(1,NCOLZ), NFREE, 1,
     *                               ZY(1,ISWAP), NFREE, 1, ZERO, ZERO )
C
      IF (ITRANS .EQ. 0) GO TO 420
C
      DO 410 J = 1, NCOLZ1
         D     = WRK1(J)
         IF (D .EQ. ZERO) GO TO 410
         CALL AXPY( NFREE, D, ZY(1,NCOLZ), NFREE, 1, ZY(1,J), NFREE, 1 )
  410 CONTINUE
C
  420 IF (.NOT. MODFYG) GO TO 600
      IF (ISWAP .GT. 0)
     *   CALL ELM   ( ORTHOG, 1, QTG(NCOLZ), 1, 1,
     *                           QTG(ISWAP), 1, 1, ZERO, ZERO )
C
      IF (ITRANS .GT. 0)
     *   CALL AXPY( NCOLZ1, QTG(NCOLZ), WRK1, NCOLZ1, 1,
     *                                  QTG , NCOLZ1, 1 )
      GO TO 600
C
C  *********************************************************************
C  ORTHOGONAL TRANSFORMATION.
C  WE USE A HOUSEHOLDER REFLECTION,   I  -  1/BETA  V V(T).
C
C  THERE ARE TWO WAYS OF APPLYING THE REFLECTION.  THE UPDATE TO  Z
C  IS DONE VIA   W  =  Z * V,   Z  =  Z  -  1/BETA  W V(T),
C  WHERE  V = WRK1 (FROM HOUSEHOLDER), AND  W = WRK2 (WORKSPACE).
C
C  THE UPDATE TO  QTG  IS THE MORE USUAL  D =  - QTG(T)*V / BETA,
C  QTG  =  QTG  +  D * V.
C
C  NOTE THAT  DELTA  HAS TO BE STORED AFTER THE REFLECTION IS USED.
C  *********************************************************************
  440 CALL REFGEN( NCOLZ1, WRK1(NCOLZ), WRK1, NCOLZ1, 1, BETA, DELTA )
      IF (BETA .LE. ZERO) GO TO 600
      CALL ZEROVC( NFREE, WRK2, NFREE, 1 )
C
      DO 450 J = 1, NCOLZ
         D     = WRK1(J)
         IF (D .EQ. ZERO) GO TO 450
         CALL AXPY( NFREE, D, ZY(1,J), NFREE, 1, WRK2, NFREE, 1 )
  450 CONTINUE
C
      DO 460 J = 1, NCOLZ
         D     = WRK1(J)
         IF (D .EQ. ZERO) GO TO 460
         D     = - D/BETA
         CALL AXPY( NFREE, D, WRK2, NFREE, 1, ZY(1,J), NFREE, 1 )
  460 CONTINUE
C
      IF (.NOT. MODFYG) GO TO 470
      D      = DOT( NCOLZ, WRK1, NCOLZ, 1, QTG, NCOLZ, 1 )
      D      = - D/BETA
      CALL AXPY( NCOLZ, D, WRK1, NCOLZ, 1, QTG, NCOLZ, 1 )
C
  470 WRK1(NCOLZ) = DELTA
      GO TO 600
C
C  ---------------------------------------------------------------------
C  R  HAS TO BE MODIFIED.  USE A SEQUENCE OF 2*2 TRANSFORMATIONS.
C  ---------------------------------------------------------------------
  500 LROWR  = NCOLR
C
      DO 510 K = 1, NCOLZ1
C
C        COMPUTE THE TRANSFORMATION THAT REDUCES WRK1(K) TO ZERO,
C        THEN APPLY IT TO THE RELEVANT COLUMNS OF  Z  AND  GRAD(T)Q.
C
         KP1    = K + 1
         CALL ELMGEN( ORTHOG, WRK1(KP1), WRK1(K), CS, SN )
         IF (.NOT. UNITQ)
     *   CALL ELM   ( ORTHOG, NFREE, ZY(1,KP1), NFREE, 1,
     *                               ZY(1,K  ), NFREE, 1, CS, SN )
         IF (MODFYG)
     *   CALL ELM   ( ORTHOG, 1, QTG(KP1), 1, 1, QTG(K), 1, 1, CS, SN )
C
C        APPLY THE SAME TRANSFORMATION TO THE COLS OF  R  IF RELEVANT.
C        THIS GENERATES A SUBDIAGONAL ELEMENT IN  R  WHICH MUST BE
C        ELIMINATED BY A ROW ROTATION.  THE LAST SUCH ROW ROTATION
C        IS NEEDED BY  QPCORE.
C
         IF (.NOT. (MODFYR  .AND.  K .LT. NCOLR)) GO TO 510
         RT(KP1,K) = ZERO
         CALL ELM   ( ORTHOG, KP1, RT(1,KP1), KP1, 1,
     *                             RT(1,K  ), KP1, 1, CS, SN )
         CALL ROTGEN( RT(K,K), RT(KP1,K), CSLAST, SNLAST )
         LROWR  = LROWR - 1
         LENRT  = NROWRT*(LROWR - 1) + 1
         CALL ROT3  ( LROWR, RT(K,KP1), LENRT, NROWRT,
     *                     RT(KP1,KP1), LENRT, NROWRT, CSLAST, SNLAST )
  510 CONTINUE
C
C  ---------------------------------------------------------------------
C  IF ADDING A GENERAL CONSTRAINT, INSERT THE NEW ROW OF  T  AND EXIT.
C  ---------------------------------------------------------------------
  600 IF (JADD .LE. N) GO TO 700
      LENRT  = NROWRT*NACTIV + 1
      CALL COPYVC( NACT1, WRK1(NCOLZ), NACT1, 1,
     *                RT(NACT1,NCOLZ), LENRT, NROWRT )
      GO TO 900
C
C  ---------------------------------------------------------------------
C  WE ARE ADDING A BOUND.  CONTINUE REDUCING THE ELEMENTS OF  WRK1
C  TO ZERO.  THIS AFFECTS  Y,  T  AND  QTG.
C  ---------------------------------------------------------------------
C  FIRST, SET THE SUPER-DIAGONAL ELEMENTS OF T TO ZERO.
C
  700 IF (NACTIV .EQ. 0) GO TO 790
      LENRT = NROWRT*(NACTIV - 1) + 1
      CALL ZEROVC( NACTIV, RT(NACTIV,NCOLZ), LENRT, (NROWRT - 1) )
      NELM   = 1
      LDIAG  = NACTIV
C
      DO 710 K = NCOLZ, NFREE1
         CALL ELMGEN( ORTHOG, WRK1(K+1), WRK1(K), CS, SN )
         CALL ELM   ( ORTHOG, NFREE, ZY(1,K+1), NQ, 1,
     *                               ZY(1,K  ), NQ, 1, CS, SN )
         CALL ELM   ( ORTHOG, NELM, RT(LDIAG,K+1), NROWRT, 1,
     *                              RT(LDIAG,K  ), NROWRT, 1, CS, SN )
         IF (MODFYG)
     *   CALL ELM   ( ORTHOG, 1, QTG(K+1), 1, 1, QTG(K), 1, 1, CS, SN )
         NELM  = NELM  + 1
         LDIAG = LDIAG - 1
  710 CONTINUE
C
C  ---------------------------------------------------------------------
C  THE DIAGONALS OF  T  HAVE BEEN ALTERED.  RECOMPUTE THE LARGEST AND
C  SMALLEST VALUES.
C  ---------------------------------------------------------------------
      LENRT  = NROWRT*(NACTIV - 1) + 1
      INCT   = NROWRT - 1
      CALL CONDVC( NACTIV, RT(NACTIV,NCOLZ1+1), LENRT, INCT,
     *             DTMAX, DTMIN )
      IF ((DTMIN/DTMAX)*CONDMX .LT. ONE) GO TO 910
      IF ((DTMIN/DTMAX)*CONDBD .LT. ONE) WRITE (NOUT, 2000) JADD
C
C  THE LAST ROW OF  ZY  HAS BEEN TRANSFORMED TO A MULTIPLE OF THE
C  UNIT VECTOR  E(NFREE).  IF ORTHOGONAL TRANSFORMATIONS HAVE BEEN
C  USED THROUGHOUT, THE LAST COLUMN OF  ZY  IS THE SAME.   WE CAN
C  THEREFORE RESURRECT THE GRADIENT ELEMENT OF THE NEWLY-FIXED VARIABLE.
C
  790 IF (ORTHOG .AND. MODFYG)
     *QTG(NFREE) = QTG(NFREE)/WRK1(NFREE)
C
C  ---------------------------------------------------------------------
C  THE FACTORIZATION HAS BEEN SUCCESSFULLY UPDATED.
C  ---------------------------------------------------------------------
  900 INFORM = 0
      RETURN
C
C  THE PROPOSED WORKING SET APPEARS TO BE LINEARLY DEPENDENT.
C
  910 INFORM = 1
      IF (.NOT. MSG .GE. 80) RETURN
C
      WRITE (NOUT, 3000)
      IF (JADD .LE. N) WRITE (NOUT, 3010) ASIZE, DTMAX, DTMIN
      IF (JADD .GT. N) WRITE (NOUT, 3020) ASIZE, DTMAX, DTMIN, DTNEW
      RETURN
C
 1010 FORMAT(/ 32H //ADDCON//  SIMPLE BOUND ADDED.
     *       / 49H //ADDCON//  NACTIV NCOLZ NFREE  IFIX  JADD UNITQ
     *       / 13H //ADDCON//  , 5I6, L6 )
 1020 FORMAT(/ 38H //ADDCON//  GENERAL CONSTRAINT ADDED.
     *       / 49H //ADDCON//  NACTIV NCOLZ NFREE  IADD  JADD UNITQ
     *       / 13H //ADDCON//  , 5I6, L6 )
 2000 FORMAT(/ 12H *** WARNING
     *       / 48H *** SERIOUS ILL-CONDITIONING IN THE WORKING SET,
     *         25H AFTER ADDING CONSTRAINT ,  I5
     *       / 48H *** OVERFLOW MAY OCCUR IN SUBSEQUENT ITERATIONS //)
 3000 FORMAT(/ 42H //ADDCON//  DEPENDENT CONSTRAINT REJECTED )
 3010 FORMAT(/ 41H //ADDCON//     ASIZE     DTMAX     DTMIN
     *       / 11H //ADDCON//, 1P3E10.2 )
 3020 FORMAT(/ 51H //ADDCON//     ASIZE     DTMAX     DTMIN     DTNEW
     *       / 11H //ADDCON//, 1P4E10.2 )
C
C  END OF ADDCON
      END