*DECK DLSEI 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