--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/npsol/lsadd.f	Fri Jul 19 01:29:55 1996 +0000
@@ -0,0 +1,301 @@
+*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+*     File  LSSUBS FORTRAN
+*
+*     LSADD    LSADDS   LSBNDS   LSCHOL   LSCORE   LSCRSH   LSDEL
+*     LSDFLT   LSFEAS   LSFILE   LSGETP   LSGSET   LSKEY    LSLOC
+*     LSMOVE   LSMULS   LSOPTN   LSPRT    LSSETX   LSSOL
+*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+      SUBROUTINE LSADD ( UNITQ,
+     $                   INFORM, IFIX, IADD, JADD,
+     $                   NACTIV, NZ, NFREE, NRANK, NRES, NGQ,
+     $                   N, NROWA, NQ, NROWR, NROWT,
+     $                   KX, CONDMX,
+     $                   A, R, T, RES, GQ, ZY, WRK1, WRK2 )
+
+      IMPLICIT           DOUBLE PRECISION(A-H,O-Z)
+      LOGICAL            UNITQ
+      INTEGER            KX(N)
+      DOUBLE PRECISION   A(NROWA,*), R(NROWR,*), T(NROWT,*),
+     $                   RES(N,*), GQ(N,*), ZY(NQ,*)
+      DOUBLE PRECISION   WRK1(N), WRK2(N)
+************************************************************************
+*  LSADD   updates the factorization,  A(free) * (Z Y) = (0 T),  when a
+*  constraint is added to the working set.  If  NRANK .gt. 0, the
+*  factorization  ( R ) = PWQ  is also updated,  where  W  is the
+*                 ( 0 )
+*  least squares matrix,  R  is upper-triangular,  and  P  is an
+*  orthogonal matrix.  The matrices  W  and  P  are not stored.
+*
+*  There are three separate cases to consider (although each case
+*  shares code with another)...
+*
+*  (1) A free variable becomes fixed on one of its bounds when there
+*      are already some general constraints in the working set.
+*
+*  (2) A free variable becomes fixed on one of its bounds when there
+*      are only bound constraints in the working set.
+*
+*  (3) A general constraint (corresponding to row  IADD  of  A) is
+*      added to the working set.
+*
+*  In cases (1) and (2), we assume that  KX(IFIX) = JADD.
+*  In all cases,  JADD  is the index of the constraint being added.
+*
+*  If there are no general constraints in the working set,  the
+*  matrix  Q = (Z Y)  is the identity and will not be touched.
+*
+*  If  NRES .GT. 0,  the row transformations are applied to the rows of
+*  the  (N by NRES)  matrix  RES.
+*  If  NGQ .GT. 0,  the column transformations are applied to the
+*  columns of the  (NGQ by N)  matrix  GQ'.
+*
+*  Systems Optimization Laboratory, Stanford University.
+*  Original version written 31-October--1984.
+*  This version of LSADD dated 29-December-1985.
+************************************************************************
+      COMMON    /SOL1CM/ NOUT
+      COMMON    /SOL4CM/ EPSPT3, EPSPT5, EPSPT8, EPSPT9
+      COMMON    /SOL5CM/ ASIZE, DTMAX, DTMIN
+
+      LOGICAL            LSDBG
+      PARAMETER         (LDBG = 5)
+      COMMON    /LSDEBG/ ILSDBG(LDBG), LSDBG
+
+      LOGICAL            BOUND , OVERFL
+      EXTERNAL           DDOT  , DDIV  , DNRM2
+      INTRINSIC          MAX   , MIN
+      PARAMETER         (ZERO = 0.0D+0, ONE = 1.0D+0)
+
+*     If the condition estimator of the updated factors is greater than
+*     CONDBD,  a warning message is printed.
+
+      CONDBD = ONE / EPSPT9
+
+      OVERFL = .FALSE.
+      BOUND  = JADD .LE. N
+      IF (BOUND) THEN
+*        ===============================================================
+*        A simple bound has entered the working set.  IADD  is not used.
+*        ===============================================================
+         IF (LSDBG  .AND.  ILSDBG(1) .GT. 0)
+     $      WRITE (NOUT, 1010) NACTIV, NZ, NFREE, IFIX, JADD, UNITQ
+         NANEW = NACTIV
+
+         IF (UNITQ) THEN
+
+*           Q  is not stored, but KX defines an ordering of the columns
+*           of the identity matrix that implicitly define  Q.
+*           Reorder KX so that variable IFIX is moved to position
+*           NFREE+1 and variables IFIX+1,...,NFREE+1 are moved one
+*           position to the left.
+
+            CALL DLOAD ( NFREE, (ZERO), WRK1, 1 )
+            WRK1(IFIX) = ONE
+
+            DO 100 I = IFIX, NFREE-1
+               KX(I) = KX(I+1)
+  100       CONTINUE
+         ELSE
+*           ------------------------------------------------------------
+*           Q  is stored explicitly.
+*           ------------------------------------------------------------
+*           Set  WRK1 = the  (IFIX)-th  row of  Q.
+*           Move the  (NFREE)-th  row of  Q  to position  IFIX.
+
+            CALL DCOPY ( NFREE, ZY(IFIX,1), NQ, WRK1, 1 )
+            IF (IFIX .LT. NFREE) THEN
+               CALL DCOPY ( NFREE, ZY(NFREE,1), NQ, ZY(IFIX,1), NQ )
+               KX(IFIX) = KX(NFREE)
+            END IF
+         END IF
+         KX(NFREE) = JADD
+      ELSE
+*        ===============================================================
+*        A general constraint has entered the working set.
+*        IFIX  is not used.
+*        ===============================================================
+         IF (LSDBG  .AND.  ILSDBG(1) .GT. 0)
+     $      WRITE (NOUT, 1020) NACTIV, NZ, NFREE, IADD, JADD, UNITQ
+
+         NANEW  = NACTIV + 1
+
+*        Transform the incoming row of  A  by  Q'.
+
+         CALL DCOPY ( N, A(IADD,1), NROWA, WRK1, 1 )
+         CALL CMQMUL( 8, N, NZ, NFREE, NQ, UNITQ, KX, WRK1, ZY, WRK2)
+
+*        Check that the incoming row is not dependent upon those
+*        already in the working set.
+
+         DTNEW  = DNRM2 ( NZ, WRK1, 1 )
+         IF (NACTIV .EQ. 0) THEN
+
+*           This is the only general constraint in the working set.
+
+            COND   = DDIV  ( ASIZE, DTNEW, OVERFL )
+            TDTMAX = DTNEW
+            TDTMIN = DTNEW
+         ELSE
+
+*           There are already some general constraints in the working
+*           set. Update the estimate of the condition number.
+
+            TDTMAX = MAX( DTNEW, DTMAX )
+            TDTMIN = MIN( DTNEW, DTMIN )
+            COND   = DDIV  ( TDTMAX, TDTMIN, OVERFL )
+         END IF
+
+         IF (COND .GT. CONDMX  .OR.  OVERFL) GO TO 900
+
+         IF (UNITQ) THEN
+
+*           This is the first general constraint to be added.
+*           Set  Q = I.
+
+            DO 200 J = 1, NFREE
+               CALL DLOAD ( NFREE, (ZERO), ZY(1,J), 1 )
+               ZY(J,J) = ONE
+  200       CONTINUE
+            UNITQ  = .FALSE.
+         END IF
+      END IF
+
+      NZERO  = NZ - 1
+      IF (BOUND) NZERO = NFREE - 1
+
+*     ------------------------------------------------------------------
+*     Use a sequence of 2*2 column transformations to reduce the
+*     first NZERO elements of WRK1 to zero.  This affects ZY, except
+*     when UNITQ is true.  The transformations may also be applied
+*     to R, T and GQ'.
+*     ------------------------------------------------------------------
+      LROWR  = N
+      NELM   = 1
+      IROWT  = NACTIV
+
+      DO 300 K = 1, NZERO
+
+*        Compute the transformation that reduces WRK1(K) to zero,
+*        then apply it to the relevant columns of  Z  and  GQ'.
+
+         CALL DROT3G( WRK1(K+1), WRK1(K), CS, SN )
+         IF (.NOT. UNITQ)
+     $      CALL DROT3 ( NFREE, ZY(1,K+1), 1, ZY(1,K), 1, CS, SN )
+         IF (NGQ .GT. 0)
+     $      CALL DROT3 ( NGQ  , GQ(K+1,1), N, GQ(K,1), N, CS, SN )
+
+         IF (K .GE. NZ  .AND.  NACTIV .GT. 0) THEN
+
+*           Apply the rotation to  T.
+
+            T(IROWT,K) = ZERO
+            CALL DROT3 ( NELM, T(IROWT,K+1), 1, T(IROWT,K), 1, CS, SN )
+            NELM  = NELM  + 1
+            IROWT = IROWT - 1
+         END IF
+
+         IF (NRANK .GT. 0) THEN
+
+*           Apply the same transformation to the columns of R.
+*           This generates a subdiagonal element in R that must be
+*           eliminated by a row rotation.
+
+            IF (K .LT. NRANK) R(K+1,K) = ZERO
+            LCOL   = MIN( K+1, NRANK )
+
+            CALL DROT3 ( LCOL, R(1,K+1), 1, R(1,K), 1, CS, SN )
+            IF (K .LT. NRANK) THEN
+               CALL DROT3G( R(K,K), R(K+1,K), CS, SN )
+               LROWR  = LROWR - 1
+               CALL DROT3 ( LROWR,   R(K,K+1)  , NROWR,
+     $                               R(K+1,K+1), NROWR, CS, SN )
+
+               IF (NRES .GT. 0)
+     $            CALL DROT3 ( NRES, RES(K,1)  , N    ,
+     $                               RES(K+1,1), N    , CS, SN )
+            END IF
+         END IF
+  300 CONTINUE
+
+      IF (BOUND) THEN
+
+*        The last row and column of ZY has been transformed to plus
+*        or minus the unit vector E(NFREE).  We can reconstitute the
+*        columns of GQ and R corresponding to the new fixed variable.
+
+         IF (WRK1(NFREE) .LT. ZERO) THEN
+            NFMIN = MIN( NRANK, NFREE )
+            IF (NFMIN .GT. 0) CALL DSCAL ( NFMIN, -ONE, R(1,NFREE) , 1 )
+            IF (NGQ   .GT. 0) CALL DSCAL ( NGQ  , -ONE, GQ(NFREE,1), N )
+         END IF
+
+*        ---------------------------------------------------------------
+*        The diagonals of T have been altered.  Recompute the
+*        largest and smallest values.
+*        ---------------------------------------------------------------
+         IF (NACTIV .GT. 0) THEN
+            CALL DCOND( NACTIV, T(NACTIV,NZ), NROWT-1, TDTMAX, TDTMIN )
+            COND   = DDIV  ( TDTMAX, TDTMIN, OVERFL )
+         END IF
+      ELSE
+*        ---------------------------------------------------------------
+*        General constraint.  Install the new row of T.
+*        ---------------------------------------------------------------
+         CALL DCOPY ( NANEW, WRK1(NZ), 1, T(NANEW,NZ), NROWT )
+      END IF
+
+*     ==================================================================
+*     Prepare to exit.  Check the magnitude of the condition estimator.
+*     ==================================================================
+  900 IF (NANEW .GT. 0) THEN
+         IF (COND .LT. CONDMX  .AND.  .NOT. OVERFL) THEN
+
+*           The factorization has been successfully updated.
+
+            INFORM = 0
+            DTMAX  = TDTMAX
+            DTMIN  = TDTMIN
+            IF (COND .GE. CONDBD) WRITE (NOUT, 2000) JADD
+         ELSE
+
+*           The proposed working set appears to be linearly dependent.
+
+            INFORM = 1
+            IF (LSDBG  .AND.  ILSDBG(1) .GT. 0) THEN
+               WRITE( NOUT, 3000 )
+               IF (BOUND) THEN
+                  WRITE (NOUT, 3010) ASIZE, DTMAX, DTMIN
+               ELSE
+                  IF (NACTIV .GT. 0) THEN
+                     WRITE (NOUT, 3020) ASIZE, DTMAX, DTMIN, DTNEW
+                  ELSE
+                     WRITE (NOUT, 3030) ASIZE, DTNEW
+                  END IF
+               END IF
+            END IF
+         END IF
+      END IF
+
+      RETURN
+
+ 1010 FORMAT(/ ' //LSADD //  Simple bound added.'
+     $       / ' //LSADD //  NACTIV    NZ NFREE  IFIX  JADD UNITQ'
+     $       / ' //LSADD //  ', 5I6, L6 )
+ 1020 FORMAT(/ ' //LSADD //  General constraint added.           '
+     $       / ' //LSADD //  NACTIV    NZ NFREE  IADD  JADD UNITQ'
+     $       / ' //LSADD //  ', 5I6, L6 )
+ 2000 FORMAT(/ ' XXX  Serious ill-conditioning in the working set',
+     $         ' after adding constraint ',  I5
+     $       / ' XXX  Overflow may occur in subsequent iterations.'//)
+ 3000 FORMAT(/ ' //LSADD //  Dependent constraint rejected.' )
+ 3010 FORMAT(/ ' //LSADD //     ASIZE     DTMAX     DTMIN        '
+     $       / ' //LSADD //', 1P3E10.2 )
+ 3020 FORMAT(/ ' //LSADD //     ASIZE     DTMAX     DTMIN     DTNEW'
+     $       / ' //LSADD //', 1P4E10.2 )
+ 3030 FORMAT(/ ' //LSADD //     ASIZE     DTNEW'
+     $       / ' //LSADD //', 1P2E10.2 )
+
+*     End of  LSADD .
+
+      END