*DECK DBOLSM SUBROUTINE DBOLSM (W, MDW, MINPUT, NCOLS, BL, BU, IND, IOPT, X, + RNORM, MODE, RW, WW, SCL, IBASIS, IBB) C***BEGIN PROLOGUE DBOLSM C***SUBSIDIARY C***PURPOSE Subsidiary to DBOCLS and DBOLS C***LIBRARY SLATEC C***TYPE DOUBLE PRECISION (SBOLSM-S, DBOLSM-D) C***AUTHOR (UNKNOWN) C***DESCRIPTION C C **** Double Precision Version of SBOLSM **** C **** All INPUT and OUTPUT real variables are DOUBLE PRECISION **** 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 DBOLSM. 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 DBOLSM 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 DBOLSM 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 DBOLSM 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 DBOLSM 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 DBOLSM 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 DBOLSM 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 DBOLSM C C Generally, if LP is the processing pointer for IOPT(*), C C IOPT(LP)=4 C IOPT(LP+1)=ITMAX C . C CALL DBOLSM 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 DBOLSM 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 DBOLSM 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 DBOLSM 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 DBOCLS, DBOLS C***ROUTINES CALLED D1MACH, DAXPY, DCOPY, DDOT, DMOUT, DNRM2, DROT, C DROTG, DSWAP, DVOUT, IVOUT, 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 DBOLSM