SLATEC Routines --- LSEI ---
*DECK LSEI
SUBROUTINE LSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML,
+ MODE, WS, IP)
C***BEGIN PROLOGUE LSEI
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 SINGLE 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 LSEI( ) are
C
C Input..
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(SRELPR). This quantity can
C be no smaller than SRELPR, the arithmetic-
C storage precision. The quantity SRELPR is the
C largest positive number such that T=1.+SRELPR
C satisfies T .EQ. 1. The quantity used
C here is internally restricted to be at
C least SRELPR. 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(SRELPR). This quantity can
C be no smaller than SRELPR, the arithmetic-
C storage precision. The quantity used
C here is internally restricted to be at
C least SRELPR. 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 SCOPY (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 LSEI( ).
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 LSEI( ). Insufficient storage
C allocated for either WS(*) or IP(*) is an
C error. This feature was included in LSEI( )
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..
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 H12, LSI, R1MACH, SASUM, SAXPY, SCOPY, SDOT, SNRM2,
C SSCAL, SSWAP, 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 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE LSEI