SUBROUTINE PDA_DBOLS (W, MDW, MROWS, NCOLS, BL, BU, IND, IOPT, X,
+ RNORM, MODE, RW, IW, STATUS)
***BEGIN PROLOGUE PDA_DBOLS
***PURPOSE Solve the problem
E*X = F (in the least squares sense)
with bounds on selected X values.
***LIBRARY SLATEC
***CATEGORY K1A2A, G2E, G2H1, G2H2
***TYPE DOUBLE PRECISION (SBOLS-S, PDA_DBOLS-D)
***KEYWORDS BOUNDS, CONSTRAINTS, INEQUALITY, LEAST SQUARES, LINEAR
***AUTHOR Hanson, R. J., (SNLA)
***DESCRIPTION
**** All INPUT and OUTPUT real variables are DOUBLE PRECISION ****
The user must have dimension statements of the form:
DIMENSION W(MDW,NCOLS+1), BL(NCOLS), BU(NCOLS),
* X(NCOLS+NX), RW(5*NCOLS)
INTEGER IND(NCOLS), IOPT(1+NI), IW(2*NCOLS)
(Here NX=number of extra locations required for option 4; NX=0
for no options; NX=NCOLS if this option is in use. Here NI=number
of extra locations required for options 1-6; NI=0 for no
options.)
INPUT
-----
--------------------
W(MDW,*),MROWS,NCOLS
--------------------
The array W(*,*) contains the matrix [E:F] on entry. The matrix
[E:F] has MROWS 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. MROWS.
Other values of MDW are errors. The values of MROWS and NCOLS
must be positive. Other values are errors. There is an exception
to this when using option 1 for accumulation of blocks of
equations. In that case MROWS is an OUTPUT variable ONLY, and the
matrix data for [E:F] is placed in W(*,*), one block of rows at a
time. MROWS contains the number of rows in the matrix after
triangularizing several blocks of equations. This is an OUTPUT
parameter ONLY when option 1 is used. See IOPT(*) CONTENTS
for details about option 1.
------------------
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).
(the value of BU(J) is not used.)
2. For IND(J)=2, require X(J) .le. BU(J).
(the value of BL(J) is not used.)
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(J) and BU(J) are not used.)
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 PDA_DBOLSM( ). 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 Return to user for accumulation of blocks
of least squares equations.
2 Check lengths of all arrays used in the
subprogram.
3 Standard scaling of the data matrix, E.
4 User provides column scaling for matrix E.
5 Provide option array to the low-level
subprogram PDA_DBOLSM( ).
6 Move the IOPT(*) processing pointer.
99 No more options to change.
----
X(*)
----
This array is used to pass data associated with option 4. Ignore
this parameter if this option is not used. Otherwise see below:
IOPT(*) CONTENTS.
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 cases -37,-36,...,-22, or -17,...,-2. Values .lt. -1
correspond to an abnormal completion of the subprogram. To
understand the abnormal completion codes see below: ERROR
MESSAGES for PDA_DBOLS( ). AN approximate solution will be returned
to the user only when max. iterations is reached, MODE=-22.
Values for MODE=-37,...,-22 come from the low-level subprogram
PDA_DBOLSM(). See the section ERROR MESSAGES for PDA_DBOLSM() in the
documentation for PDA_DBOLSM().
------
STATUS
------
Returned error status.
The status must be zero on entry. This
routine does not check the status on entry.
-----------
RW(*),IW(*)
-----------
These are working arrays with 5*NCOLS and 2*NCOLS entries.
(normally the user can ignore the contents of these arrays,
but they must be dimensioned properly.)
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. This value is updated as each option is 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 option 4 uses the
locations X(NCOLS+IOFF),...,X(NCOLS+IOFF+NCOLS-1) for passing
scaling data. The user must manage the allocation of these
locations.
1
-
This option allows the user to solve problems with a large number
of rows compared to the number of variables. The idea is that the
subprogram returns to the user (perhaps many times) and receives
new least squares equations from the calling program unit.
Eventually the user signals "that's all" and then computes the
solution with one final call to subprogram PDA_DBOLS( ). The value of
MROWS is an OUTPUT variable when this option is used. Its value
is always in the range 0 .le. MROWS .le. NCOLS+1. It is equal to
the number of rows after the triangularization of the entire set
of equations. If LP is the processing pointer for IOPT(*), the
usage for the sequential processing of blocks of equations is
IOPT(LP)=1
Move block of equations to W(*,*) starting at
the first row of W(*,*).
IOPT(LP+3)=# of rows in the block; user defined
The user now calls PDA_DBOLS( ) in a loop. The value of IOPT(LP+1)
directs the user's action. The value of IOPT(LP+2) points to
where the subsequent rows are to be placed in W(*,*).
.<LOOP
. CALL PDA_DBOLS()
. IF(IOPT(LP+1) .EQ. 1) THEN
. IOPT(LP+3)=# OF ROWS IN THE NEW BLOCK; USER DEFINED
. PLACE NEW BLOCK OF IOPT(LP+3) ROWS IN
. W(*,*) STARTING AT ROW IOPT(LP+2).
.
. IF( THIS IS THE LAST BLOCK OF EQUATIONS ) THEN
. IOPT(LP+1)=2
.<------CYCLE LOOP
. ELSE IF (IOPT(LP+1) .EQ. 2) THEN
<-------EXIT LOOP SOLUTION COMPUTED IF MODE .GE. 0
. ELSE
. ERROR CONDITION; SHOULD NOT HAPPEN.
.<END LOOP
Use of this option adds 4 to the required length of IOPT(*).
2
-
This option is useful for checking the lengths of all arrays used
by PDA_DBOLS() against their actual requirements for this problem.
The idea is simple: the user's program unit passes the declared
dimension information of the arrays. These values are compared
against the problem-dependent needs within the subprogram. If any
of the dimensions are too small an error message is printed and a
negative value of MODE is returned, -11 to -17. The printed error
message tells how long the dimension should be. If LP is the
processing pointer for IOPT(*),
IOPT(LP)=2
IOPT(LP+1)=Row dimension of W(*,*)
IOPT(LP+2)=Col. dimension of W(*,*)
IOPT(LP+3)=Dimensions of BL(*),BU(*),IND(*)
IOPT(LP+4)=Dimension of X(*)
IOPT(LP+5)=Dimension of RW(*)
IOPT(LP+6)=Dimension of IW(*)
IOPT(LP+7)=Dimension of IOPT(*)
.
CALL PDA_DBOLS()
Use of this option adds 8 to the required length of IOPT(*).
3
-
This option changes the type of scaling for the data matrix E.
Nominally each nonzero column of E is scaled so that the
magnitude of its largest entry is equal to the value ONE. If LP
is the processing pointer for IOPT(*),
IOPT(LP)=3
IOPT(LP+1)=1,2 or 3
1= Nominal scaling as noted;
2= Each nonzero column scaled to have length ONE;
3= Identity scaling; scaling effectively suppressed.
.
CALL PDA_DBOLS()
Use of this option adds 2 to the required length of IOPT(*).
4
-
This option allows the user to provide arbitrary (positive)
column scaling for the matrix E. If LP is the processing pointer
for IOPT(*),
IOPT(LP)=4
IOPT(LP+1)=IOFF
X(NCOLS+IOFF),...,X(NCOLS+IOFF+NCOLS-1)
= Positive scale factors for cols. of E.
.
CALL PDA_DBOLS()
Use of this option adds 2 to the required length of IOPT(*) and
NCOLS to the required length of X(*).
5
-
This option allows the user to provide an option array to the
low-level subprogram PDA_DBOLSM(). If LP is the processing pointer
for IOPT(*),
IOPT(LP)=5
IOPT(LP+1)= Position in IOPT(*) where option array
data for PDA_DBOLSM() begins.
.
CALL PDA_DBOLS()
Use of this option adds 2 to the required length of IOPT(*).
6
-
Move the processing pointer (either forward or backward) to the
location IOPT(LP+1). The processing point is moved to entry
LP+2 of IOPT(*) if the option is left with -6 in IOPT(LP). For
example to skip over locations 3,...,NCOLS+2 of IOPT(*),
IOPT(1)=6
IOPT(2)=NCOLS+3
(IOPT(I), I=3,...,NCOLS+2 are not defined here.)
IOPT(NCOLS+3)=99
CALL PDA_DBOLS()
CAUTION: Misuse of this option can yield some very hard
-to-find bugs. Use it with care.
99
--
There are no more options to change.
Only option numbers -99, -6,-5,...,-1, 1,2,...,6, and 99 are
permitted. Other values are errors. Options -99,-1,...,-6 mean
that the respective options 99,1,...,6 are left at their default
values. An example is the option to modify the (rank) tolerance:
IOPT(1)=-3 Option is recognized but not changed
IOPT(2)=2 Scale nonzero cols. to have length ONE
IOPT(3)=99
ERROR MESSAGES for PDA_DBOLS()
----- -------- --- -------
WARNING IN...
PDA_DBOLS(). MDW=(I1) MUST BE POSITIVE.
IN ABOVE MESSAGE, I1= 0
ERROR NUMBER = 2
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS(). NCOLS=(I1) THE NO. OF VARIABLES MUST BE POSITIVE.
IN ABOVE MESSAGE, I1= 0
ERROR NUMBER = 3
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS(). FOR J=(I1), IND(J)=(I2) MUST BE 1-4.
IN ABOVE MESSAGE, I1= 1
IN ABOVE MESSAGE, I2= 0
ERROR NUMBER = 4
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS(). FOR J=(I1), BOUND BL(J)=(R1) IS .GT. BU(J)=(R2).
IN ABOVE MESSAGE, I1= 1
IN ABOVE MESSAGE, R1= 0.
IN ABOVE MESSAGE, R2= ABOVE MESSAGE, I1= 0
ERROR NUMBER = 6
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS(). ISCALE OPTION=(I1) MUST BE 1-3.
IN ABOVE MESSAGE, I1= 0
ERROR NUMBER = 7
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS(). OFFSET PAST X(NCOLS) (I1) FOR USER-PROVIDED COLUMN SCALING
MUST BE POSITIVE.
IN ABOVE MESSAGE, I1= 0
ERROR NUMBER = 8
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS(). EACH PROVIDED COL. SCALE FACTOR MUST BE POSITIVE.
COMPONENT (I1) NOW = (R1).
IN ABOVE MESSAGE, I1= ND. .LE. MDW=(I2).
IN ABOVE MESSAGE, I1= 1
IN ABOVE MESSAGE, I2= 0
ERROR NUMBER = 10
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS().THE ROW DIMENSION OF W(,)=(I1) MUST BE .GE.THE NUMBER OF ROWS=
(I2).
IN ABOVE MESSAGE, I1= 0
IN ABOVE MESSAGE, I2= 1
ERROR NUMBER = 11
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS(). THE COLUMN DIMENSION OF W(,)=(I1) MUST BE .GE. NCOLS+1=(I2).
IN ABOVE MESSAGE, I1= 0
IN ABOVE MESSAGE, I2= 2
ERROR NUMBER = 12
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS().THE DIMENSIONS OF THE ARRAYS BL(),BU(), AND IND()=(I1) MUST BE
.GE. NCOLS=(I2).
IN ABOVE MESSAGE, I1= 0
IN ABOVE MESSAGE, I2= 1
ERROR NUMBER = 13
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS(). THE DIMENSION OF X()=(I1) MUST BE .GE. THE REQD. LENGTH=(I2).
IN ABOVE MESSAGE, I1= 0
IN ABOVE MESSAGE, I2= 2
ERROR NUMBER = 14
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS(). THE DIMENSION OF RW()=(I1) MUST BE .GE. 5*NCOLS=(I2).
IN ABOVE MESSAGE, I1= 0
IN ABOVE MESSAGE, I2= 3
ERROR NUMBER = 15
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS() THE DIMENSION OF IW()=(I1) MUST BE .GE. 2*NCOLS=(I2).
IN ABOVE MESSAGE, I1= 0
IN ABOVE MESSAGE, I2= 2
ERROR NUMBER = 16
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
WARNING IN...
PDA_DBOLS() THE DIMENSION OF IOPT()=(I1) MUST BE .GE. THE REQD. LEN.=(I2).
IN ABOVE MESSAGE, I1= 0
IN ABOVE MESSAGE, I2= 1
ERROR NUMBER = 17
(NORMALLY A RETURN TO THE USER TAKES PLACE FOLLOWING THIS MESSAGE.)
***REFERENCES R. J. Hanson, Linear least squares with bounds and
linear constraints, Report SAND82-1517, Sandia
Laboratories, August 1982.
***ROUTINES CALLED PDA_DBOLSM, PDA_DCOPY, PDA_DNRM2, PDA_DROT, PDA_DROTG,
PDA_IDAMAX, PDA_XERMSG
***REVISION HISTORY (YYMMDD)
821220 DATE WRITTEN
891006 Cosmetic changes to prologue. (WRB)
891006 REVISION DATE from Version 3.2
891214 Prologue converted to Version 4.0 format. (BAB)
900510 Convert XERRWV calls to PDA_XERMSG calls. (RWC)
920501 Reformatted the REFERENCES section. (WRB)
950404 Implement status. (HME)
***END PROLOGUE PDA_DBOLS
PDA [1ex