*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 C C PURPOSE C ------- C THIS IS THE MAIN SUBPROGRAM THAT SOLVES THE BOUNDED C LEAST SQUARES PROBLEM. THE PROBLEM SOLVED HERE IS: C C SOLVE E*X = F (LEAST SQUARES SENSE) C WITH BOUNDS ON SELECTED X VALUES. C C TO CHANGE THIS SUBPROGRAM FROM SINGLE TO DOUBLE PRECISION BEGIN C EDITING AT THE CARD 'C++'. C CHANGE THE SUBPROGRAM NAME TO DBOLSM AND THE STRINGS C /SAXPY/ TO /DAXPY/, /SCOPY/ TO /DCOPY/, C /SDOT/ TO /DDOT/, /SNRM2/ TO /DNRM2/, C /SROT/ TO /DROT/, /SROTG/ TO /DROTG/, /R1MACH/ TO /D1MACH/, C /SVOUT/ TO /DVOUT/, /SMOUT/ TO /DMOUT/, C /SSWAP/ TO /DSWAP/, /E0/ TO /D0/, C /REAL / TO /DOUBLE PRECISION/. C++ C DOUBLE PRECISION W(MDW,*),BL(*),BU(*) DOUBLE PRECISION X(*),RW(*),WW(*),SCL(*) DOUBLE PRECISION ALPHA,BETA,BOU,COLABV,COLBLO DOUBLE PRECISION CL1,CL2,CL3,ONE,BIG DOUBLE PRECISION FAC,RNORM,SC,SS,T,TOLIND,WT DOUBLE PRECISION TWO,T1,T2,WBIG,WLARGE,WMAG,XNEW DOUBLE PRECISION ZERO,DDOT,DNRM2 DOUBLE PRECISION D1MACH,TOLSZE INTEGER IBASIS(*),IBB(*),IND(*),IOPT(*) LOGICAL FOUND,CONSTR CHARACTER*8 XERN1, XERN2 CHARACTER*16 XERN3, XERN4 C PARAMETER (ZERO=0.0D0, ONE=1.0D0, TWO=2.0D0) C INEXT(IDUM) = MIN(IDUM+1,MROWS) C***FIRST EXECUTABLE STATEMENT DBOLSM C C Verify that the problem dimensions are defined properly. C IF (MINPUT.LE.0) THEN WRITE (XERN1, '(I8)') MINPUT CALL XERMSG ('SLATEC', 'DBOLSM', 'THE NUMBER OF ROWS = ' // * XERN1 // ' MUST BE POSITIVE.', 31, 1) MODE = -31 RETURN ENDIF C IF (NCOLS.LE.0) THEN WRITE (XERN1, '(I8)') NCOLS CALL XERMSG ('SLATEC', 'DBOLSM', 'THE NUMBER OF COLUMNS = ' // * XERN1 // ' MUST BE POSITIVE.', 32, 1) MODE = -32 RETURN ENDIF C IF (MDW.LT.MINPUT) THEN WRITE (XERN1, '(I8)') MDW WRITE (XERN2, '(I8)') MINPUT CALL XERMSG ('SLATEC', 'DBOLSM', * 'THE ROW DIMENSION OF W(,) = ' // XERN1 // * ' MUST BE .GE. THE NUMBER OF ROWS = ' // XERN2, 33, 1) MODE = -33 RETURN ENDIF C C Verify that bound information is correct. C DO 10 J = 1,NCOLS IF (IND(J).LT.1 .OR. IND(J).GT.4) THEN WRITE (XERN1, '(I8)') J WRITE (XERN2, '(I8)') IND(J) CALL XERMSG ('SLATEC', 'DBOLSM', 'FOR J = ' // XERN1 // * ' THE CONSTRAINT INDICATOR MUST BE 1-4', 34, 1) MODE = -34 RETURN ENDIF 10 CONTINUE C DO 20 J = 1,NCOLS IF (IND(J).EQ.3) THEN IF (BU(J).LT.BL(J)) THEN WRITE (XERN1, '(I8)') J WRITE (XERN3, '(1PD15.6)') BL(J) WRITE (XERN4, '(1PD15.6)') BU(J) CALL XERMSG ('SLATEC', 'DBOLSM', 'FOR J = ' // XERN1 * // ' THE LOWER BOUND = ' // XERN3 // * ' IS .GT. THE UPPER BOUND = ' // XERN4, 35, 1) MODE = -35 RETURN ENDIF ENDIF 20 CONTINUE C C Check that permutation and polarity arrays have been set. C DO 30 J = 1,NCOLS IF (IBASIS(J).LT.1 .OR. IBASIS(J).GT.NCOLS) THEN WRITE (XERN1, '(I8)') IBASIS(J) WRITE (XERN2, '(I8)') NCOLS CALL XERMSG ('SLATEC', 'DBOLSM', * 'THE INPUT ORDER OF COLUMNS = ' // XERN1 // * ' IS NOT BETWEEN 1 AND NCOLS = ' // XERN2, 36, 1) MODE = -36 RETURN ENDIF C IF (IBB(J).LE.0) THEN WRITE (XERN1, '(I8)') J WRITE (XERN2, '(I8)') IBB(J) CALL XERMSG ('SLATEC', 'DBOLSM', * 'THE BOUND POLARITY FLAG IN COMPONENT J = ' // XERN1 // * ' MUST BE POSITIVE.$$NOW = ' // XERN2, 37, 1) MODE = -37 RETURN ENDIF 30 CONTINUE C C Process the option array. C FAC = 0.75D0 TOLIND = SQRT(D1MACH(4)) TOLSZE = SQRT(D1MACH(4)) ITMAX = 5*MAX(MINPUT,NCOLS) WT = ONE MVAL = 0 IPRINT = 0 C C Changes to some parameters can occur through the option array, C IOPT(*). Process this array looking carefully for input data C errors. C LP = 0 LDS = 0 C C Test for no more options. C 590 LP = LP + LDS IP = IOPT(LP+1) JP = ABS(IP) IF (IP.EQ.99) THEN GO TO 470 ELSE IF (JP.EQ.99) THEN LDS = 1 ELSE IF (JP.EQ.1) THEN C C Move the IOPT(*) processing pointer. C IF (IP.GT.0) THEN LP = IOPT(LP+2) - 1 LDS = 0 ELSE LDS = 2 ENDIF ELSE IF (JP.EQ.2) THEN C C Change tolerance for rank determination. C IF (IP.GT.0) THEN IOFF = IOPT(LP+2) IF (IOFF.LE.0) THEN WRITE (XERN1, '(I8)') IOFF WRITE (XERN2, '(I8)') NCOLS CALL XERMSG ('SLATEC', 'DBOLSM', 'THE OFFSET = ' // * XERN1 // ' BEYOND POSITION NCOLS = ' // XERN2 // * ' MUST BE POSITIVE FOR OPTION NUMBER 2.', 24, 1) MODE = -24 RETURN ENDIF C TOLIND = X(NCOLS+IOFF) IF (TOLIND.LT.D1MACH(4)) THEN WRITE (XERN3, '(1PD15.6)') TOLIND WRITE (XERN4, '(1PD15.6)') D1MACH(4) CALL XERMSG ('SLATEC', 'DBOLSM', * 'THE TOLERANCE FOR RANK DETERMINATION = ' // XERN3 * // ' IS LESS THAN MACHINE PRECISION = ' // XERN4, * 25, 0) MODE = -25 ENDIF ENDIF LDS = 2 ELSE IF (JP.EQ.3) THEN C C Change blowup factor for allowing variables to become C inactive. C IF (IP.GT.0) THEN IOFF = IOPT(LP+2) IF (IOFF.LE.0) THEN WRITE (XERN1, '(I8)') IOFF WRITE (XERN2, '(I8)') NCOLS CALL XERMSG ('SLATEC', 'DBOLSM', 'THE OFFSET = ' // * XERN1 // ' BEYOND POSITION NCOLS = ' // XERN2 // * ' MUST BE POSITIVE FOR OPTION NUMBER 3.', 26, 1) MODE = -26 RETURN ENDIF C TOLSZE = X(NCOLS+IOFF) IF (TOLSZE.LE.ZERO) THEN WRITE (XERN3, '(1PD15.6)') TOLSZE CALL XERMSG ('SLATEC', 'DBOLSM', 'THE RECIPROCAL ' // * 'OF THE BLOW-UP FACTOR FOR REJECTING VARIABLES ' // * 'MUST BE POSITIVE.$$NOW = ' // XERN3, 27, 1) MODE = -27 RETURN ENDIF ENDIF LDS = 2 ELSE IF (JP.EQ.4) THEN C C Change the maximum number of iterations allowed. C IF (IP.GT.0) THEN ITMAX = IOPT(LP+2) IF (ITMAX.LE.0) THEN WRITE (XERN1, '(I8)') ITMAX CALL XERMSG ('SLATEC', 'DBOLSM', * 'THE MAXIMUM NUMBER OF ITERATIONS = ' // XERN1 // * ' MUST BE POSITIVE.', 28, 1) MODE = -28 RETURN ENDIF ENDIF LDS = 2 ELSE IF (JP.EQ.5) THEN C C Change the factor for pretriangularizing the data matrix. C IF (IP.GT.0) THEN IOFF = IOPT(LP+2) IF (IOFF.LE.0) THEN WRITE (XERN1, '(I8)') IOFF WRITE (XERN2, '(I8)') NCOLS CALL XERMSG ('SLATEC', 'DBOLSM', 'THE OFFSET = ' // * XERN1 // ' BEYOND POSITION NCOLS = ' // XERN2 // * ' MUST BE POSITIVE FOR OPTION NUMBER 5.', 29, 1) MODE = -29 RETURN ENDIF C FAC = X(NCOLS+IOFF) IF (FAC.LT.ZERO) THEN WRITE (XERN3, '(1PD15.6)') FAC CALL XERMSG ('SLATEC', 'DBOLSM', * 'THE FACTOR (NCOLS/MINPUT) WHERE PRE-' // * 'TRIANGULARIZING IS PERFORMED MUST BE NON-' // * 'NEGATIVE.$$NOW = ' // XERN3, 30, 0) MODE = -30 RETURN ENDIF ENDIF LDS = 2 ELSE IF (JP.EQ.6) THEN C C Change the weighting factor (from 1.0) to apply to components C numbered .gt. MVAL (initially set to 1.) This trick is needed C for applications of this subprogram to the heavily weighted C least squares problem that come from equality constraints. C IF (IP.GT.0) THEN IOFF = IOPT(LP+2) MVAL = IOPT(LP+3) WT = X(NCOLS+IOFF) ENDIF C IF (MVAL.LT.0 .OR. MVAL.GT.MINPUT .OR. WT.LE.ZERO) THEN WRITE (XERN1, '(I8)') MVAL WRITE (XERN2, '(I8)') MINPUT WRITE (XERN3, '(1PD15.6)') WT CALL XERMSG ('SLATEC', 'DBOLSM', * 'THE ROW SEPARATOR TO APPLY WEIGHTING (' // XERN1 // * ') MUST LIE BETWEEN 0 AND MINPUT = ' // XERN2 // * '.$$WEIGHT = ' // XERN3 // ' MUST BE POSITIVE.', 38, 0) MODE = -38 RETURN ENDIF LDS = 3 ELSE IF (JP.EQ.7) THEN C C Turn on debug output. C IF (IP.GT.0) IPRINT = 1 LDS = 2 ELSE WRITE (XERN1, '(I8)') IP CALL XERMSG ('SLATEC', 'DBOLSM', 'THE OPTION NUMBER = ' // * XERN1 // ' IS NOT DEFINED.', 23, 1) MODE = -23 RETURN ENDIF GO TO 590 C C Pretriangularize rectangular arrays of certain sizes for C increased efficiency. C 470 IF (FAC*MINPUT.GT.NCOLS) THEN DO 490 J = 1,NCOLS+1 DO 480 I = MINPUT,J+MVAL+1,-1 CALL DROTG(W(I-1,J),W(I,J),SC,SS) W(I,J) = ZERO CALL DROT(NCOLS-J+1,W(I-1,J+1),MDW,W(I,J+1),MDW,SC,SS) 480 CONTINUE 490 CONTINUE MROWS = NCOLS + MVAL + 1 ELSE MROWS = MINPUT ENDIF C C Set the X(*) array to zero so all components are defined. C CALL DCOPY(NCOLS,ZERO,0,X,1) C C The arrays IBASIS(*) and IBB(*) are initialized by the calling C program and the column scaling is defined in the calling program. C 'BIG' is plus infinity on this machine. C BIG = D1MACH(2) DO 550 J = 1,NCOLS IF (IND(J).EQ.1) THEN BU(J) = BIG ELSE IF (IND(J).EQ.2) THEN BL(J) = -BIG ELSE IF (IND(J).EQ.4) THEN BL(J) = -BIG BU(J) = BIG ENDIF 550 CONTINUE C DO 570 J = 1,NCOLS IF ((BL(J).LE.ZERO.AND.ZERO.LE.BU(J).AND.ABS(BU(J)).LT. * ABS(BL(J))) .OR. BU(J).LT.ZERO) THEN T = BU(J) BU(J) = -BL(J) BL(J) = -T SCL(J) = -SCL(J) DO 560 I = 1,MROWS W(I,J) = -W(I,J) 560 CONTINUE ENDIF C C Indices in set T(=TIGHT) are denoted by negative values C of IBASIS(*). C IF (BL(J).GE.ZERO) THEN IBASIS(J) = -IBASIS(J) T = -BL(J) BU(J) = BU(J) + T CALL DAXPY(MROWS,T,W(1,J),1,W(1,NCOLS+1),1) ENDIF 570 CONTINUE C NSETB = 0 ITER = 0 C IF (IPRINT.GT.0) THEN CALL DMOUT(MROWS,NCOLS+1,MDW,W,'('' PRETRI. INPUT MATRIX'')', * -4) CALL DVOUT(NCOLS,BL,'('' LOWER BOUNDS'')',-4) CALL DVOUT(NCOLS,BU,'('' UPPER BOUNDS'')',-4) ENDIF C 580 ITER = ITER + 1 IF (ITER.GT.ITMAX) THEN WRITE (XERN1, '(I8)') ITMAX CALL XERMSG ('SLATEC', 'DBOLSM', 'MORE THAN ITMAX = ' // XERN1 * // ' ITERATIONS SOLVING BOUNDED LEAST SQUARES PROBLEM.', * 22, 1) MODE = -22 C C Rescale and translate variables. C IGOPR = 1 GO TO 130 ENDIF C C Find a variable to become non-active. C T C Compute (negative) of gradient vector, W = E *(F-E*X). C CALL DCOPY(NCOLS,ZERO,0,WW,1) DO 200 J = NSETB+1,NCOLS JCOL = ABS(IBASIS(J)) WW(J) = DDOT(MROWS-NSETB,W(INEXT(NSETB),J),1, * W(INEXT(NSETB),NCOLS+1),1)*ABS(SCL(JCOL)) 200 CONTINUE C IF (IPRINT.GT.0) THEN CALL DVOUT(NCOLS,WW,'('' GRADIENT VALUES'')',-4) CALL IVOUT(NCOLS,IBASIS,'('' INTERNAL VARIABLE ORDER'')',-4) CALL IVOUT(NCOLS,IBB,'('' BOUND POLARITY'')',-4) ENDIF C C If active set = number of total rows, quit. C 210 IF (NSETB.EQ.MROWS) THEN FOUND = .FALSE. GO TO 120 ENDIF C C Choose an extremal component of gradient vector for a candidate C to become non-active. C WLARGE = -BIG WMAG = -BIG DO 220 J = NSETB+1,NCOLS T = WW(J) IF (T.EQ.BIG) GO TO 220 ITEMP = IBASIS(J) JCOL = ABS(ITEMP) T1 = DNRM2(MVAL-NSETB,W(INEXT(NSETB),J),1) IF (ITEMP.LT.0) THEN IF (MOD(IBB(JCOL),2).EQ.0) T = -T IF (T.LT.ZERO) GO TO 220 IF (MVAL.GT.NSETB) T = T1 IF (T.GT.WLARGE) THEN WLARGE = T JLARGE = J ENDIF ELSE IF (MVAL.GT.NSETB) T = T1 IF (ABS(T).GT.WMAG) THEN WMAG = ABS(T) JMAG = J ENDIF ENDIF 220 CONTINUE C C Choose magnitude of largest component of gradient for candidate. C JBIG = 0 WBIG = ZERO IF (WLARGE.GT.ZERO) THEN JBIG = JLARGE WBIG = WLARGE ENDIF C IF (WMAG.GE.WBIG) THEN JBIG = JMAG WBIG = WMAG ENDIF C IF (JBIG.EQ.0) THEN FOUND = .FALSE. IF (IPRINT.GT.0) THEN CALL IVOUT(0,I,'('' FOUND NO VARIABLE TO ENTER'')',-4) ENDIF GO TO 120 ENDIF C C See if the incoming column is sufficiently independent. This C test is made before an elimination is performed. C IF (IPRINT.GT.0) * CALL IVOUT(1,JBIG,'('' TRY TO BRING IN THIS COL.'')',-4) C IF (MVAL.LE.NSETB) THEN CL1 = DNRM2(MVAL,W(1,JBIG),1) CL2 = ABS(WT)*DNRM2(NSETB-MVAL,W(INEXT(MVAL),JBIG),1) CL3 = ABS(WT)*DNRM2(MROWS-NSETB,W(INEXT(NSETB),JBIG),1) CALL DROTG(CL1,CL2,SC,SS) COLABV = ABS(CL1) COLBLO = CL3 ELSE CL1 = DNRM2(NSETB,W(1,JBIG),1) CL2 = DNRM2(MVAL-NSETB,W(INEXT(NSETB),JBIG),1) CL3 = ABS(WT)*DNRM2(MROWS-MVAL,W(INEXT(MVAL),JBIG),1) COLABV = CL1 CALL DROTG(CL2,CL3,SC,SS) COLBLO = ABS(CL2) ENDIF C IF (COLBLO.LE.TOLIND*COLABV) THEN WW(JBIG) = BIG IF (IPRINT.GT.0) * CALL IVOUT(0,I,'('' VARIABLE IS DEPENDENT, NOT USED.'')', * -4) GO TO 210 ENDIF C C Swap matrix columns NSETB+1 and JBIG, plus pointer information, C and gradient values. C NSETB = NSETB + 1 IF (NSETB.NE.JBIG) THEN CALL DSWAP(MROWS,W(1,NSETB),1,W(1,JBIG),1) CALL DSWAP(1,WW(NSETB),1,WW(JBIG),1) ITEMP = IBASIS(NSETB) IBASIS(NSETB) = IBASIS(JBIG) IBASIS(JBIG) = ITEMP ENDIF C C Eliminate entries below the pivot line in column NSETB. C IF (MROWS.GT.NSETB) THEN DO 230 I = MROWS,NSETB+1,-1 IF (I.EQ.MVAL+1) GO TO 230 CALL DROTG(W(I-1,NSETB),W(I,NSETB),SC,SS) W(I,NSETB) = ZERO CALL DROT(NCOLS-NSETB+1,W(I-1,NSETB+1),MDW,W(I,NSETB+1), * MDW,SC,SS) 230 CONTINUE C IF (MVAL.GE.NSETB .AND. MVAL.LT.MROWS) THEN CALL DROTG(W(NSETB,NSETB),W(MVAL+1,NSETB),SC,SS) W(MVAL+1,NSETB) = ZERO CALL DROT(NCOLS-NSETB+1,W(NSETB,NSETB+1),MDW, * W(MVAL+1,NSETB+1),MDW,SC,SS) ENDIF ENDIF C IF (W(NSETB,NSETB).EQ.ZERO) THEN WW(NSETB) = BIG NSETB = NSETB - 1 IF (IPRINT.GT.0) THEN CALL IVOUT(0,I,'('' PIVOT IS ZERO, NOT USED.'')',-4) ENDIF GO TO 210 ENDIF C C Check that new variable is moving in the right direction. C ITEMP = IBASIS(NSETB) JCOL = ABS(ITEMP) XNEW = (W(NSETB,NCOLS+1)/W(NSETB,NSETB))/ABS(SCL(JCOL)) IF (ITEMP.LT.0) THEN C C IF(WW(NSETB).GE.ZERO.AND.XNEW.LE.ZERO) exit(quit) C IF(WW(NSETB).LE.ZERO.AND.XNEW.GE.ZERO) exit(quit) C IF ((WW(NSETB).GE.ZERO.AND.XNEW.LE.ZERO) .OR. * (WW(NSETB).LE.ZERO.AND.XNEW.GE.ZERO)) GO TO 240 ENDIF FOUND = .TRUE. GO TO 120 C 240 WW(NSETB) = BIG NSETB = NSETB - 1 IF (IPRINT.GT.0) * CALL IVOUT(0,I,'('' VARIABLE HAS BAD DIRECTION, NOT USED.'')', * -4) GO TO 210 C C Solve the triangular system. C 270 CALL DCOPY(NSETB,W(1,NCOLS+1),1,RW,1) DO 280 J = NSETB,1,-1 RW(J) = RW(J)/W(J,J) JCOL = ABS(IBASIS(J)) T = RW(J) IF (MOD(IBB(JCOL),2).EQ.0) RW(J) = -RW(J) CALL DAXPY(J-1,-T,W(1,J),1,RW,1) RW(J) = RW(J)/ABS(SCL(JCOL)) 280 CONTINUE C IF (IPRINT.GT.0) THEN CALL DVOUT(NSETB,RW,'('' SOLN. VALUES'')',-4) CALL IVOUT(NSETB,IBASIS,'('' COLS. USED'')',-4) ENDIF C IF (LGOPR.EQ.2) THEN CALL DCOPY(NSETB,RW,1,X,1) DO 450 J = 1,NSETB ITEMP = IBASIS(J) JCOL = ABS(ITEMP) IF (ITEMP.LT.0) THEN BOU = ZERO ELSE BOU = BL(JCOL) ENDIF C IF ((-BOU).NE.BIG) BOU = BOU/ABS(SCL(JCOL)) IF (X(J).LE.BOU) THEN JDROP1 = J GO TO 340 ENDIF C BOU = BU(JCOL) IF (BOU.NE.BIG) BOU = BOU/ABS(SCL(JCOL)) IF (X(J).GE.BOU) THEN JDROP2 = J GO TO 340 ENDIF 450 CONTINUE GO TO 340 ENDIF C C See if the unconstrained solution (obtained by solving the C triangular system) satisfies the problem bounds. C ALPHA = TWO BETA = TWO X(NSETB) = ZERO DO 310 J = 1,NSETB ITEMP = IBASIS(J) JCOL = ABS(ITEMP) T1 = TWO T2 = TWO IF (ITEMP.LT.0) THEN BOU = ZERO ELSE BOU = BL(JCOL) ENDIF IF ((-BOU).NE.BIG) BOU = BOU/ABS(SCL(JCOL)) IF (RW(J).LE.BOU) T1 = (X(J)-BOU)/ (X(J)-RW(J)) BOU = BU(JCOL) IF (BOU.NE.BIG) BOU = BOU/ABS(SCL(JCOL)) IF (RW(J).GE.BOU) T2 = (BOU-X(J))/ (RW(J)-X(J)) C C If not, then compute a step length so that the variables remain C feasible. C IF (T1.LT.ALPHA) THEN ALPHA = T1 JDROP1 = J ENDIF C IF (T2.LT.BETA) THEN BETA = T2 JDROP2 = J ENDIF 310 CONTINUE C CONSTR = ALPHA .LT. TWO .OR. BETA .LT. TWO IF (.NOT.CONSTR) THEN C C Accept the candidate because it satisfies the stated bounds C on the variables. C CALL DCOPY(NSETB,RW,1,X,1) GO TO 580 ENDIF C C Take a step that is as large as possible with all variables C remaining feasible. C DO 330 J = 1,NSETB X(J) = X(J) + MIN(ALPHA,BETA)* (RW(J)-X(J)) 330 CONTINUE C IF (ALPHA.LE.BETA) THEN JDROP2 = 0 ELSE JDROP1 = 0 ENDIF C 340 IF (JDROP1+JDROP2.LE.0 .OR. NSETB.LE.0) GO TO 580 350 JDROP = JDROP1 + JDROP2 ITEMP = IBASIS(JDROP) JCOL = ABS(ITEMP) IF (JDROP2.GT.0) THEN C C Variable is at an upper bound. Subtract multiple of this C column from right hand side. C T = BU(JCOL) IF (ITEMP.GT.0) THEN BU(JCOL) = T - BL(JCOL) BL(JCOL) = -T ITEMP = -ITEMP SCL(JCOL) = -SCL(JCOL) DO 360 I = 1,JDROP W(I,JDROP) = -W(I,JDROP) 360 CONTINUE ELSE IBB(JCOL) = IBB(JCOL) + 1 IF (MOD(IBB(JCOL),2).EQ.0) T = -T ENDIF C C Variable is at a lower bound. C ELSE IF (ITEMP.LT.ZERO) THEN T = ZERO ELSE T = -BL(JCOL) BU(JCOL) = BU(JCOL) + T ITEMP = -ITEMP ENDIF ENDIF C CALL DAXPY(JDROP,T,W(1,JDROP),1,W(1,NCOLS+1),1) C C Move certain columns left to achieve upper Hessenberg form. C CALL DCOPY(JDROP,W(1,JDROP),1,RW,1) DO 370 J = JDROP+1,NSETB IBASIS(J-1) = IBASIS(J) X(J-1) = X(J) CALL DCOPY(J,W(1,J),1,W(1,J-1),1) 370 CONTINUE C IBASIS(NSETB) = ITEMP W(1,NSETB) = ZERO CALL DCOPY(MROWS-JDROP,W(1,NSETB),0,W(JDROP+1,NSETB),1) CALL DCOPY(JDROP,RW,1,W(1,NSETB),1) C C Transform the matrix from upper Hessenberg form to upper C triangular form. C NSETB = NSETB - 1 DO 390 I = JDROP,NSETB C C Look for small pivots and avoid mixing weighted and C nonweighted rows. C IF (I.EQ.MVAL) THEN T = ZERO DO 380 J = I,NSETB JCOL = ABS(IBASIS(J)) T1 = ABS(W(I,J)*SCL(JCOL)) IF (T1.GT.T) THEN JBIG = J T = T1 ENDIF 380 CONTINUE GO TO 400 ENDIF CALL DROTG(W(I,I),W(I+1,I),SC,SS) W(I+1,I) = ZERO CALL DROT(NCOLS-I+1,W(I,I+1),MDW,W(I+1,I+1),MDW,SC,SS) 390 CONTINUE GO TO 430 C C The triangularization is completed by giving up the Hessenberg C form and triangularizing a rectangular matrix. C 400 CALL DSWAP(MROWS,W(1,I),1,W(1,JBIG),1) CALL DSWAP(1,WW(I),1,WW(JBIG),1) CALL DSWAP(1,X(I),1,X(JBIG),1) ITEMP = IBASIS(I) IBASIS(I) = IBASIS(JBIG) IBASIS(JBIG) = ITEMP JBIG = I DO 420 J = JBIG,NSETB DO 410 I = J+1,MROWS CALL DROTG(W(J,J),W(I,J),SC,SS) W(I,J) = ZERO CALL DROT(NCOLS-J+1,W(J,J+1),MDW,W(I,J+1),MDW,SC,SS) 410 CONTINUE 420 CONTINUE C C See if the remaining coefficients are feasible. They should be C because of the way MIN(ALPHA,BETA) was chosen. Any that are not C feasible will be set to their bounds and appropriately translated. C 430 JDROP1 = 0 JDROP2 = 0 LGOPR = 2 GO TO 270 C C Find a variable to become non-active. C 120 IF (FOUND) THEN LGOPR = 1 GO TO 270 ENDIF C C Rescale and translate variables. C IGOPR = 2 130 CALL DCOPY(NSETB,X,1,RW,1) CALL DCOPY(NCOLS,ZERO,0,X,1) DO 140 J = 1,NSETB JCOL = ABS(IBASIS(J)) X(JCOL) = RW(J)*ABS(SCL(JCOL)) 140 CONTINUE C DO 150 J = 1,NCOLS IF (MOD(IBB(J),2).EQ.0) X(J) = BU(J) - X(J) 150 CONTINUE C DO 160 J = 1,NCOLS JCOL = IBASIS(J) IF (JCOL.LT.0) X(-JCOL) = BL(-JCOL) + X(-JCOL) 160 CONTINUE C DO 170 J = 1,NCOLS IF (SCL(J).LT.ZERO) X(J) = -X(J) 170 CONTINUE C I = MAX(NSETB,MVAL) RNORM = DNRM2(MROWS-I,W(INEXT(I),NCOLS+1),1) C IF (IGOPR.EQ.2) MODE = NSETB RETURN END