*DECK SBOCLS SUBROUTINE SBOCLS (W, MDW, MCON, MROWS, NCOLS, BL, BU, IND, IOPT, + X, RNORMC, RNORM, MODE, RW, IW) C***BEGIN PROLOGUE SBOCLS C***PURPOSE Solve the bounded and constrained least squares C problem consisting of solving the equation C E*X = F (in the least squares sense) C subject to the linear constraints C C*X = Y. C***LIBRARY SLATEC C***CATEGORY K1A2A, G2E, G2H1, G2H2 C***TYPE SINGLE PRECISION (SBOCLS-S, DBOCLS-D) C***KEYWORDS BOUNDS, CONSTRAINTS, INEQUALITY, LEAST SQUARES, LINEAR C***AUTHOR Hanson, R. J., (SNLA) C***DESCRIPTION C C This subprogram solves the bounded and constrained least squares C problem. The problem statement is: C C Solve E*X = F (least squares sense), subject to constraints C C*X=Y. C C In this formulation both X and Y are unknowns, and both may C have bounds on any of their components. This formulation C of the problem allows the user to have equality and inequality C constraints as well as simple bounds on the solution components. C C This constrained linear least squares subprogram solves E*X=F C subject to C*X=Y, where E is MROWS by NCOLS, C is MCON by NCOLS. C C The user must have dimension statements of the form C C DIMENSION W(MDW,NCOLS+MCON+1), BL(NCOLS+MCON), BU(NCOLS+MCON), C * X(2*(NCOLS+MCON)+2+NX), RW(6*NCOLS+5*MCON) C INTEGER IND(NCOLS+MCON), IOPT(17+NI), IW(2*(NCOLS+MCON)) C C (here NX=number of extra locations required for the options; NX=0 C if no options are in use. Also NI=number of extra locations C for options 1-9.) C C INPUT C ----- C C ------------------------- C W(MDW,*),MCON,MROWS,NCOLS C ------------------------- C The array W contains the (possibly null) matrix [C:*] followed by C [E:F]. This must be placed in W as follows: C [C : *] C W = [ ] C [E : F] C The (*) after C indicates that this data can be undefined. The C matrix [E:F] has MROWS rows and NCOLS+1 columns. The matrix C is C placed in the first MCON rows of W(*,*) while [E:F] C follows in rows MCON+1 through MCON+MROWS of W(*,*). The vector F C is placed in rows MCON+1 through MCON+MROWS, column NCOLS+1. The C values of MDW and NCOLS must be positive; the value of MCON must C be nonnegative. An exception to this occurs when using option 1 C for accumulation of blocks of equations. In that case MROWS is an C OUTPUT variable only, and the matrix data for [E:F] is placed in C W(*,*), one block of rows at a time. See IOPT(*) contents, option C number 1, for further details. The row dimension, MDW, of the C array W(*,*) must satisfy the inequality: C C If using option 1, C MDW .ge. MCON + max(max. number of C rows accumulated, NCOLS) + 1. C If using option 8, C MDW .ge. MCON + MROWS. C Else C MDW .ge. MCON + max(MROWS, NCOLS). C C Other values are errors, but this is checked only when using C option=2. The value of MROWS is an output parameter when C using option number 1 for accumulating large blocks of least C squares equations before solving the problem. C See IOPT(*) contents for details about option 1. C C ------------------ C BL(*),BU(*),IND(*) C ------------------ C These arrays contain the information about the bounds that the C solution values are to satisfy. The value of IND(J) tells the C type of bound and BL(J) and BU(J) give the explicit values for C the respective upper and lower bounds on the unknowns X and Y. C The first NVARS entries of IND(*), BL(*) and BU(*) specify C bounds on X; the next MCON entries specify bounds on Y. C C 1. For IND(J)=1, require X(J) .ge. BL(J); C IF J.gt.NCOLS, Y(J-NCOLS) .ge. BL(J). C (the value of BU(J) is not used.) C 2. For IND(J)=2, require X(J) .le. BU(J); C IF J.gt.NCOLS, Y(J-NCOLS) .le. BU(J). C (the value of BL(J) is not used.) C 3. For IND(J)=3, require X(J) .ge. BL(J) and C X(J) .le. BU(J); C IF J.gt.NCOLS, Y(J-NCOLS) .ge. BL(J) and C Y(J-NCOLS) .le. BU(J). C (to impose equality constraints have BL(J)=BU(J)= C constraining value.) C 4. For IND(J)=4, no bounds on X(J) or Y(J-NCOLS) are required. C (the values of BL(J) and BU(J) are not used.) C C Values other than 1,2,3 or 4 for IND(J) are errors. In the case C IND(J)=3 (upper and lower bounds) the condition BL(J) .gt. BU(J) C is an error. The values BL(J), BU(J), J .gt. NCOLS, will be C changed. Significant changes mean that the constraints are C infeasible. (Users must make this decision themselves.) C The new values for BL(J), BU(J), J .gt. NCOLS, define a C region such that the perturbed problem is feasible. If users C know that their problem is feasible, this step can be skipped C by using option number 8 described below. C C See IOPT(*) description. C C C ------- C IOPT(*) C ------- C This is the array where the user can specify nonstandard options C for SBOCLS( ). Most of the time this feature can be ignored by C setting the input value IOPT(1)=99. Occasionally users may have C needs that require use of the following subprogram options. For C details about how to use the options see below: IOPT(*) CONTENTS. C C Option Number Brief Statement of Purpose C ------ ------ ----- --------- -- ------- C 1 Return to user for accumulation of blocks C of least squares equations. The values C of IOPT(*) are changed with this option. C The changes are updates to pointers for C placing the rows of equations into position C for processing. C 2 Check lengths of all arrays used in the C subprogram. C 3 Column scaling of the data matrix, [C]. C [E] C 4 User provides column scaling for matrix [C]. C [E] C 5 Provide option array to the low-level C subprogram SBOLS( ). C 6 Provide option array to the low-level C subprogram SBOLSM( ). C 7 Move the IOPT(*) processing pointer. C 8 Do not preprocess the constraints to C resolve infeasibilities. C 9 Do not pretriangularize the least squares matrix. C 99 No more options to change. C C ---- C X(*) C ---- C This array is used to pass data associated with options 4,5 and C 6. Ignore this parameter (on input) if no options are used. C Otherwise see below: IOPT(*) CONTENTS. C C C OUTPUT C ------ C C ----------------- C X(*),RNORMC,RNORM C ----------------- C The array X(*) contains a solution (if MODE .ge.0 or .eq.-22) for C the constrained least squares problem. The value RNORMC is the C minimum residual vector length for the constraints C*X - Y = 0. C The value RNORM is the minimum residual vector length for the C least squares equations. Normally RNORMC=0, but in the case of C inconsistent constraints this value will be nonzero. C The values of X are returned in the first NVARS entries of X(*). C The values of Y are returned in the last MCON entries of X(*). C C ---- C MODE C ---- C The sign of MODE determines whether the subprogram has completed C normally, or encountered an error condition or abnormal status. A C value of MODE .ge. 0 signifies that the subprogram has completed C normally. The value of mode (.ge. 0) is the number of variables C in an active status: not at a bound nor at the value zero, for C the case of free variables. A negative value of MODE will be one C of the cases (-57)-(-41), (-37)-(-22), (-19)-(-2). Values .lt. -1 C correspond to an abnormal completion of the subprogram. These C error messages are in groups for the subprograms SBOCLS(), C SBOLSM(), and SBOLS(). An approximate solution will be returned C to the user only when max. iterations is reached, MODE=-22. C C ----------- C RW(*),IW(*) C ----------- C These are working arrays. (normally the user can ignore the C contents of these arrays.) C C IOPT(*) CONTENTS C ------- -------- C The option array allows a user to modify some internal variables C in the subprogram without recompiling the source code. A central C goal of the initial software design was to do a good job for most C people. Thus the use of options will be restricted to a select C group of users. The processing of the option array proceeds as C follows: a pointer, here called LP, is initially set to the value C 1. At the pointer position the option number is extracted and C used for locating other information that allows for options to be C changed. The portion of the array IOPT(*) that is used for each C option is fixed; the user and the subprogram both know how many C locations are needed for each option. The value of LP is updated C for each option based on the amount of storage in IOPT(*) that is C required. A great deal of error checking is done by the C subprogram on the contents of the option array. Nevertheless it C is still possible to give the subprogram optional input that is C meaningless. For example option 4 uses the locations C X(NCOLS+IOFF),...,X(NCOLS+IOFF+NCOLS-1) for passing scaling data. C The user must manage the allocation of these locations. C C 1 C - C This option allows the user to solve problems with a large number C of rows compared to the number of variables. The idea is that the C subprogram returns to the user (perhaps many times) and receives C new least squares equations from the calling program unit. C Eventually the user signals "that's all" and a solution is then C computed. The value of MROWS is an output variable when this C option is used. Its value is always in the range 0 .le. MROWS C .le. NCOLS+1. It is the number of rows after the C triangularization of the entire set of equations. If LP is the C processing pointer for IOPT(*), the usage for the sequential C processing of blocks of equations is C C C IOPT(LP)=1 C Move block of equations to W(*,*) starting at C the first row of W(*,*). C IOPT(LP+3)=# of rows in the block; user defined C C The user now calls SBOCLS( ) in a loop. The value of IOPT(LP+1) C directs the user's action. The value of IOPT(LP+2) points to C where the subsequent rows are to be placed in W(*,*). Both of C these values are first defined in the subprogram. The user C changes the value of IOPT(LP+1) (to 2) as a signal that all of C the rows have been processed. C C C .