SLATEC Routines --- SBOCLS ---
*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 .R1MACH, SASUM, SBOLS, SCOPY, SDOT, SNRM2, SSCAL,
C XERMSG
C***REVISION HISTORY (YYMMDD)
C 821220 DATE WRITTEN
C 870803 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900510 Convert XERRWV calls to XERMSG calls. (RWC)
C 910819 Added variable M for MOUT+MCON in reference to SBOLS. (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE SBOCLS