# dbocls.f

          SUBROUTINE DBOCLS (W, MDW, MCON, MROWS, NCOLS, BL, BU, IND, IOPT,
         +   X, RNORMC, RNORM, MODE, RW, IW)
    C***BEGIN PROLOGUE  DBOCLS
    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      DOUBLE PRECISION (SBOCLS-S, DBOCLS-D)
    C***KEYWORDS  BOUNDS, CONSTRAINTS, INEQUALITY, LEAST SQUARES, LINEAR
    C***AUTHOR  Hanson, R. J., (SNLA)
    C***DESCRIPTION
    C
    C   **** All INPUT and OUTPUT real variables are DOUBLE PRECISION ****
    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     See IOPT(*) description.
    C
    C
    C    -------
    C    IOPT(*)
    C    -------
    C     This is the array where the user can specify nonstandard options
    C     for DBOCLS( ). 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 DBOCLS(),
    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 DBOCLS( ) 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      .