dlsei Subroutine

  Input.. All TYPE REAL variables are DOUBLE PRECISION

  W(*,*),MDW,   The array W(*,*) is doubly subscripted with
  ME,MA,MG,N    first dimensioning parameter equal to MDW.
                For this discussion let us call M = ME+MA+MG.  Then
                MDW must satisfy MDW >= M.  The condition
                MDW < M is an error.

                The array W(*,*) contains the matrices and vectors

                               (E  F)
                               (A  B)
                               (G  H)

                in rows and columns 1,...,M and 1,...,N+1
                respectively.

                The integers ME, MA, and MG are the
                respective matrix row dimensions
                of E, A and G.  Each matrix has N columns.

  PRGOPT(*)    This real-valued array is the option vector.
               If the user is satisfied with the nominal
               subprogram features set

               PRGOPT(1)=1 (or PRGOPT(1)=1.0)

               Otherwise PRGOPT(*) is a linked list consisting of
               groups of data of the following form

               LINK
               KEY
               DATA SET

               The parameters LINK and KEY are each one word.
               The DATA SET can be comprised of several words.
               The number of items depends on the value of KEY.
               The value of LINK points to the first
               entry of the next group of data within
               PRGOPT(*).  The exception is when there are
               no more options to change.  In that
               case, LINK=1 and the values KEY and DATA SET
               are not referenced.  The general layout of
               PRGOPT(*) is as follows.

            ...PRGOPT(1) = LINK1 (link to first entry of next group)
            .  PRGOPT(2) = KEY1 (key to the option change)
            .  PRGOPT(3) = data value (data value for this change)
            .       .
            .       .
            .       .
            ...PRGOPT(LINK1)   = LINK2 (link to the first entry of
            .                       next group)
            .  PRGOPT(LINK1+1) = KEY2 (key to the option change)
            .  PRGOPT(LINK1+2) = data value
            ...     .
            .       .
            .       .
            ...PRGOPT(LINK) = 1 (no more options to change)

               Values of LINK that are nonpositive are errors.
               A value of LINK > NLINK=100000 is also an error.
               This helps prevent using invalid but positive
               values of LINK that will probably extend
               beyond the program limits of PRGOPT(*).
               Unrecognized values of KEY are ignored.  The
               order of the options is arbitrary and any number
               of options can be changed with the following
               restriction.  To prevent cycling in the
               processing of the option array, a count of the
               number of options changed is maintained.
               Whenever this count exceeds NOPT=1000, an error
               message is printed and the subprogram returns.

               Options..

               KEY=1
                      Compute in W(*,*) the N by N
               covariance matrix of the solution variables
               as an output parameter.  Nominally the
               covariance matrix will not be computed.
               (This requires no user input.)
               The data set for this option is a single value.
               It must be nonzero when the covariance matrix
               is desired.  If it is zero, the covariance
               matrix is not computed.  When the covariance matrix
               is computed, the first dimensioning parameter
               of the array W(*,*) must satisfy MDW >= MAX(M,N).

               KEY=10
                      Suppress scaling of the inverse of the
               normal matrix by the scale factor RNORM**2/
               MAX(1, no. of degrees of freedom).  This option
               only applies when the option for computing the
               covariance matrix (KEY=1) is used.  With KEY=1 and
               KEY=10 used as options the unscaled inverse of the
               normal matrix is returned in W(*,*).
               The data set for this option is a single value.
               When it is nonzero no scaling is done.  When it is
               zero scaling is done.  The nominal case is to do
               scaling so if option (KEY=1) is used alone, the
               matrix will be scaled on output.

               KEY=2
                      Scale the nonzero columns of the
                      entire data matrix.
               (E)
               (A)
               (G)

               to have length one.  The data set for this
               option is a single value.  It must be
               nonzero if unit length column scaling
               is desired.

               KEY=3
                      Scale columns of the entire data matrix
               (E)
               (A)
               (G)

               with a user-provided diagonal matrix.
               The data set for this option consists
               of the N diagonal scaling factors, one for
               each matrix column.

               KEY=4
                      Change the rank determination tolerance for
               the equality constraint equations from
               the nominal value of SQRT(DRELPR).  This quantity can
               be no smaller than DRELPR, the arithmetic-
               storage precision.  The quantity DRELPR is the
               largest positive number such that T=1.+DRELPR
               satisfies T == 1.  The quantity used
               here is internally restricted to be at
               least DRELPR.  The data set for this option
               is the new tolerance.

               KEY=5
                      Change the rank determination tolerance for
               the reduced least squares equations from
               the nominal value of SQRT(DRELPR).  This quantity can
               be no smaller than DRELPR, the arithmetic-
               storage precision.  The quantity used
               here is internally restricted to be at
               least DRELPR.  The data set for this option
               is the new tolerance.

               For example, suppose we want to change
               the tolerance for the reduced least squares
               problem, compute the covariance matrix of
               the solution parameters, and provide
               column scaling for the data matrix.  For
               these options the dimension of PRGOPT(*)
               must be at least N+9.  The Fortran statements
               defining these options would be as follows:

               PRGOPT(1)=4 (link to entry 4 in PRGOPT(*))
               PRGOPT(2)=1 (covariance matrix key)
               PRGOPT(3)=1 (covariance matrix wanted)

               PRGOPT(4)=7 (link to entry 7 in PRGOPT(*))
               PRGOPT(5)=5 (least squares equas.  tolerance key)
               PRGOPT(6)=... (new value of the tolerance)

               PRGOPT(7)=N+9 (link to entry N+9 in PRGOPT(*))
               PRGOPT(8)=3 (user-provided column scaling key)

               CALL DCOPY (N, D, 1, PRGOPT(9), 1)  (Copy the N
                 scaling factors from the user array D(*)
                 to PRGOPT(9)-PRGOPT(N+8))

               PRGOPT(N+9)=1 (no more options to change)

               The contents of PRGOPT(*) are not modified
               by the subprogram.
               The options for WNNLS( ) can also be included
               in this array.  The values of KEY recognized
               by WNNLS( ) are 6, 7 and 8.  Their functions
               are documented in the usage instructions for
               subroutine WNNLS( ).  Normally these options
               do not need to be modified when using [[DLSEI]].

  IP(1),       The amounts of working storage actually
  IP(2)        allocated for the working arrays WS(*) and
               IP(*), respectively.  These quantities are
               compared with the actual amounts of storage
               needed by [[DLSEI]].  Insufficient storage
               allocated for either WS(*) or IP(*) is an
               error.  This feature was included in [[DLSEI]]
               because miscalculating the storage formulas
               for WS(*) and IP(*) might very well lead to
               subtle and hard-to-find execution errors.

               The length of WS(*) must be at least

               LW = 2*(ME+N)+K+(MG+2)*(N+7)

               where K = max(MA+MG,N)
               This test will not be made if IP(1)<=0.

               The length of IP(*) must be at least

               LIP = MG+2*N+2
               This test will not be made if IP(2)<=0.

  Output.. All TYPE REAL variables are DOUBLE PRECISION

  X(*),RNORME,  The array X(*) contains the solution parameters
  RNORML        if the integer output flag MODE = 0 or 1.
                The definition of MODE is given directly below.
                When MODE = 0 or 1, RNORME and RNORML
                respectively contain the residual vector
                Euclidean lengths of F - EX and B - AX.  When
                MODE=1 the equality constraint equations EX=F
                are contradictory, so RNORME /= 0.  The residual
                vector F-EX has minimal Euclidean length.  For
                MODE >= 2, none of these parameters is defined.

  MODE          Integer flag that indicates the subprogram
                status after completion.  If MODE >= 2, no
                solution has been computed.

                MODE =

                0  Both equality and inequality constraints
                   are compatible and have been satisfied.

                1  Equality constraints are contradictory.
                   A generalized inverse solution of EX=F was used
                   to minimize the residual vector length F-EX.
                   In this sense, the solution is still meaningful.

                2  Inequality constraints are contradictory.

                3  Both equality and inequality constraints
                   are contradictory.

                The following interpretation of
                MODE=1,2 or 3 must be made.  The
                sets consisting of all solutions
                of the equality constraints EX=F
                and all vectors satisfying GX >= H
                have no points in common.  (In
                particular this does not say that
                each individual set has no points
                at all, although this could be the
                case.)

                4  Usage error occurred.  The value
                   of MDW is < ME+MA+MG, MDW is
                   < N and a covariance matrix is
                   requested, or the option vector
                   PRGOPT(*) is not properly defined,
                   or the lengths of the working arrays
                   WS(*) and IP(*), when specified in
                   IP(1) and IP(2) respectively, are not
                   long enough.

  W(*,*)        The array W(*,*) contains the N by N symmetric
                covariance matrix of the solution parameters,
                provided this was requested on input with
                the option vector PRGOPT(*) and the output
                flag is returned with MODE = 0 or 1.

  IP(*)         The integer working array has three entries
                that provide rank and working array length
                information after completion.

                   IP(1) = rank of equality constraint
                           matrix.  Define this quantity
                           as KRANKE.

                   IP(2) = rank of reduced least squares
                           problem.

                   IP(3) = the amount of storage in the
                           working array WS(*) that was
                           actually used by the subprogram.
                           The formula given above for the length
                           of WS(*) is a necessary overestimate.
                           If exactly the same problem matrices
                           are used in subsequent executions,
                           the declared dimension of WS(*) can
                           be reduced to this output value.
  User Designated
  Working Arrays..

  WS(*),IP(*)              These are respectively type real
                           and type integer working arrays.
                           Their required minimal lengths are
                           given above.