| [w,mdw,mrows,ncols,bl,bu,ind,iopt,x,rnorm,mode,rw,iw]=sbols(w,mdw,mrows,ncols,bl,bu,ind,iopt,x,rnorm,mode,rw,iw); |
function [w,mdw,mrows,ncols,bl,bu,ind,iopt,x,rnorm,mode,rw,iw]=sbols(w,mdw,mrows,ncols,bl,bu,ind,iopt,x,rnorm,mode,rw,iw);
persistent checkl firstCall gt100 gt200 i ibig igo inrows ip iscale j jp lds lenx liopt llb lliw llrw llx lmdw lndw locacc locdim lopt lp mnew nerr one sc ss xern1 xern2 xern3 xern4 zero ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(ibig), ibig=0; end;
if isempty(igo), igo=0; end;
if isempty(inrows), inrows=0; end;
if isempty(ip), ip=0; end;
if isempty(iscale), iscale=0; end;
if isempty(j), j=0; end;
if isempty(jp), jp=0; end;
if isempty(lds), lds=0; end;
if isempty(lenx), lenx=0; end;
if isempty(liopt), liopt=0; end;
if isempty(llb), llb=0; end;
if isempty(lliw), lliw=0; end;
if isempty(llrw), llrw=0; end;
if isempty(llx), llx=0; end;
if isempty(lmdw), lmdw=0; end;
if isempty(lndw), lndw=0; end;
if isempty(locacc), locacc=0; end;
if isempty(locdim), locdim=0; end;
if isempty(gt100), gt100=0; end;
if isempty(gt200), gt200=0; end;
%***BEGIN PROLOGUE SBOLS
%***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 SINGLE PRECISION (SBOLS-S, DBOLS-D)
%***KEYWORDS BOUNDS, CONSTRAINTS, INEQUALITY, LEAST SQUARES, LINEAR
%***AUTHOR Hanson, R. J., (SNLA)
%***DESCRIPTION
%
% 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 >= 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) >= BL(J).
% (the value of BU(J) is not used.)
% 2. For IND(J)=2, require X(J) <= BU(J).
% (the value of BL(J) is not used.)
% 3. For IND(J)=3, require X(J) >= BL(J) and
% X(J) <= 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) > 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 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 SBOLSM( ).
% 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 >=0 or ==-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 >= 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 < -1
% correspond to an abnormal completion of the subprogram. To
% understand the abnormal completion codes see below: ERROR
% MESSAGES for SBOLS( ). 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
% SBOLSM(). See the section ERROR MESSAGES for SBOLSM() in the
% documentation for SBOLSM().
%
% -----------
% 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 SBOLS( ). The value of
% MROWS is an OUTPUT variable when this option is used. Its value
% is always in the range 0 <= MROWS <= 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 SBOLS( ) 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 SBOLS()
% . 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 SBOLS() 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 SBOLS()
%
% 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 SBOLS()
%
% 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 SBOLS()
%
% 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 SBOLSM(). If LP is the processing pointer
% for IOPT(*),
%
% IOPT(LP)=5
% IOPT(LP+1)= Position in IOPT(*) where option array
% data for SBOLSM() begins.
% .
% CALL SBOLS()
%
% 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 SBOLS()
%
% 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 SBOLS()
% ----- -------- --- -------
%
% WARNING IN...
% SBOLS(). 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...
% SBOLS(). 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...
% SBOLS(). 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...
% SBOLS(). 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...
% SBOLS(). 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...
% SBOLS(). 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...
% SBOLS(). 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...
% SBOLS().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...
% SBOLS(). 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...
% SBOLS().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...
% SBOLS(). 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...
% SBOLS(). 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...
% SBOLS() 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...
% SBOLS() 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 ISAMAX, SBOLSM, SCOPY, SNRM2, SROT, SROTG, XERMSG
%***REVISION HISTORY (YYMMDD)
% 821220 DATE WRITTEN
if isempty(lopt), lopt=0; end;
if isempty(lp), lp=0; end;
if isempty(mnew), mnew=0; end;
if isempty(nerr), nerr=0; end;
% 861211 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900510 Convert XERRWV calls to XERMSG calls. (RWC)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE SBOLS
%
% SOLVE LINEAR LEAST SQUARES SYSTEM WITH BOUNDS ON
% SELECTED VARIABLES.
% REVISED 850329-1400
% REVISED YYMMDD-HHMM
% TO CHANGE THIS SUBPROGRAM FROM SINGLE TO doubleprecision BEGIN
% EDITING AT THE CARD 'C++'.
% CHANGE THIS SUBPROGRAM NAME TO DBOLS AND THE STRINGS
% /SCOPY/ TO /DCOPY/, /SBOL/ TO /DBOL/,
% /SNRM2/ TO /DNRM2/, /ISAMAX/ TO /IDAMAX/,
% /SROTG/ TO /DROTG/, /SROT/ TO /DROT/, /E0/ TO /D0/,
% /REAL / TO /doubleprecision/.
% ++
w_shape=size(w);w=reshape([w(:).',zeros(1,ceil(numel(w)./prod([mdw])).*prod([mdw])-numel(w))],mdw,[]);
bl_shape=size(bl);bl=reshape(bl,1,[]);
bu_shape=size(bu);bu=reshape(bu,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
rw_shape=size(rw);rw=reshape(rw,1,[]);
if isempty(sc), sc=0; end;
if isempty(ss), ss=0; end;
if isempty(one), one=0; end;
if isempty(zero), zero=0; end;
%
% THIS VARIABLE SHOULD REMAIN TYPE REAL.
ind_shape=size(ind);ind=reshape(ind,1,[]);
iopt_shape=size(iopt);iopt=reshape(iopt,1,[]);
iw_shape=size(iw);iw=reshape(iw,1,[]);
if isempty(checkl), checkl=false; end;
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern2), xern2=repmat(' ',1,8); end;
if isempty(xern3), xern3=repmat(' ',1,16); end;
if isempty(xern4), xern4=repmat(' ',1,16); end;
if firstCall, igo=[0]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT SBOLS
nerr = 0;
mode = 0;
while (1);
gt100=0;
gt200=0;
if( igo==0 )
% DO(CHECK VALIDITY OF INPUT DATA)
% PROCEDURE(CHECK VALIDITY OF INPUT DATA)
%
% SEE THAT MDW IS .GT.0. GROSS CHECK ONLY.
if( mdw<=0 )
xern1=sprintf(['%8i'], mdw);
% DO(RETURN TO USER PROGRAM UNIT)
xermsg('SLATEC','SBOLS',['MDW = ',[xern1,' MUST BE POSITIVE.']],2,1);
%
% SEE THAT NUMBER OF UNKNOWNS IS POSITIVE.
elseif( ncols<=0 ) ;
xern1=sprintf(['%8i'], ncols);
% DO(RETURN TO USER PROGRAM UNIT)
xermsg('SLATEC','SBOLS',['NCOLS = ',[xern1,' THE NO. OF VARIABLES MUST BE POSITIVE.']],3,1);
else;
%
% SEE THAT CONSTRAINT INDICATORS ARE ALL WELL-DEFINED.
for j = 1 : ncols;
if( ind(j)<1 || ind(j)>4 )
xern1=sprintf(['%8i'], j);
xern2=sprintf(['%8i'], ind(j));
xermsg('SLATEC','SBOLS',['IND(',[xern1,[') = ',[xern2,' MUST BE 1-4.']]]],4,1);
% DO(RETURN TO USER PROGRAM UNIT)
gt200=1;
break;
end;
end;
if(gt200==1)
break;
end;
%
% SEE THAT BOUNDS ARE CONSISTENT.
for j = 1 : ncols;
if( ind(j)==3 )
if( bl(j)>bu(j) )
xern1=sprintf(['%8i'], j);
xern3=sprintf([repmat('%15.6f',1,1)], bl(j));
xern4=sprintf([repmat('%15.6f',1,1)], bu(j));
xermsg('SLATEC','SBOLS',['BOUND BL(',[xern1,[') = ',[xern3,[' IS .GT. BU(',[xern1,[') = ',xern4]]]]]]],5,1);
% DO(RETURN TO USER PROGRAM UNIT)
gt200=1;
break;
end;
end;
end;
if(gt200==1)
break;
end;
% end PROCEDURE
% DO(PROCESS OPTION ARRAY)
% PROCEDURE(PROCESS OPTION ARRAY)
zero = 0.0e0;
one = 1.0e0;
checkl = false;
lenx = fix(ncols);
iscale = 1;
igo = 2;
lopt = 0;
lp = 0;
lds = 0;
while( true );
lp = fix(lp + lds);
ip = fix(iopt(lp+1));
jp = fix(abs(ip));
%
% TEST FOR NO MORE OPTIONS.
if( ip==99 )
if( lopt==0 )
lopt = fix(lp + 1);
end;
% end PROCEDURE
if( ~checkl )
gt100=1;
break;
end;
% DO(CHECK LENGTHS OF ARRAYS)
% PROCEDURE(CHECK LENGTHS OF ARRAYS)
%
% THIS FEATURE ALLOWS THE USER TO MAKE SURE THAT THE
% ARRAYS ARE LONG ENOUGH FOR THE INTENDED PROBLEM SIZE AND USE.
if( lmdw<mrows )
xern1=sprintf(['%8i'], lmdw);
xern2=sprintf(['%8i'], mrows);
xermsg('SLATEC','SBOLS',['THE ROW DIMENSION OF W() = ',[xern1,[' MUST BE .GE. THE NUMBER OF ROWS = ',xern2]]],11,1);
elseif( lndw<ncols+1 ) ;
xern1=sprintf(['%8i'], lndw);
xern2=sprintf(['%8i'], ncols + 1);
% DO(RETURN TO USER PROGRAM UNIT)
xermsg('SLATEC','SBOLS',['THE COLUMN DIMENSION OF W() = ',[xern1,[' MUST BE .GE. NCOLS+1 = ',xern2]]],12,1);
elseif( llb<ncols ) ;
xern1=sprintf(['%8i'], llb);
xern2=sprintf(['%8i'], ncols);
% DO(RETURN TO USER PROGRAM UNIT)
xermsg('SLATEC','SBOLS',['THE DIMENSIONS OF THE ARRAYS BL(), BU(), AND IND() = ',[xern1,[' MUST BE .GE. NCOLS = ',xern2]]],13,1);
elseif( llx<lenx ) ;
xern1=sprintf(['%8i'], llx);
xern2=sprintf(['%8i'], lenx);
% DO(RETURN TO USER PROGRAM UNIT)
xermsg('SLATEC','SBOLS',['THE DIMENSION OF X() = ',[xern1,[' MUST BE .GE. THE REQUIRED LENGTH = ',xern2]]],14,1);
elseif( llrw<5.*ncols ) ;
xern1=sprintf(['%8i'], llrw);
xern2=sprintf(['%8i'], 5.*ncols);
% DO(RETURN TO USER PROGRAM UNIT)
xermsg('SLATEC','SBOLS',['THE DIMENSION OF RW() = ',[xern1,[' MUST BE .GE. 5*NCOLS = ',xern2]]],15,1);
elseif( lliw<2.*ncols ) ;
xern1=sprintf(['%8i'], lliw);
xern2=sprintf(['%8i'], 2.*ncols);
% DO(RETURN TO USER PROGRAM UNIT)
xermsg('SLATEC','SBOLS',['THE DIMENSION OF IW() = ',[xern1,[' MUST BE .GE. 2*NCOLS = ',xern2]]],16,1);
elseif( liopt<lp+1 ) ;
% end PROCEDURE
xern1=sprintf(['%8i'], liopt);
xern2=sprintf(['%8i'], lp + 1);
% DO(RETURN TO USER PROGRAM UNIT)
xermsg('SLATEC','SBOLS',['THE DIMENSION OF IOPT() = ',[xern1,[' MUST BE .GE. THE REQD. LEN = ',xern2]]],17,1);
else;
gt100=1;
break;
end;
gt200=1;
break;
elseif( jp==99 ) ;
lds = 1;
elseif( jp==1 ) ;
if( ip>0 )
%
% SET UP DIRECTION FLAG, ROW STACKING POINTER
% LOCATION, AND LOCATION FOR NUMBER OF NEW ROWS.
locacc = fix(lp + 2);
%
% IOPT(LOCACC-1)=OPTION NUMBER FOR SEQ. ACCUMULATION.
% CONTENTS.. IOPT(LOCACC )=USER DIRECTION FLAG, 1 OR 2.
% IOPT(LOCACC+1)=ROW STACKING POINTER.
% IOPT(LOCACC+2)=NUMBER OF NEW ROWS TO PROCESS.
% USER ACTION WITH THIS OPTION..
% (SET UP OPTION DATA FOR SEQ. ACCUMULATION IN IOPT(*).
% MUST ALSO START PROCESS WITH IOPT(LOCACC)=1.)
% (MOVE BLOCK OF EQUATIONS INTO W(*,*) STARTING AT FIRST
% ROW OF W(*,*). SET IOPT(LOCACC+2)=NO. OF ROWS IN BLOCK.)
% LOOP
% CALL SBOLS()
%
% IF(IOPT(LOCACC) .EQ. 1) THEN
% STACK EQUAS., STARTING AT ROW IOPT(LOCACC+1),
% INTO W(*,*).
% SET IOPT(LOCACC+2)=NO. OF EQUAS.
% IF LAST BLOCK OF EQUAS., SET IOPT(LOCACC)=2.
% ELSE IF IOPT(LOCACC) .EQ. 2) THEN
% (PROCESS IS OVER. EXIT LOOP.)
% ELSE
% (ERROR CONDITION. SHOULD NOT HAPPEN.)
% endif
% end LOOP
% SET IOPT(LOCACC-1)=-OPTION NUMBER FOR SEQ. ACCUMULATION.
% CALL SBOLS( )
iopt(locacc+1) = 1;
igo = 1;
end;
lds = 4;
elseif( jp==2 ) ;
if( ip>0 )
%
% GET ACTUAL LENGTHS OF ARRAYS FOR CHECKING AGAINST NEEDS.
locdim = fix(lp + 2);
%
% LMDW.GE.MROWS
% LNDW.GE.NCOLS+1
% LLB .GE.NCOLS
% LLX .GE.NCOLS+EXTRA REQD. IN OPTIONS.
% LLRW.GE.5*NCOLS
% LLIW.GE.2*NCOLS
% LIOP.GE. AMOUNT REQD. FOR IOPTION ARRAY.
lmdw = fix(iopt(locdim));
lndw = fix(iopt(locdim+1));
llb = fix(iopt(locdim+2));
llx = fix(iopt(locdim+3));
llrw = fix(iopt(locdim+4));
lliw = fix(iopt(locdim+5));
liopt = fix(iopt(locdim+6));
checkl = true;
end;
lds = 8;
%
% OPTION TO MODIFY THE COLUMN SCALING.
elseif( jp==3 ) ;
if( ip>0 )
iscale = fix(iopt(lp+2));
%
% SEE THAT ISCALE IS 1 THRU 3.
if( iscale<1 || iscale>3 )
xern1=sprintf(['%8i'], iscale);
xermsg('SLATEC','SBOLS',['ISCALE OPTION = ',[xern1,' MUST BE 1-3']],7,1);
% DO(RETURN TO USER PROGRAM UNIT)
gt200=1;
break;
end;
end;
% CYCLE FOREVER
lds = 2;
%
% IN THIS OPTION THE USER HAS PROVIDED SCALING. THE
% SCALE FACTORS FOR THE COLUMNS BEGIN IN X(NCOLS+IOPT(LP+2)).
elseif( jp==4 ) ;
if( ip>0 )
iscale = 4;
if( iopt(lp+2)<=0 )
xern1=sprintf(['%8i'], iopt(lp+2));
xermsg('SLATEC','SBOLS',['OFFSET PAST X(NCOLS) (',[xern1,') FOR USER-PROVIDED COLUMN SCALING MUST BE POSITIVE.']],8,1);
% DO(RETURN TO USER PROGRAM UNIT)
gt200=1;
break;
else;
[ncols,x(sub2ind(size(x),max(ncols+iopt(lp+2),1)):end),dumvar3,rw]=scopy(ncols,x(sub2ind(size(x),max(ncols+iopt(lp+2),1)):end),1,rw,1);
lenx = fix(lenx + ncols);
for j = 1 : ncols;
if( rw(j)<=zero )
xern1=sprintf(['%8i'], j);
xern3=sprintf([repmat('%15.6f',1,1)], rw(j));
xermsg('SLATEC','SBOLS',['EACH PROVIDED COLUMN SCALE FACTOR ',['MUST BE POSITIVE.$$COMPONENT ',[xern1,[' NOW = ',xern3]]]],9,1);
% DO(RETURN TO USER PROGRAM UNIT)
gt200=1;
break;
end;
end;
if(gt200==1)
break;
end;
end;
end;
% CYCLE FOREVER
lds = 2;
%
% IN THIS OPTION AN OPTION ARRAY IS PROVIDED TO SBOLSM().
elseif( jp==5 ) ;
if( ip>0 )
lopt = fix(iopt(lp+2));
end;
% CYCLE FOREVER
lds = 2;
%
% THIS OPTION USES THE NEXT LOC OF IOPT(*) AS AN
% INCREMENT TO SKIP.
elseif( jp~=6 ) ;
break;
elseif( ip>0 ) ;
%
% NO VALID OPTION NUMBER WAS NOTED. THIS IS AN ERROR CONDITION.
lp = fix(iopt(lp+2) - 1);
lds = 0;
else;
lds = 2;
% CYCLE FOREVER
end;
end;
if(gt200==1 || gt100==1)
break;
end;
xern1=sprintf(['%8i'], jp);
% DO(RETURN TO USER PROGRAM UNIT)
xermsg('SLATEC','SBOLS',['THE OPTION NUMBER = ',[xern1,' IS NOT DEFINED.']],6,1);
end;
gt200=1;
break;
end;
break;
end;
while(gt200==0);
if( igo==1 )
%
% GO BACK TO THE USER FOR ACCUMULATION OF LEAST SQUARES
% EQUATIONS AND DIRECTIONS TO QUIT PROCESSING.
% CASE 1
% DO(ACCUMULATE LEAST SQUARES EQUATIONS)
% PROCEDURE(ACCUMULATE LEAST SQUARES EQUATIONS)
mrows = fix(iopt(locacc+1) - 1);
inrows = fix(iopt(locacc+2));
mnew = fix(mrows + inrows);
if( mnew<0 || mnew>mdw )
xern1=sprintf(['%8i'], mnew);
xern2=sprintf(['%8i'], mdw);
xermsg('SLATEC','SBOLS',['NO. OF ROWS = ',[xern1,[' MUST BE .GE. 0 .AND. .LE. MDW = ',xern2]]],10,1);
% DO(RETURN TO USER PROGRAM UNIT)
break;
else;
for j = 1 : min(ncols+1,mnew);
for i = mnew : -1: max(mrows,j) + 1 ;
ibig = fix(isamax(i-j,w(sub2ind(size(w),j,j):end),1) + j - 1);
%
% PIVOT FOR INCREASED STABILITY.
[w(ibig,j),w(i,j),sc,ss]=srotg(w(ibig,j),w(i,j),sc,ss);
mdw_orig=mdw; [dumvar1,dumvar2,mdw,dumvar4,dumvar5,sc,ss]=srot(ncols+1-j,w(sub2ind(size(w),ibig,j+1):end),mdw,w(sub2ind(size(w),i,j+1):end),mdw,sc,ss); mdw(dumvar5~=mdw_orig)=dumvar5(dumvar5~=mdw_orig); w(sub2ind(size(w),ibig,j+1):end)=dumvar2; w(sub2ind(size(w),i,j+1):end)=dumvar4;
w(i,j) = zero;
end; i = fix(max(mrows,j) + 1 -1);
end; j = fix(min(ncols+1,mnew)+1);
mrows = fix(min(ncols+1,mnew));
iopt(locacc+1) = fix(mrows + 1);
igo = fix(iopt(locacc));
% end PROCEDURE
if( igo==2 )
igo = 0;
end;
end;
elseif( igo==2 ) ;
% CASE 2
% DO(INITIALIZE VARIABLES AND DATA VALUES)
% PROCEDURE(INITIALIZE VARIABLES AND DATA VALUES)
for j = 1 : ncols;
if( iscale==1 )
% CASE 1
%
% THIS IS THE NOMINAL SCALING. EACH NONZERO
% COL. HAS MAX. NORM EQUAL TO ONE.
[ibig ,mrows,w(sub2ind(size(w),1,j):end)]=isamax(mrows,w(sub2ind(size(w),1,j):end),1);
rw(j) = abs(w(ibig,j));
if( rw(j)==zero )
rw(j) = one;
else;
rw(j) = one./rw(j);
end;
elseif( iscale==2 ) ;
% CASE 2
%
% THIS CHOICE OF SCALING MAKES EACH NONZERO COLUMN
% HAVE EUCLIDEAN LENGTH EQUAL TO ONE.
[rw(j) ,mrows,w(sub2ind(size(w),1,j):end)]=snrm2(mrows,w(sub2ind(size(w),1,j):end),1);
if( rw(j)==zero )
rw(j) = one;
else;
rw(j) = one./rw(j);
end;
elseif( iscale==3 ) ;
% CASE 3
%
% THIS CASE EFFECTIVELY SUPPRESSES SCALING BY SETTING
% THE SCALING MATRIX TO THE IDENTITY MATRIX.
rw(1) = one;
rw_orig=rw; [ncols,rw,dumvar3,dumvar4]=scopy(ncols,rw,0,rw,1); rw(dumvar4~=rw_orig)=dumvar4(dumvar4~=rw_orig);
break;
elseif( iscale==4 ) ;
break;
end;
end;
% end PROCEDURE
% DO(SOLVE BOUNDED LEAST SQUARES PROBLEM)
% PROCEDURE(SOLVE BOUNDED LEAST SQUARES PROBLEM)
%
% INITIALIZE IBASIS(*), J=1,NCOLS, AND IBB(*), J=1,NCOLS,
% TO =J,AND =1, FOR USE IN SBOLSM( ).
for j = 1 : ncols;
iw(j) = fix(j);
iw(j+ncols) = 1;
rw(3.*ncols+j) = bl(j);
rw(4.*ncols+j) = bu(j);
end; j = fix(ncols+1);
[w,mdw,mrows,ncols,dumvar5,dumvar6,ind,iopt(sub2ind(size(iopt),max(lopt,1)):end),x,rnorm,mode,dumvar12,dumvar13,rw,iw,dumvar16]=sbolsm(w,mdw,mrows,ncols,rw(sub2ind(size(rw),max(3.*ncols+1,1)):end),rw(sub2ind(size(rw),max(4.*ncols+1,1)):end),ind,iopt(sub2ind(size(iopt),max(lopt,1)):end),x,rnorm,mode,rw(sub2ind(size(rw),max(ncols+1,1)):end),rw(sub2ind(size(rw),max(2.*ncols+1,1)):end),rw,iw,iw(sub2ind(size(iw),max(ncols+1,1)):end)); dumvar5i=find((rw(sub2ind(size(rw),max(3.*ncols+1,1)):end))~=(dumvar5));dumvar6i=find((rw(sub2ind(size(rw),max(4.*ncols+1,1)):end))~=(dumvar6));dumvar12i=find((rw(sub2ind(size(rw),max(ncols+1,1)):end))~=(dumvar12));dumvar13i=find((rw(sub2ind(size(rw),max(2.*ncols+1,1)):end))~=(dumvar13));dumvar16i=find((iw(sub2ind(size(iw),max(ncols+1,1)):end))~=(dumvar16)); rw(3.*ncols+1-1+dumvar5i)=dumvar5(dumvar5i); rw(4.*ncols+1-1+dumvar6i)=dumvar6(dumvar6i); rw(ncols+1-1+dumvar12i)=dumvar12(dumvar12i); rw(2.*ncols+1-1+dumvar13i)=dumvar13(dumvar13i); iw(ncols+1-1+dumvar16i)=dumvar16(dumvar16i);
% end PROCEDURE
igo = 0;
end;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
bl_shape=zeros(bl_shape);bl_shape(:)=bl(1:numel(bl_shape));bl=bl_shape;
bu_shape=zeros(bu_shape);bu_shape(:)=bu(1:numel(bu_shape));bu=bu_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
rw_shape=zeros(rw_shape);rw_shape(:)=rw(1:numel(rw_shape));rw=rw_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
iopt_shape=zeros(iopt_shape);iopt_shape(:)=iopt(1:numel(iopt_shape));iopt=iopt_shape;
iw_shape=zeros(iw_shape);iw_shape(:)=iw(1:numel(iw_shape));iw=iw_shape;
return;
% PROCEDURE(RETURN TO USER PROGRAM UNIT)
break;
end;
if( mode>=0 )
mode = fix(-nerr);
end;
igo = 0;
% end PROCEDURE
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
bl_shape=zeros(bl_shape);bl_shape(:)=bl(1:numel(bl_shape));bl=bl_shape;
bu_shape=zeros(bu_shape);bu_shape(:)=bu(1:numel(bu_shape));bu=bu_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
rw_shape=zeros(rw_shape);rw_shape(:)=rw(1:numel(rw_shape));rw=rw_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
iopt_shape=zeros(iopt_shape);iopt_shape(:)=iopt(1:numel(iopt_shape));iopt=iopt_shape;
iw_shape=zeros(iw_shape);iw_shape(:)=iw(1:numel(iw_shape));iw=iw_shape;
end %subroutine sbols
%DECK SBOLSM
|
|
