SLATEC Routines --- SBOLSM ---


*DECK SBOLSM
      SUBROUTINE SBOLSM (W, MDW, MINPUT, NCOLS, BL, BU, IND, IOPT, X,
     +   RNORM, MODE, RW, WW, SCL, IBASIS, IBB)
C***BEGIN PROLOGUE  SBOLSM
C***SUBSIDIARY
C***PURPOSE  Subsidiary to SBOCLS and SBOLS
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (SBOLSM-S, DBOLSM-D)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C          Solve E*X = F (least squares sense) with bounds on
C            selected X values.
C     The user must have DIMENSION statements of the form:
C
C       DIMENSION W(MDW,NCOLS+1), BL(NCOLS), BU(NCOLS),
C      * X(NCOLS+NX), RW(NCOLS), WW(NCOLS), SCL(NCOLS)
C       INTEGER IND(NCOLS), IOPT(1+NI), IBASIS(NCOLS), IBB(NCOLS)
C
C     (Here NX=number of extra locations required for options 1,...,7;
C     NX=0 for no options; here NI=number of extra locations possibly
C     required for options 1-7; NI=0 for no options; NI=14 if all the
C     options are simultaneously in use.)
C
C    INPUT
C    -----
C
C    --------------------
C    W(MDW,*),MINPUT,NCOLS
C    --------------------
C     The array W(*,*) contains the matrix [E:F] on entry. The matrix
C     [E:F] has MINPUT rows and NCOLS+1 columns. This data is placed in
C     the array W(*,*) with E occupying the first NCOLS columns and the
C     right side vector F in column NCOLS+1. The row dimension, MDW, of
C     the array W(*,*) must satisfy the inequality MDW .ge. MINPUT.
C     Other values of MDW are errors. The values of MINPUT and NCOLS
C     must be positive. Other values are errors.
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.
C
C    1.    For IND(J)=1, require X(J) .ge. BL(J).
C    2.    For IND(J)=2, require X(J) .le. BU(J).
C    3.    For IND(J)=3, require X(J) .ge. BL(J) and
C                                X(J) .le. BU(J).
C    4.    For IND(J)=4, no bounds on X(J) are required.
C     The values of BL(*),BL(*) are modified by the subprogram. Values
C     other than 1,2,3 or 4 for IND(J) are errors. In the case IND(J)=3
C     (upper and lower bounds) the condition BL(J) .gt. BU(J) is an
C     error.
C
C    -------
C    IOPT(*)
C    -------
C     This is the array where the user can specify nonstandard options
C     for SBOLSM. 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         Move the IOPT(*) processing pointer.
C           2         Change rank determination tolerance.
C           3         Change blow-up factor that determines the
C                     size of variables being dropped from active
C                     status.
C           4         Reset the maximum number of iterations to use
C                     in solving the problem.
C           5         The data matrix is triangularized before the
C                     problem is solved whenever (NCOLS/MINPUT) .lt.
C                     FAC. Change the value of FAC.
C           6         Redefine the weighting matrix used for
C                     linear independence checking.
C           7         Debug output is desired.
C          99         No more options to change.
C
C    ----
C    X(*)
C    ----
C     This array is used to pass data associated with options 1,2,3 and
C     5. Ignore this input parameter if none of these options are used.
C     Otherwise see below: IOPT(*) CONTENTS.
C
C    ----------------
C    IBASIS(*),IBB(*)
C    ----------------
C     These arrays must be initialized by the user. The values
C         IBASIS(J)=J, J=1,...,NCOLS
C         IBB(J)   =1, J=1,...,NCOLS
C     are appropriate except when using nonstandard features.
C
C    ------
C    SCL(*)
C    ------
C     This is the array of scaling factors to use on the columns of the
C     matrix E. These values must be defined by the user. To suppress
C     any column scaling set SCL(J)=1.0, J=1,...,NCOLS.
C
C    OUTPUT
C    ------
C
C    ----------
C    X(*),RNORM
C    ----------
C     The array X(*) contains a solution (if MODE .ge. 0 or .eq. -22)
C     for the constrained least squares problem. The value RNORM is the
C     minimum residual vector length.
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.
C     A 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 18 cases -38,-37,...,-22, or -1. Values .lt. -1 correspond
C     to an abnormal completion of the subprogram. To understand the
C     abnormal completion codes see below: ERROR MESSAGES for SBOLSM
C     An approximate solution will be returned to the user only when
C     maximum iterations is reached, MODE=-22.
C
C    -----------
C    RW(*),WW(*)
C    -----------
C     These are working arrays each with NCOLS entries. The array RW(*)
C     contains the working (scaled, nonactive) solution values. The
C     array WW(*) contains the working (scaled, active) gradient vector
C     values.
C
C    ----------------
C    IBASIS(*),IBB(*)
C    ----------------
C     These arrays contain information about the status of the solution
C     when MODE .ge. 0. The indices IBASIS(K), K=1,...,MODE, show the
C     nonactive variables; indices IBASIS(K), K=MODE+1,..., NCOLS are
C     the active variables. The value (IBB(J)-1) is the number of times
C     variable J was reflected from its upper bound. (Normally the user
C     can ignore these parameters.)
C
C    IOPT(*) CONTENTS
C    ------- --------
C     The option array allows a user to modify internal variables in
C     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. The value is updated as the options are processed.  At the
C     pointer position the option number is extracted and used for
C     locating other information that allows for options to be changed.
C     The portion of the array IOPT(*) that is used for each option is
C     fixed; the user and the subprogram both know how many locations
C     are needed for each option. A great deal of error checking is
C     done by the subprogram on the contents of the option array.
C     Nevertheless it is still possible to give the subprogram optional
C     input that is meaningless. For example, some of the options use
C     the location X(NCOLS+IOFF) for passing data. The user must manage
C     the allocation of these locations when more than one piece of
C     option data is being passed to the subprogram.
C
C   1
C   -
C     Move the processing pointer (either forward or backward) to the
C     location IOPT(LP+1). The processing pointer is moved to location
C     LP+2 of IOPT(*) in case IOPT(LP)=-1.  For example to skip over
C     locations 3,...,NCOLS+2 of IOPT(*),
C
C       IOPT(1)=1
C       IOPT(2)=NCOLS+3
C       (IOPT(I), I=3,...,NCOLS+2 are not defined here.)
C       IOPT(NCOLS+3)=99
C       CALL SBOLSM
C
C     CAUTION: Misuse of this option can yield some very hard-to-find
C     bugs.  Use it with care.
C
C   2
C   -
C     The algorithm that solves the bounded least squares problem
C     iteratively drops columns from the active set. This has the
C     effect of joining a new column vector to the QR factorization of
C     the rectangular matrix consisting of the partially triangularized
C     nonactive columns. After triangularizing this matrix a test is
C     made on the size of the pivot element. The column vector is
C     rejected as dependent if the magnitude of the pivot element is
C     .le. TOL* magnitude of the column in components strictly above
C     the pivot element. Nominally the value of this (rank) tolerance
C     is TOL = SQRT(R1MACH(4)). To change only the value of TOL, for
C     example,
C
C       X(NCOLS+1)=TOL
C       IOPT(1)=2
C       IOPT(2)=1
C       IOPT(3)=99
C       CALL SBOLSM
C
C     Generally, if LP is the processing pointer for IOPT(*),
C
C       X(NCOLS+IOFF)=TOL
C       IOPT(LP)=2
C       IOPT(LP+1)=IOFF
C        .
C       CALL SBOLSM
C
C     The required length of IOPT(*) is increased by 2 if option 2 is
C     used; The required length of X(*) is increased by 1. A value of
C     IOFF .le. 0 is an error. A value of TOL .le. R1MACH(4) gives a
C     warning message; it is not considered an error.
C
C   3
C   -
C     A solution component is left active (not used) if, roughly
C     speaking, it seems too large. Mathematically the new component is
C     left active if the magnitude is .ge.((vector norm of F)/(matrix
C     norm of E))/BLOWUP. Nominally the factor BLOWUP = SQRT(R1MACH(4)).
C     To change only the value of BLOWUP, for example,
C
C       X(NCOLS+2)=BLOWUP
C       IOPT(1)=3
C       IOPT(2)=2
C       IOPT(3)=99
C       CALL SBOLSM
C
C     Generally, if LP is the processing pointer for IOPT(*),
C
C       X(NCOLS+IOFF)=BLOWUP
C       IOPT(LP)=3
C       IOPT(LP+1)=IOFF
C        .
C       CALL SBOLSM
C
C     The required length of IOPT(*) is increased by 2 if option 3 is
C     used; the required length of X(*) is increased by 1. A value of
C     IOFF .le. 0 is an error. A value of BLOWUP .le. 0.0 is an error.
C
C   4
C   -
C     Normally the algorithm for solving the bounded least squares
C     problem requires between NCOLS/3 and NCOLS drop-add steps to
C     converge. (this remark is based on examining a small number of
C     test cases.) The amount of arithmetic for such problems is
C     typically about twice that required for linear least squares if
C     there are no bounds and if plane rotations are used in the
C     solution method. Convergence of the algorithm, while
C     mathematically certain, can be much slower than indicated. To
C     avoid this potential but unlikely event ITMAX drop-add steps are
C     permitted. Nominally ITMAX=5*(MAX(MINPUT,NCOLS)). To change the
C     value of ITMAX, for example,
C
C       IOPT(1)=4
C       IOPT(2)=ITMAX
C       IOPT(3)=99
C       CALL SBOLSM
C
C     Generally, if LP is the processing pointer for IOPT(*),
C
C       IOPT(LP)=4
C       IOPT(LP+1)=ITMAX
C        .
C       CALL SBOLSM
C
C     The value of ITMAX must be .gt. 0. Other values are errors. Use
C     of this option increases the required length of IOPT(*) by 2.
C
C   5
C   -
C     For purposes of increased efficiency the MINPUT by NCOLS+1 data
C     matrix [E:F] is triangularized as a first step whenever MINPUT
C     satisfies FAC*MINPUT .gt. NCOLS. Nominally FAC=0.75. To change the
C     value of FAC,
C
C       X(NCOLS+3)=FAC
C       IOPT(1)=5
C       IOPT(2)=3
C       IOPT(3)=99
C       CALL SBOLSM
C
C     Generally, if LP is the processing pointer for IOPT(*),
C
C       X(NCOLS+IOFF)=FAC
C       IOPT(LP)=5
C       IOPT(LP+1)=IOFF
C        .
C       CALL SBOLSM
C
C     The value of FAC must be nonnegative. Other values are errors.
C     Resetting FAC=0.0 suppresses the initial triangularization step.
C     Use of this option increases the required length of IOPT(*) by 2;
C     The required length of of X(*) is increased by 1.
C
C   6
C   -
C     The norm used in testing the magnitudes of the pivot element
C     compared to the mass of the column above the pivot line can be
C     changed. The type of change that this option allows is to weight
C     the components with an index larger than MVAL by the parameter
C     WT. Normally MVAL=0 and WT=1. To change both the values MVAL and
C     WT, where LP is the processing pointer for IOPT(*),
C
C       X(NCOLS+IOFF)=WT
C       IOPT(LP)=6
C       IOPT(LP+1)=IOFF
C       IOPT(LP+2)=MVAL
C
C     Use of this option increases the required length of IOPT(*) by 3.
C     The length of X(*) is increased by 1. Values of MVAL must be
C     nonnegative and not greater than MINPUT. Other values are errors.
C     The value of WT must be positive. Any other value is an error. If
C     either error condition is present a message will be printed.
C
C   7
C   -
C     Debug output, showing the detailed add-drop steps for the
C     constrained least squares problem, is desired. This option is
C     intended to be used to locate suspected bugs.
C
C   99
C   --
C     There are no more options to change.
C
C     The values for options are 1,...,7,99, and are the only ones
C     permitted. Other values are errors. Options -99,-1,...,-7 mean
C     that the repective options 99,1,...,7 are left at their default
C     values. An example is the option to modify the (rank) tolerance:
C
C       X(NCOLS+1)=TOL
C       IOPT(1)=-2
C       IOPT(2)=1
C       IOPT(3)=99
C
C    Error Messages for SBOLSM
C    ----- -------- --- ---------
C    -22    MORE THAN ITMAX = ... ITERATIONS SOLVING BOUNDED LEAST
C           SQUARES PROBLEM.
C
C    -23    THE OPTION NUMBER = ... IS NOT DEFINED.
C
C    -24    THE OFFSET = ... BEYOND POSTION NCOLS = ... MUST BE POSITIVE
C           FOR OPTION NUMBER 2.
C
C    -25    THE TOLERANCE FOR RANK DETERMINATION = ... IS LESS THAN
C           MACHINE PRECISION = ....
C
C    -26    THE OFFSET = ... BEYOND POSITION NCOLS = ... MUST BE POSTIVE
C           FOR OPTION NUMBER 3.
C
C    -27    THE RECIPROCAL OF THE BLOW-UP FACTOR FOR REJECTING VARIABLES
C           MUST BE POSITIVE. NOW = ....
C
C    -28    THE MAXIMUM NUMBER OF ITERATIONS = ... MUST BE POSITIVE.
C
C    -29    THE OFFSET = ... BEYOND POSITION NCOLS = ... MUST BE POSTIVE
C           FOR OPTION NUMBER 5.
C
C    -30    THE FACTOR (NCOLS/MINPUT) WHERE PRETRIANGULARIZING IS
C           PERFORMED MUST BE NONNEGATIVE. NOW = ....
C
C    -31    THE NUMBER OF ROWS = ... MUST BE POSITIVE.
C
C    -32    THE NUMBER OF COLUMNS = ... MUST BE POSTIVE.
C
C    -33    THE ROW DIMENSION OF W(,) = ... MUST BE .GE. THE NUMBER OF
C           ROWS = ....
C
C    -34    FOR J = ... THE CONSTRAINT INDICATOR MUST BE 1-4.
C
C    -35    FOR J = ... THE LOWER BOUND = ... IS .GT. THE UPPER BOUND =
C           ....
C
C    -36    THE INPUT ORDER OF COLUMNS = ... IS NOT BETWEEN 1 AND NCOLS
C           = ....
C
C    -37    THE BOUND POLARITY FLAG IN COMPONENT J = ... MUST BE
C           POSITIVE. NOW = ....
C
C    -38    THE ROW SEPARATOR TO APPLY WEIGHTING (...) MUST LIE BETWEEN
C           0 AND MINPUT = .... WEIGHT = ... MUST BE POSITIVE.
C
C***SEE ALSO  SBOCLS, SBOLS
C***ROUTINES CALLED  IVOUT, R1MACH, SAXPY, SCOPY, SDOT, SMOUT, SNRM2,
C                    SROT, SROTG, SSWAP, SVOUT, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   821220  DATE WRITTEN
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
C   920422  Fixed usage of MINPUT.  (WRB)
C   901009  Editorial changes, code now reads from top to bottom.  (RWC)
C***END PROLOGUE  SBOLSM