SBOLSM(3)

SYNOPSIS
SUBROUTINE SBOLSM (W, MDW, MINPUT, NCOLS, BL, BU, IND, IOPT, X, + RNORM, MODE, RW, WW, SCL, IBASIS, IBB)

NAME
sbolsm - Subsidiary to SBOCLS and SBOLS

LIBRARY
SLATEC

TYPE
SINGLE PRECISION (SBOLSM-S, DBOLSM-D)

AUTHOR
(UNKNOWN)

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

SEE ALSO
SBOCLS, SBOLS

ROUTINES CALLED
IVOUT, R1MACH, SAXPY, SCOPY, SDOT, SMOUT, SNRM2, SROT, SROTG, SSWAP, SVOUT, XERMSG

REVISION HISTORY (YYMMDD)
821220 DATE WRITTEN
891214 Prologue converted to Version 4.0 format. (BAB)
900328 Added TYPE section. (WRB)
900510 Convert XERRWV calls to XERMSG calls. (RWC)
920422 Fixed usage of MINPUT. (WRB)
901009 Editorial changes, code now reads from top to bottom. (RWC)