# dlsei.f

          SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME,
         +   RNORML, MODE, WS, IP)
    C***BEGIN PROLOGUE  DLSEI
    C***PURPOSE  Solve a linearly constrained least squares problem with
    C            equality and inequality constraints, and optionally compute
    C            a covariance matrix.
    C***LIBRARY   SLATEC
    C***CATEGORY  K1A2A, D9
    C***TYPE      DOUBLE PRECISION (LSEI-S, DLSEI-D)
    C***KEYWORDS  CONSTRAINED LEAST SQUARES, CURVE FITTING, DATA FITTING,
    C             EQUALITY CONSTRAINTS, INEQUALITY CONSTRAINTS,
    C             QUADRATIC PROGRAMMING
    C***AUTHOR  Hanson, R. J., (SNLA)
    C           Haskell, K. H., (SNLA)
    C***DESCRIPTION
    C
    C     Abstract
    C
    C     This subprogram solves a linearly constrained least squares
    C     problem with both equality and inequality constraints, and, if the
    C     user requests, obtains a covariance matrix of the solution
    C     parameters.
    C
    C     Suppose there are given matrices E, A and G of respective
    C     dimensions ME by N, MA by N and MG by N, and vectors F, B and H of
    C     respective lengths ME, MA and MG.  This subroutine solves the
    C     linearly constrained least squares problem
    C
    C                   EX = F, (E ME by N) (equations to be exactly
    C                                       satisfied)
    C                   AX = B, (A MA by N) (equations to be
    C                                       approximately satisfied,
    C                                       least squares sense)
    C                   GX .GE. H,(G MG by N) (inequality constraints)
    C
    C     The inequalities GX .GE. H mean that every component of the
    C     product GX must be .GE. the corresponding component of H.
    C
    C     In case the equality constraints cannot be satisfied, a
    C     generalized inverse solution residual vector length is obtained
    C     for F-EX.  This is the minimal length possible for F-EX.
    C
    C     Any values ME .GE. 0, MA .GE. 0, or MG .GE. 0 are permitted.  The
    C     rank of the matrix E is estimated during the computation.  We call
    C     this value KRANKE.  It is an output parameter in IP(1) defined
    C     below.  Using a generalized inverse solution of EX=F, a reduced
    C     least squares problem with inequality constraints is obtained.
    C     The tolerances used in these tests for determining the rank
    C     of E and the rank of the reduced least squares problem are
    C     given in Sandia Tech. Rept. SAND-78-1290.  They can be
    C     modified by the user if new values are provided in
    C     the option list of the array PRGOPT(*).
    C
    C     The user must dimension all arrays appearing in the call list..
    C     W(MDW,N+1),PRGOPT(*),X(N),WS(2*(ME+N)+K+(MG+2)*(N+7)),IP(MG+2*N+2)
    C     where K=MAX(MA+MG,N).  This allows for a solution of a range of
    C     problems in the given working space.  The dimension of WS(*)
    C     given is a necessary overestimate.  Once a particular problem
    C     has been run, the output parameter IP(3) gives the actual
    C     dimension required for that problem.
    C
    C     The parameters for DLSEI( ) are
    C
    C     Input.. All TYPE REAL variables are DOUBLE PRECISION
    C
    C     W(*,*),MDW,   The array W(*,*) is doubly subscripted with
    C     ME,MA,MG,N    first dimensioning parameter equal to MDW.
    C                   For this discussion let us call M = ME+MA+MG.  Then
    C                   MDW must satisfy MDW .GE. M.  The condition
    C                   MDW .LT. M is an error.
    C
    C                   The array W(*,*) contains the matrices and vectors
    C
    C                                  (E  F)
    C                                  (A  B)
    C                                  (G  H)
    C
    C                   in rows and columns 1,...,M and 1,...,N+1
    C                   respectively.
    C
    C                   The integers ME, MA, and MG are the
    C                   respective matrix row dimensions
    C                   of E, A and G.  Each matrix has N columns.
    C
    C     PRGOPT(*)    This real-valued array is the option vector.
    C                  If the user is satisfied with the nominal
    C                  subprogram features set
    C
    C                  PRGOPT(1)=1 (or PRGOPT(1)=1.0)
    C
    C                  Otherwise PRGOPT(*) is a linked list consisting of
    C                  groups of data of the following form
    C
    C                  LINK
    C                  KEY
    C                  DATA SET
    C
    C                  The parameters LINK and KEY are each one word.
    C                  The DATA SET can be comprised of several words.
    C                  The number of items depends on the value of KEY.
    C                  The value of LINK points to the first
    C                  entry of the next group of data within
    C                  PRGOPT(*).  The exception is when there are
    C                  no more options to change.  In that
    C                  case, LINK=1 and the values KEY and DATA SET
    C                  are not referenced.  The general layout of
    C                  PRGOPT(*) is as follows.
    C
    C               ...PRGOPT(1) = LINK1 (link to first entry of next group)
    C               .  PRGOPT(2) = KEY1 (key to the option change)
    C               .  PRGOPT(3) = data value (data value for this change)
    C               .       .
    C               .       .
    C               .       .
    C               ...PRGOPT(LINK1)   = LINK2 (link to the first entry of
    C               .                       next group)
    C               .  PRGOPT(LINK1+1) = KEY2 (key to the option change)
    C               .  PRGOPT(LINK1+2) = data value
    C               ...     .
    C               .       .
    C               .       .
    C               ...PRGOPT(LINK) = 1 (no more options to change)
    C
    C                  Values of LINK that are nonpositive are errors.
    C                  A value of LINK .GT. NLINK=100000 is also an error.
    C                  This helps prevent using invalid but positive
    C                  values of LINK that will probably extend
    C                  beyond the program limits of PRGOPT(*).
    C                  Unrecognized values of KEY are ignored.  The
    C                  order of the options is arbitrary and any number
    C                  of options can be changed with the following
    C                  restriction.  To prevent cycling in the
    C                  processing of the option array, a count of the
    C                  number of options changed is maintained.
    C                  Whenever this count exceeds NOPT=1000, an error
    C                  message is printed and the subprogram returns.
    C
    C                  Options..
    C
    C                  KEY=1
    C                         Compute in W(*,*) the N by N
    C                  covariance matrix of the solution variables
    C                  as an output parameter.  Nominally the
    C                  covariance matrix will not be computed.
    C                  (This requires no user input.)
    C                  The data set for this option is a single value.
    C                  It must be nonzero when the covariance matrix
    C                  is desired.  If it is zero, the covariance
    C                  matrix is not computed.  When the covariance matrix
    C                  is computed, the first dimensioning parameter
    C                  of the array W(*,*) must satisfy MDW .GE. MAX(M,N).
    C
    C                  KEY=10
    C                         Suppress scaling of the inverse of the
    C                  normal matrix by the scale factor RNORM**2/
    C                  MAX(1, no. of degrees of freedom).  This option
    C                  only applies when the option for computing the
    C                  covariance matrix (KEY=1) is used.  With KEY=1 and
    C                  KEY=10 used as options the unscaled inverse of the
    C                  normal matrix is returned in W(*,*).
    C                  The data set for this option is a single value.
    C                  When it is nonzero no scaling is done.  When it is
    C                  zero scaling is done.  The nominal case is to do
    C                  scaling so if option (KEY=1) is used alone, the
    C                  matrix will be scaled on output.
    C
    C                  KEY=2
    C                         Scale the nonzero columns of the
    C                         entire data matrix.
    C                  (E)
    C                  (A)
    C                  (G)
    C
    C                  to have length one.  The data set for this
    C                  option is a single value.  It must be
    C                  nonzero if unit length column scaling
    C                  is desired.
    C
    C                  KEY=3
    C                         Scale columns of the entire data matrix
    C                  (E)
    C                  (A)
    C                  (G)
    C
    C                  with a user-provided diagonal matrix.
    C                  The data set for this option consists
    C                  of the N diagonal scaling factors, one for
    C                  each matrix column.
    C
    C                  KEY=4
    C                         Change the rank determination tolerance for
    C                  the equality constraint equations from
    C                  the nominal value of SQRT(DRELPR).  This quantity can
    C                  be no smaller than DRELPR, the arithmetic-
    C                  storage precision.  The quantity DRELPR is the
    C                  largest positive number such that T=1.+DRELPR
    C                  satisfies T .EQ. 1.  The quantity used
    C                  here is internally restricted to be at
    C                  least DRELPR.  The data set for this option
    C                  is the new tolerance.
    C
    C                  KEY=5
    C                         Change the rank determination tolerance for
    C                  the reduced least squares equations from
    C                  the nominal value of SQRT(DRELPR).  This quantity can
    C                  be no smaller than DRELPR, the arithmetic-
    C                  storage precision.  The quantity used
    C                  here is internally restricted to be at
    C                  least DRELPR.  The data set for this option
    C                  is the new tolerance.
    C
    C                  For example, suppose we want to change
    C                  the tolerance for the reduced least squares
    C                  problem, compute the covariance matrix of
    C                  the solution parameters, and provide
    C                  column scaling for the data matrix.  For
    C                  these options the dimension of PRGOPT(*)
    C                  must be at least N+9.  The Fortran statements
    C                  defining these options would be as follows:
    C
    C                  PRGOPT(1)=4 (link to entry 4 in PRGOPT(*))
    C                  PRGOPT(2)=1 (covariance matrix key)
    C                  PRGOPT(3)=1 (covariance matrix wanted)
    C
    C                  PRGOPT(4)=7 (link to entry 7 in PRGOPT(*))
    C                  PRGOPT(5)=5 (least squares equas.  tolerance key)
    C                  PRGOPT(6)=... (new value of the tolerance)
    C
    C                  PRGOPT(7)=N+9 (link to entry N+9 in PRGOPT(*))
    C                  PRGOPT(8)=3 (user-provided column scaling key)
    C
    C                  CALL DCOPY (N, D, 1, PRGOPT(9), 1)  (Copy the N
    C                    scaling factors from the user array D(*)
    C                    to PRGOPT(9)-PRGOPT(N+8))
    C
    C                  PRGOPT(N+9)=1 (no more options to change)
    C
    C                  The contents of PRGOPT(*) are not modified
    C                  by the subprogram.
    C                  The options for WNNLS( ) can also be included
    C                  in this array.  The values of KEY recognized
    C                  by WNNLS( ) are 6, 7 and 8.  Their functions
    C                  are documented in the usage instructions for
    C                  subroutine WNNLS( ).  Normally these options
    C                  do not need to be modified when using DLSEI( ).
    C
    C     IP(1),       The amounts of working storage actually
    C     IP(2)        allocated for the working arrays WS(*) and
    C                  IP(*), respectively.  These quantities are
    C                  compared with the actual amounts of storage
    C                  needed by DLSEI( ).  Insufficient storage
    C                  allocated for either WS(*) or IP(*) is an
    C                  error.  This feature was included in DLSEI( )
    C                  because miscalculating the storage formulas
    C                  for WS(*) and IP(*) might very well lead to
    C                  subtle and hard-to-find execution errors.
    C
    C                  The length of WS(*) must be at least
    C
    C                  LW = 2*(ME+N)+K+(MG+2)*(N+7)
    C
    C                  where K = max(MA+MG,N)
    C                  This test will not be made if IP(1).LE.0.
    C
    C                  The length of IP(*) must be at least
    C
    C                  LIP = MG+2*N+2
    C                  This test will not be made if IP(2).LE.0.
    C
    C     Output.. All TYPE REAL variables are DOUBLE PRECISION
    C
    C     X(*),RNORME,  The array X(*) contains the solution parameters
    C     RNORML        if the integer output flag MODE = 0 or 1.
    C                   The definition of MODE is given directly below.
    C                   When MODE = 0 or 1, RNORME and RNORML
    C                   respectively contain the residual vector
    C                   Euclidean lengths of F - EX and B - AX.  When
    C                   MODE=1 the equality constraint equations EX=F
    C                   are contradictory, so RNORME .NE. 0.  The residual
    C                   vector F-EX has minimal Euclidean length.  For
    C                   MODE .GE. 2, none of these parameters is defined.
    C
    C     MODE          Integer flag that indicates the subprogram
    C                   status after completion.  If MODE .GE. 2, no
    C                   solution has been computed.
    C
    C                   MODE =
    C
    C                   0  Both equality and inequality constraints
    C                      are compatible and have been satisfied.
    C
    C                   1  Equality constraints are contradictory.
    C                      A generalized inverse solution of EX=F was used
    C                      to minimize the residual vector length F-EX.
    C                      In this sense, the solution is still meaningful.
    C
    C                   2  Inequality constraints are contradictory.
    C
    C                   3  Both equality and inequality constraints
    C                      are contradictory.
    C
    C                   The following interpretation of
    C                   MODE=1,2 or 3 must be made.  The
    C                   sets consisting of all solutions
    C                   of the equality constraints EX=F
    C                   and all vectors satisfying GX .GE. H
    C                   have no points in common.  (In
    C                   particular this does not say that
    C                   each individual set has no points
    C                   at all, although this could be the
    C                   case.)
    C
    C                   4  Usage error occurred.  The value
    C                      of MDW is .LT. ME+MA+MG, MDW is
    C                      .LT. N and a covariance matrix is
    C                      requested, or the option vector
    C                      PRGOPT(*) is not properly defined,
    C                      or the lengths of the working arrays
    C                      WS(*) and IP(*), when specified in
    C                      IP(1) and IP(2) respectively, are not
    C                      long enough.
    C
    C     W(*,*)        The array W(*,*) contains the N by N symmetric
    C                   covariance matrix of the solution parameters,
    C                   provided this was requested on input with
    C                   the option vector PRGOPT(*) and the output
    C                   flag is returned with MODE = 0 or 1.
    C
    C     IP(*)         The integer working array has three entries
    C                   that provide rank and working array length
    C                   information after completion.
    C
    C                      IP(1) = rank of equality constraint
    C                              matrix.  Define this quantity
    C                              as KRANKE.
    C
    C                      IP(2) = rank of reduced least squares
    C                              problem.
    C
    C                      IP(3) = the amount of storage in the
    C                              working array WS(*) that was
    C                              actually used by the subprogram.
    C                              The formula given above for the length
    C                              of WS(*) is a necessary overestimate.
    C                              If exactly the same problem matrices
    C                              are used in subsequent executions,
    C                              the declared dimension of WS(*) can
    C                              be reduced to this output value.
    C     User Designated
    C     Working Arrays..
    C
    C     WS(*),IP(*)              These are respectively type real
    C                              and type integer working arrays.
    C                              Their required minimal lengths are
    C                              given above.
    C
    C***REFERENCES  K. H. Haskell and R. J. Hanson, An algorithm for
    C                 linear least squares problems with equality and
    C                 nonnegativity constraints, Report SAND77-0552, Sandia
    C                 Laboratories, June 1978.
    C               K. H. Haskell and R. J. Hanson, Selected algorithms for
    C                 the linearly constrained least squares problem - a
    C                 users guide, Report SAND78-1290, Sandia Laboratories,
    C                 August 1979.
    C               K. H. Haskell and R. J. Hanson, An algorithm for
    C                 linear least squares problems with equality and
    C                 nonnegativity constraints, Mathematical Programming
    C                 21 (1981), pp. 98-118.
    C               R. J. Hanson and K. H. Haskell, Two algorithms for the
    C                 linearly constrained least squares problem, ACM
    C                 Transactions on Mathematical Software, September 1982.
    C***ROUTINES CALLED  D1MACH, DASUM, DAXPY, DCOPY, DDOT, DH12, DLSI,
    C                    DNRM2, DSCAL, DSWAP, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   790701  DATE WRITTEN
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    C   890618  Completely restructured and extensively revised (WRB & RWC)
    C   890831  REVISION DATE from Version 3.2
    C   891214  Prologue converted to Version 4.0 format.  (BAB)
    C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
    C   900604  DP version created from SP version.  (RWC)
    C   920501  Reformatted the REFERENCES section.  (WRB)
    C***END PROLOGUE  DLSEI