| [w,mdw,minput,ncols,bl,bu,ind,iopt,x,rnorm,mode,rw,ww,scl,ibasis,ibb]=dbolsm(w,mdw,minput,ncols,bl,bu,ind,iopt,x,rnorm,mode,rw,ww,scl,ibasis,ibb); |
function [w,mdw,minput,ncols,bl,bu,ind,iopt,x,rnorm,mode,rw,ww,scl,ibasis,ibb]=dbolsm(w,mdw,minput,ncols,bl,bu,ind,iopt,x,rnorm,mode,rw,ww,scl,ibasis,ibb);
persistent alpha beta big bou cl1 cl2 cl3 colabv colblo constr fac found gt i idum igo igopr ii ioff ip iprint itemp iter itmax j jbig jcol jdrop jdrop1 jdrop2 jlarge jmag jp lds lgopr lp mrows mval nsetb one sc ss t t1 t2 tolind tolsze two wbig wlarge wmag wt xern1 xern2 xern3 xern4 xnew zero ;
if isempty(i), i=0; end;
if isempty(idum), idum=0; end;
if isempty(igopr), igopr=0; end;
% inext= @(idum) min(idum+1,mrows);integer :: inext ;
if isempty(ioff), ioff=0; end;
if isempty(ip), ip=0; end;
if isempty(iprint), iprint=0; end;
if isempty(itemp), itemp=0; end;
if isempty(iter), iter=0; end;
if isempty(itmax), itmax=0; end;
if isempty(j), j=0; end;
if isempty(jbig), jbig=0; end;
if isempty(jcol), jcol=0; end;
if isempty(jdrop), jdrop=0; end;
if isempty(jdrop1), jdrop1=0; end;
if isempty(jdrop2), jdrop2=0; end;
if isempty(jlarge), jlarge=0; end;
if isempty(jmag), jmag=0; end;
if isempty(jp), jp=0; end;
if isempty(lds), lds=0; end;
if isempty(ii), ii=0; end;
if isempty(igo), igo=0; end;
if isempty(gt), gt=zeros(1,6); end;
%***BEGIN PROLOGUE DBOLSM
%***SUBSIDIARY
%***PURPOSE Subsidiary to DBOCLS and DBOLS
%***LIBRARY SLATEC
%***TYPE doubleprecision (SBOLSM-S, DBOLSM-D)
%***AUTHOR (UNKNOWN)
%***DESCRIPTION
%
% **** doubleprecision Version of SBOLSM ****
% **** All INPUT and OUTPUT real variables are doubleprecision ****
%
% 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 >= 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) >= BL(J).
% 2. For IND(J)=2, require X(J) <= BU(J).
% 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(*),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) > BU(J) is an
% error.
%
% -------
% IOPT(*)
% -------
% This is the array where the user can specify nonstandard options
% for 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 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) <
% 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 >= 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 (>= 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 < -1 correspond
% to an abnormal completion of the subprogram. To understand the
% abnormal completion codes see below: ERROR MESSAGES for DBOLSM
% 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 >= 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 DBOLSM
%
% 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
% <= TOL* magnitude of the column in components strictly above
% 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 DBOLSM
%
% Generally, if LP is the processing pointer for IOPT(*),
%
% X(NCOLS+IOFF)=TOL
% IOPT(LP)=2
% IOPT(LP+1)=IOFF
% .
% CALL DBOLSM
%
% 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 <= 0 is an error. A value of TOL <= 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 >=((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 DBOLSM
%
% Generally, if LP is the processing pointer for IOPT(*),
%
% X(NCOLS+IOFF)=BLOWUP
% IOPT(LP)=3
% IOPT(LP+1)=IOFF
% .
% CALL DBOLSM
%
% 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 <= 0 is an error. A value of BLOWUP <= 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 DBOLSM
%
% Generally, if LP is the processing pointer for IOPT(*),
%
% IOPT(LP)=4
% IOPT(LP+1)=ITMAX
% .
% CALL DBOLSM
%
% The value of ITMAX must be > 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 > 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 DBOLSM
%
% Generally, if LP is the processing pointer for IOPT(*),
%
% X(NCOLS+IOFF)=FAC
% IOPT(LP)=5
% IOPT(LP+1)=IOFF
% .
% CALL DBOLSM
%
% 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 DBOLSM
% ----- -------- --- ---------
% -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 DBOCLS, DBOLS
%***ROUTINES CALLED D1MACH, DAXPY, DCOPY, DDOT, DMOUT, DNRM2, DROT,
% DROTG, DSWAP, DVOUT, IVOUT, 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)
%***end PROLOGUE DBOLSM
%
% PURPOSE
% -------
% THIS IS THE MAIN SUBPROGRAM THAT SOLVES THE BOUNDED
if isempty(lgopr), lgopr=0; end;
if isempty(lp), lp=0; end;
if isempty(mrows), mrows=0; end;
if isempty(mval), mval=0; end;
if isempty(nsetb), nsetb=0; end;
% LEAST SQUARES PROBLEM. THE PROBLEM SOLVED HERE IS:
%
% SOLVE E*X = F (LEAST SQUARES SENSE)
% WITH BOUNDS ON SELECTED X VALUES.
%
% TO CHANGE THIS SUBPROGRAM FROM SINGLE TO doubleprecision BEGIN
% EDITING AT THE CARD 'C++'.
% CHANGE THE SUBPROGRAM NAME TO DBOLSM AND THE STRINGS
% /SAXPY/ TO /DAXPY/, /SCOPY/ TO /DCOPY/,
% /SDOT/ TO /DDOT/, /SNRM2/ TO /DNRM2/,
% /SROT/ TO /DROT/, /SROTG/ TO /DROTG/, /R1MACH/ TO /D1MACH/,
% /SVOUT/ TO /DVOUT/, /SMOUT/ TO /DMOUT/,
% /SSWAP/ TO /DSWAP/, /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,[]);
ww_shape=size(ww);ww=reshape(ww,1,[]);
scl_shape=size(scl);scl=reshape(scl,1,[]);
if isempty(alpha), alpha=0; end;
if isempty(beta), beta=0; end;
if isempty(bou), bou=0; end;
if isempty(colabv), colabv=0; end;
if isempty(colblo), colblo=0; end;
if isempty(cl1), cl1=0; end;
if isempty(cl2), cl2=0; end;
if isempty(cl3), cl3=0; end;
if isempty(big), big=0; end;
if isempty(fac), fac=0; end;
if isempty(sc), sc=0; end;
if isempty(ss), ss=0; end;
if isempty(t), t=0; end;
if isempty(tolind), tolind=0; end;
if isempty(wt), wt=0; end;
if isempty(t1), t1=0; end;
if isempty(t2), t2=0; end;
if isempty(wbig), wbig=0; end;
if isempty(wlarge), wlarge=0; end;
if isempty(wmag), wmag=0; end;
if isempty(xnew), xnew=0; end;
if isempty(tolsze), tolsze=0; end;
ibasis_shape=size(ibasis);ibasis=reshape(ibasis,1,[]);
ibb_shape=size(ibb);ibb=reshape(ibb,1,[]);
ind_shape=size(ind);ind=reshape(ind,1,[]);
iopt_shape=size(iopt);iopt=reshape(iopt,1,[]);
if isempty(found), found=false; end;
if isempty(constr), constr=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 isempty(zero), zero=0.0d0; end;
if isempty(one), one=1.0d0; end;
if isempty(two), two=2.0d0 ; end;
%
inext= @(idum) min(idum+1,mrows);
%***FIRST EXECUTABLE STATEMENT DBOLSM
%
% Verify that the problem dimensions are defined properly.
%
if( minput<=0 )
xern1=sprintf(['%8i'], minput);
xermsg('SLATEC','DBOLSM',['THE NUMBER OF ROWS = ',[xern1,' MUST BE POSITIVE.']],31,1);
mode = -31;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
%
if( ncols<=0 )
xern1=sprintf(['%8i'], ncols);
xermsg('SLATEC','DBOLSM',['THE NUMBER OF COLUMNS = ',[xern1,' MUST BE POSITIVE.']],32,1);
mode = -32;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
%
if( mdw<minput )
xern1=sprintf(['%8i'], mdw);
xern2=sprintf(['%8i'], minput);
xermsg('SLATEC','DBOLSM',['THE ROW DIMENSION OF W() = ',[xern1,[' MUST BE .GE. THE NUMBER OF ROWS = ',xern2]]],33,1);
mode = -33;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
%
% Verify that bound information is correct.
%
for j = 1 : ncols;
if( ind(j)<1 || ind(j)>4 )
xern1=sprintf(['%8i'], j);
xern2=sprintf(['%8i'], ind(j));
xermsg('SLATEC','DBOLSM',['FOR J = ',[xern1,' THE CONSTRAINT INDICATOR MUST BE 1-4']],34,1);
mode = -34;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
end; j = fix(ncols+1);
%
for j = 1 : ncols;
if( ind(j)==3 )
if( bu(j)<bl(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','DBOLSM',['FOR J = ',[xern1,[' THE LOWER BOUND = ',[xern3,[' IS .GT. THE UPPER BOUND = ',xern4]]]]],35,1);
mode = -35;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
end;
end; j = fix(ncols+1);
%
% Check that permutation and polarity arrays have been set.
%
for j = 1 : ncols;
if( ibasis(j)<1 || ibasis(j)>ncols )
xern1=sprintf(['%8i'], ibasis(j));
xern2=sprintf(['%8i'], ncols);
xermsg('SLATEC','DBOLSM',['THE INPUT ORDER OF COLUMNS = ',[xern1,[' IS NOT BETWEEN 1 AND NCOLS = ',xern2]]],36,1);
mode = -36;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
%
if( ibb(j)<=0 )
xern1=sprintf(['%8i'], j);
xern2=sprintf(['%8i'], ibb(j));
xermsg('SLATEC','DBOLSM',['THE BOUND POLARITY FLAG IN COMPONENT J = ',[xern1,[' MUST BE POSITIVE.$$NOW = ',xern2]]],37,1);
mode = -37;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
end; j = fix(ncols+1);
%
% Process the option array.
%
fac = 0.75d0;
tolind = sqrt(d1mach(4));
tolsze = sqrt(d1mach(4));
itmax = fix(5.*max(minput,ncols));
wt = one;
mval = 0;
iprint = 0;
%
% Changes to some parameters can occur through the option array,
% IOPT(*). Process this array looking carefully for input data
% errors.
%
lp = 0;
lds = 0;
%
% Test for no more options.
%
gt(:)=0;
while (1);
if(gt(5)==0)
if(gt(4)==0)
if(gt(3)==0)
if(gt(2)==0)
if(gt(1)==0)
lp = fix(lp + lds);
ip = fix(iopt(lp+1));
jp = fix(abs(ip));
if( ip==99 )
%
% Pretriangularize rectangular arrays of certain sizes for
% increased efficiency.
%
if( fac.*minput>ncols )
for j = 1 : ncols + 1;
for i = minput : -1: j + mval + 1 ;
[w(i-1,j),w(i,j),sc,ss]=drotg(w(i-1,j),w(i,j),sc,ss);
w(i,j) = zero;
mdw_orig=mdw; [dumvar1,dumvar2,mdw,dumvar4,dumvar5,sc,ss]=drot(ncols-j+1,w(sub2ind(size(w),i-1,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),i-1,j+1):end)=dumvar2; w(sub2ind(size(w),i,j+1):end)=dumvar4;
end; i = fix(j + mval + 1 -1);
end; j = fix(ncols + 1+1);
mrows = fix(ncols + mval + 1);
else;
mrows = fix(minput);
end;
%
% Set the X(*) array to zero so all components are defined.
%
[ncols,zero,dumvar3,x]=dcopy(ncols,zero,0,x,1);
%
% The arrays IBASIS(*) and IBB(*) are initialized by the calling
% program and the column scaling is defined in the calling program.
% 'BIG' is plus infinity on this machine.
%
[big ]=d1mach(2);
for j = 1 : ncols;
if( ind(j)==1 )
bu(j) = big;
elseif( ind(j)==2 ) ;
bl(j) = -big;
elseif( ind(j)==4 ) ;
bl(j) = -big;
bu(j) = big;
end;
end; j = fix(ncols+1);
%
for j = 1 : ncols;
if((bl(j)<=zero && zero<=bu(j) && abs(bu(j))<abs(bl(j)))|| bu(j)<zero )
t = bu(j);
bu(j) = -bl(j);
bl(j) = -t;
scl(j) = -scl(j);
for i = 1 : mrows;
w(i,j) = -w(i,j);
end; i = fix(mrows+1);
end;
%
% Indices in set T(=TIGHT) are denoted by negative values
% of IBASIS(*).
%
if( bl(j)>=zero )
ibasis(j) = fix(-ibasis(j));
t = -bl(j);
bu(j) = bu(j) + t;
[mrows,t,dumvar3,dumvar4,dumvar5]=daxpy(mrows,t,w(sub2ind(size(w),1,j):end),1,w(sub2ind(size(w),1,ncols+1):end),1); w(sub2ind(size(w),1,j):end)=dumvar3; w(sub2ind(size(w),1,ncols+1):end)=dumvar5;
end;
end; j = fix(ncols+1);
%
nsetb = 0;
iter = 0;
%
if( iprint>0 )
[mrows,dumvar2,mdw,w]=dmout(mrows,ncols+1,mdw,w,'('' PRETRI. INPUT MATRIX'')',-4);
[ncols,bl]=dvout(ncols,bl,'('' LOWER BOUNDS'')',-4);
[ncols,bu]=dvout(ncols,bu,'('' UPPER BOUNDS'')',-4);
end;
gt(1)=1;
continue;
elseif( jp==99 ) ;
lds = 1;
elseif( jp==1 ) ;
%
% Move the IOPT(*) processing pointer.
%
if( ip>0 )
lp = fix(iopt(lp+2) - 1);
lds = 0;
else;
lds = 2;
end;
elseif( jp==2 ) ;
%
% Change tolerance for rank determination.
%
if( ip>0 )
ioff = fix(iopt(lp+2));
if( ioff<=0 )
xern1=sprintf(['%8i'], ioff);
xern2=sprintf(['%8i'], ncols);
xermsg('SLATEC','DBOLSM',['THE OFFSET = ',[xern1,[' BEYOND POSITION NCOLS = ',[xern2,' MUST BE POSITIVE FOR OPTION NUMBER 2.']]]],24,1);
mode = -24;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
%
tolind = x(ncols+ioff);
if( tolind<d1mach(4) )
xern3=sprintf([repmat('%15.6f',1,1)], tolind);
xern4=sprintf([repmat('%15.6f',1,1)], d1mach(4));
xermsg('SLATEC','DBOLSM',['THE TOLERANCE FOR RANK DETERMINATION = ',[xern3,[' IS LESS THAN MACHINE PRECISION = ',xern4]]],25,0);
mode = -25;
end;
end;
lds = 2;
elseif( jp==3 ) ;
%
% Change blowup factor for allowing variables to become
% inactive.
%
if( ip>0 )
ioff = fix(iopt(lp+2));
if( ioff<=0 )
xern1=sprintf(['%8i'], ioff);
xern2=sprintf(['%8i'], ncols);
xermsg('SLATEC','DBOLSM',['THE OFFSET = ',[xern1,[' BEYOND POSITION NCOLS = ',[xern2,' MUST BE POSITIVE FOR OPTION NUMBER 3.']]]],26,1);
mode = -26;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
%
tolsze = x(ncols+ioff);
if( tolsze<=zero )
xern3=sprintf([repmat('%15.6f',1,1)], tolsze);
xermsg('SLATEC','DBOLSM',['THE RECIPROCAL ',['OF THE BLOW-UP FACTOR FOR REJECTING VARIABLES ',['MUST BE POSITIVE.$$NOW = ',xern3]]],27,1);
mode = -27;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
end;
lds = 2;
elseif( jp==4 ) ;
%
% Change the maximum number of iterations allowed.
%
if( ip>0 )
itmax = fix(iopt(lp+2));
if( itmax<=0 )
xern1=sprintf(['%8i'], itmax);
xermsg('SLATEC','DBOLSM',['THE MAXIMUM NUMBER OF ITERATIONS = ',[xern1,' MUST BE POSITIVE.']],28,1);
mode = -28;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
end;
lds = 2;
elseif( jp==5 ) ;
%
% Change the factor for pretriangularizing the data matrix.
%
if( ip>0 )
ioff = fix(iopt(lp+2));
if( ioff<=0 )
xern1=sprintf(['%8i'], ioff);
xern2=sprintf(['%8i'], ncols);
xermsg('SLATEC','DBOLSM',['THE OFFSET = ',[xern1,[' BEYOND POSITION NCOLS = ',[xern2,' MUST BE POSITIVE FOR OPTION NUMBER 5.']]]],29,1);
mode = -29;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
%
fac = x(ncols+ioff);
if( fac<zero )
xern3=sprintf([repmat('%15.6f',1,1)], fac);
xermsg('SLATEC','DBOLSM',['THE FACTOR (NCOLS/MINPUT) WHERE PRE-',['TRIANGULARIZING IS PERFORMED MUST BE NON-',['NEGATIVE.$$NOW = ',xern3]]],30,0);
mode = -30;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
end;
lds = 2;
elseif( jp==6 ) ;
%
% Change the weighting factor (from 1.0) to apply to components
% numbered > MVAL (initially set to 1.) This trick is needed
% for applications of this subprogram to the heavily weighted
% least squares problem that come from equality constraints.
%
if( ip>0 )
ioff = fix(iopt(lp+2));
mval = fix(iopt(lp+3));
wt = x(ncols+ioff);
end;
%
if( mval<0 || mval>minput || wt<=zero )
xern1=sprintf(['%8i'], mval);
xern2=sprintf(['%8i'], minput);
xern3=sprintf([repmat('%15.6f',1,1)], wt);
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;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
lds = 3;
elseif( jp==7 ) ;
%
% Turn on debug output.
%
if( ip>0 )
iprint = 1;
end;
lds = 2;
else;
xern1=sprintf(['%8i'], ip);
xermsg('SLATEC','DBOLSM',['THE OPTION NUMBER = ',[xern1,' IS NOT DEFINED.']],23,1);
mode = -23;
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
return;
end;
gt(:)=0;
continue;
end;
gt(1)=0;
%
iter = fix(iter + 1);
if( iter>itmax )
xern1=sprintf(['%8i'], itmax);
xermsg('SLATEC','DBOLSM',['MORE THAN ITMAX = ',[xern1,' ITERATIONS SOLVING BOUNDED LEAST SQUARES PROBLEM.']],22,1);
mode = -22;
%
% Rescale and translate variables.
%
igopr = 1;
break;
end;
%
% Find a variable to become non-active.
% T
% Compute (negative) of gradient vector, W = E *(F-E*X).
%
[ncols,zero,dumvar3,ww]=dcopy(ncols,zero,0,ww,1);
for j = nsetb + 1 : ncols;
jcol = fix(abs(ibasis(j)));
inext= @(idum) min(idum+1,mrows);
ww(j) = ddot(mrows-nsetb,w(sub2ind(size(w),inext(nsetb),j):end),1,w(sub2ind(size(w),inext(nsetb),ncols+1):end),1).*abs(scl(jcol));
end; j = fix(ncols+1);
%
if( iprint>0 )
[ncols,ww]=dvout(ncols,ww,'('' GRADIENT VALUES'')',-4);
[ncols,ibasis]=ivout(ncols,ibasis,'('' INTERNAL VARIABLE ORDER'')',-4);
[ncols,ibb]=ivout(ncols,ibb,'('' BOUND POLARITY'')',-4);
end;
%
% If active set = number of total rows, quit.
%
end;
gt(2)=0;
if( nsetb==mrows )
found = false;
gt(5)=1;
continue;
end;
%
% Choose an extremal component of gradient vector for a candidate
% to become non-active.
%
wlarge = -big;
wmag = -big;
for j = nsetb + 1 : ncols;
t = ww(j);
if( t~=big )
itemp = fix(ibasis(j));
jcol = fix(abs(itemp));
inext= @(idum) min(idum+1,mrows);
[t1 ,dumvar2,w(sub2ind(size(w),inext(nsetb),j):end)]=dnrm2(mval-nsetb,w(sub2ind(size(w),inext(nsetb),j):end),1);
if( itemp<0 )
if( rem(ibb(jcol),2)==0 )
t = -t;
end;
if( t>=zero )
if( mval>nsetb )
t = t1;
end;
if( t>wlarge )
wlarge = t;
jlarge = fix(j);
end;
end;
else;
if( mval>nsetb )
t = t1;
end;
if( abs(t)>wmag )
wmag = abs(t);
jmag = fix(j);
end;
end;
end;
end; j = fix(ncols+1);
%
% Choose magnitude of largest component of gradient for candidate.
%
jbig = 0;
wbig = zero;
if( wlarge>zero )
jbig = fix(jlarge);
wbig = wlarge;
end;
%
if( wmag>=wbig )
jbig = fix(jmag);
wbig = wmag;
end;
%
if( jbig==0 )
found = false;
if( iprint>0 )
[dumvar1,i]=ivout(0,i,'('' FOUND NO VARIABLE TO ENTER'')',-4);
end;
gt(5)=1;
continue;
end;
%
% See if the incoming column is sufficiently independent. This
% test is made before an elimination is performed.
%
if( iprint>0 )
[dumvar1,jbig]=ivout(1,jbig,'('' TRY TO BRING IN THIS COL.'')',-4);
end;
%
if( mval<=nsetb )
[cl1 ,mval,w(sub2ind(size(w),1,jbig):end)]=dnrm2(mval,w(sub2ind(size(w),1,jbig):end),1);
inext= @(idum) min(idum+1,mrows);
cl2 = abs(wt).*dnrm2(nsetb-mval,w(sub2ind(size(w),inext(mval),jbig):end),1);
inext= @(idum) min(idum+1,mrows);
cl3 = abs(wt).*dnrm2(mrows-nsetb,w(sub2ind(size(w),inext(nsetb),jbig):end),1);
[cl1,cl2,sc,ss]=drotg(cl1,cl2,sc,ss);
colabv = abs(cl1);
colblo = cl3;
else;
[cl1 ,nsetb,w(sub2ind(size(w),1,jbig):end)]=dnrm2(nsetb,w(sub2ind(size(w),1,jbig):end),1);
inext= @(idum) min(idum+1,mrows);
[cl2 ,dumvar2,w(sub2ind(size(w),inext(nsetb),jbig):end)]=dnrm2(mval-nsetb,w(sub2ind(size(w),inext(nsetb),jbig):end),1);
inext= @(idum) min(idum+1,mrows);
cl3 = abs(wt).*dnrm2(mrows-mval,w(sub2ind(size(w),inext(mval),jbig):end),1);
colabv = cl1;
[cl2,cl3,sc,ss]=drotg(cl2,cl3,sc,ss);
colblo = abs(cl2);
end;
%
if( colblo<=tolind.*colabv )
ww(jbig) = big;
if( iprint>0 )
[dumvar1,i]=ivout(0,i,'('' VARIABLE IS DEPENDENT, NOT USED.'')',-4);
end;
gt(2)=1;
continue;
end;
%
% Swap matrix columns NSETB+1 and JBIG, plus pointer information,
% and gradient values.
%
nsetb = fix(nsetb + 1);
if( nsetb~=jbig )
[mrows,dumvar2,dumvar3,dumvar4]=dswap(mrows,w(sub2ind(size(w),1,nsetb):end),1,w(sub2ind(size(w),1,jbig):end),1); w(sub2ind(size(w),1,nsetb):end)=dumvar2; w(sub2ind(size(w),1,jbig):end)=dumvar4;
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,ww(sub2ind(size(ww),max(nsetb,1)):end),1,ww(sub2ind(size(ww),max(jbig,1)):end),1); dumvar2i=find((ww(sub2ind(size(ww),max(nsetb,1)):end))~=(dumvar2));dumvar4i=find((ww(sub2ind(size(ww),max(jbig,1)):end))~=(dumvar4)); ww(nsetb-1+dumvar2i)=dumvar2(dumvar2i); ww(jbig-1+dumvar4i)=dumvar4(dumvar4i);
itemp = fix(ibasis(nsetb));
ibasis(nsetb) = fix(ibasis(jbig));
ibasis(jbig) = fix(itemp);
end;
%
% Eliminate entries below the pivot line in column NSETB.
%
if( mrows>nsetb )
for i = mrows : -1: nsetb + 1 ;
if( i~=mval+1 )
[w(i-1,nsetb),w(i,nsetb),sc,ss]=drotg(w(i-1,nsetb),w(i,nsetb),sc,ss);
w(i,nsetb) = zero;
mdw_orig=mdw; [dumvar1,dumvar2,mdw,dumvar4,dumvar5,sc,ss]=drot(ncols-nsetb+1,w(sub2ind(size(w),i-1,nsetb+1):end),mdw,w(sub2ind(size(w),i,nsetb+1):end),mdw,sc,ss); mdw(dumvar5~=mdw_orig)=dumvar5(dumvar5~=mdw_orig); w(sub2ind(size(w),i-1,nsetb+1):end)=dumvar2; w(sub2ind(size(w),i,nsetb+1):end)=dumvar4;
end;
end; i = fix(nsetb + 1 -1);
%
if( mval>=nsetb && mval<mrows )
[w(nsetb,nsetb),w(mval+1,nsetb),sc,ss]=drotg(w(nsetb,nsetb),w(mval+1,nsetb),sc,ss);
w(mval+1,nsetb) = zero;
mdw_orig=mdw; [dumvar1,dumvar2,mdw,dumvar4,dumvar5,sc,ss]=drot(ncols-nsetb+1,w(sub2ind(size(w),nsetb,nsetb+1):end),mdw,w(sub2ind(size(w),mval+1,nsetb+1):end),mdw,sc,ss); mdw(dumvar5~=mdw_orig)=dumvar5(dumvar5~=mdw_orig); w(sub2ind(size(w),nsetb,nsetb+1):end)=dumvar2; w(sub2ind(size(w),mval+1,nsetb+1):end)=dumvar4;
end;
end;
%
if( w(nsetb,nsetb)==zero )
ww(nsetb) = big;
nsetb = fix(nsetb - 1);
if( iprint>0 )
[dumvar1,i]=ivout(0,i,'('' PIVOT IS ZERO, NOT USED.'')',-4);
end;
gt(2)=1;
continue;
end;
%
% Check that new variable is moving in the right direction.
%
itemp = fix(ibasis(nsetb));
jcol = fix(abs(itemp));
xnew =(w(nsetb,ncols+1)./w(nsetb,nsetb))./abs(scl(jcol));
if( itemp<0 )
%
% IF(WW(NSETB).GE.ZERO.AND.XNEW.LE.ZERO) exit(quit)
% IF(WW(NSETB).LE.ZERO.AND.XNEW.GE.ZERO) exit(quit)
%
if((ww(nsetb)>=zero && xnew<=zero) ||(ww(nsetb)<=zero && xnew>=zero) )
%
ww(nsetb) = big;
nsetb = fix(nsetb - 1);
if( iprint>0 )
[dumvar1,i]=ivout(0,i,'('' VARIABLE HAS BAD DIRECTION, NOT USED.'')',-4);
end;
gt(2)=1;
continue;
end;
end;
found = true;
gt(5)=1;
continue;
%
% Solve the triangular system.
%
end;
gt(3)=0;
[nsetb,w(sub2ind(size(w),1,ncols+1):end),dumvar3,rw]=dcopy(nsetb,w(sub2ind(size(w),1,ncols+1):end),1,rw,1);
for j = nsetb : -1: 1 ;
rw(j) = rw(j)./w(j,j);
jcol = fix(abs(ibasis(j)));
t = rw(j);
if( rem(ibb(jcol),2)==0 )
rw(j) = -rw(j);
end;
[dumvar1,dumvar2,w(sub2ind(size(w),1,j):end),dumvar4,rw]=daxpy(j-1,-t,w(sub2ind(size(w),1,j):end),1,rw,1);
rw(j) = rw(j)./abs(scl(jcol));
end; j = fix(1 -1);
%
if( iprint>0 )
[nsetb,rw]=dvout(nsetb,rw,'('' SOLN. VALUES'')',-4);
[nsetb,ibasis]=ivout(nsetb,ibasis,'('' COLS. USED'')',-4);
end;
%
if( lgopr==2 )
[nsetb,rw,dumvar3,x]=dcopy(nsetb,rw,1,x,1);
for j = 1 : nsetb;
itemp = fix(ibasis(j));
jcol = fix(abs(itemp));
if( itemp<0 )
bou = zero;
else;
bou = bl(jcol);
end;
%
if((-bou)~=big )
bou = bou./abs(scl(jcol));
end;
if( x(j)<=bou )
jdrop1 = fix(j);
break;
end;
%
bou = bu(jcol);
if( bou~=big )
bou = bou./abs(scl(jcol));
end;
if( x(j)>=bou )
jdrop2 = fix(j);
break;
end;
end;
gt(4)=1;
continue;
end;
%
% See if the unconstrained solution (obtained by solving the
% triangular system) satisfies the problem bounds.
%
alpha = two;
beta = two;
x(nsetb) = zero;
for j = 1 : nsetb;
itemp = fix(ibasis(j));
jcol = fix(abs(itemp));
t1 = two;
t2 = two;
if( itemp<0 )
bou = zero;
else;
bou = bl(jcol);
end;
if((-bou)~=big )
bou = bou./abs(scl(jcol));
end;
if( rw(j)<=bou )
t1 =(x(j)-bou)./(x(j)-rw(j));
end;
bou = bu(jcol);
if( bou~=big )
bou = bou./abs(scl(jcol));
end;
if( rw(j)>=bou )
t2 =(bou-x(j))./(rw(j)-x(j));
end;
%
% If not, then compute a step length so that the variables remain
% feasible.
%
if( t1<alpha )
alpha = t1;
jdrop1 = fix(j);
end;
%
if( t2<beta )
beta = t2;
jdrop2 = fix(j);
end;
end; j = fix(nsetb+1);
%
constr = alpha<two | beta<two;
if( ~constr )
%
% Accept the candidate because it satisfies the stated bounds
% on the variables.
%
[nsetb,rw,dumvar3,x]=dcopy(nsetb,rw,1,x,1);
gt(1)=1;
continue;
end;
%
% Take a step that is as large as possible with all variables
% remaining feasible.
%
for j = 1 : nsetb;
x(j) = x(j) + min(alpha,beta).*(rw(j)-x(j));
end; j = fix(nsetb+1);
%
if( alpha<=beta )
jdrop2 = 0;
else;
jdrop1 = 0;
end;
%
end;
gt(4)=0;
if( jdrop1+jdrop2<=0 || nsetb<=0 )
gt(1)=1;
continue;
end;
jdrop = fix(jdrop1 + jdrop2);
itemp = fix(ibasis(jdrop));
jcol = fix(abs(itemp));
if( jdrop2>0 )
%
% Variable is at an upper bound. Subtract multiple of this
% column from right hand side.
%
t = bu(jcol);
if( itemp>0 )
bu(jcol) = t - bl(jcol);
bl(jcol) = -t;
itemp = fix(-itemp);
scl(jcol) = -scl(jcol);
for i = 1 : jdrop;
w(i,jdrop) = -w(i,jdrop);
end; i = fix(jdrop+1);
else;
ibb(jcol) = fix(ibb(jcol) + 1);
if( rem(ibb(jcol),2)==0 )
t = -t;
end;
end;
%
% Variable is at a lower bound.
%
elseif( itemp<zero ) ;
t = zero;
else;
t = -bl(jcol);
bu(jcol) = bu(jcol) + t;
itemp = fix(-itemp);
end;
%
[jdrop,t,dumvar3,dumvar4,dumvar5]=daxpy(jdrop,t,w(sub2ind(size(w),1,jdrop):end),1,w(sub2ind(size(w),1,ncols+1):end),1); w(sub2ind(size(w),1,jdrop):end)=dumvar3; w(sub2ind(size(w),1,ncols+1):end)=dumvar5;
%
% Move certain columns left to achieve upper Hessenberg form.
%
[jdrop,w(sub2ind(size(w),1,jdrop):end),dumvar3,rw]=dcopy(jdrop,w(sub2ind(size(w),1,jdrop):end),1,rw,1);
for j = jdrop + 1 : nsetb;
ibasis(j-1) = fix(ibasis(j));
x(j-1) = x(j);
[j,dumvar2,dumvar3,dumvar4]=dcopy(j,w(sub2ind(size(w),1,j):end),1,w(sub2ind(size(w),1,j-1):end),1); w(sub2ind(size(w),1,j):end)=dumvar2; w(sub2ind(size(w),1,j-1):end)=dumvar4;
end; j = fix(nsetb+1);
%
ibasis(nsetb) = fix(itemp);
w(1,nsetb) = zero;
[dumvar1,dumvar2,dumvar3,dumvar4]=dcopy(mrows-jdrop,w(sub2ind(size(w),1,nsetb):end),0,w(sub2ind(size(w),jdrop+1,nsetb):end),1); w(sub2ind(size(w),1,nsetb):end)=dumvar2; w(sub2ind(size(w),jdrop+1,nsetb):end)=dumvar4;
[jdrop,rw,dumvar3,w(sub2ind(size(w),1,nsetb):end)]=dcopy(jdrop,rw,1,w(sub2ind(size(w),1,nsetb):end),1);
%
% Transform the matrix from upper Hessenberg form to upper
% triangular form.
%
nsetb = fix(nsetb - 1);
for i = jdrop : nsetb;
%
% Look for small pivots and avoid mixing weighted and
% nonweighted rows.
%
if( i==mval )
t = zero;
for j = i : nsetb;
jcol = fix(abs(ibasis(j)));
t1 = abs(w(i,j).*scl(jcol));
if( t1>t )
jbig = fix(j);
t = t1;
end;
end; j = fix(nsetb+1);
%
% The triangularization is completed by giving up the Hessenberg
% form and triangularizing a rectangular matrix.
%
[mrows,dumvar2,dumvar3,dumvar4]=dswap(mrows,w(sub2ind(size(w),1,i):end),1,w(sub2ind(size(w),1,jbig):end),1); w(sub2ind(size(w),1,i):end)=dumvar2; w(sub2ind(size(w),1,jbig):end)=dumvar4;
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,ww(sub2ind(size(ww),max(i,1)):end),1,ww(sub2ind(size(ww),max(jbig,1)):end),1); dumvar2i=find((ww(sub2ind(size(ww),max(i,1)):end))~=(dumvar2));dumvar4i=find((ww(sub2ind(size(ww),max(jbig,1)):end))~=(dumvar4)); ww(i-1+dumvar2i)=dumvar2(dumvar2i); ww(jbig-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,x(sub2ind(size(x),max(i,1)):end),1,x(sub2ind(size(x),max(jbig,1)):end),1); dumvar2i=find((x(sub2ind(size(x),max(i,1)):end))~=(dumvar2));dumvar4i=find((x(sub2ind(size(x),max(jbig,1)):end))~=(dumvar4)); x(i-1+dumvar2i)=dumvar2(dumvar2i); x(jbig-1+dumvar4i)=dumvar4(dumvar4i);
itemp = fix(ibasis(i));
ibasis(i) = fix(ibasis(jbig));
ibasis(jbig) = fix(itemp);
jbig = fix(i);
for j = jbig : nsetb;
for ii = j + 1 : mrows;
[w(j,j),w(ii,j),sc,ss]=drotg(w(j,j),w(ii,j),sc,ss);
w(ii,j) = zero;
mdw_orig=mdw; [dumvar1,dumvar2,mdw,dumvar4,dumvar5,sc,ss]=drot(ncols-j+1,w(sub2ind(size(w),j,j+1):end),mdw,w(sub2ind(size(w),ii,j+1):end),mdw,sc,ss); mdw(dumvar5~=mdw_orig)=dumvar5(dumvar5~=mdw_orig); w(sub2ind(size(w),j,j+1):end)=dumvar2; w(sub2ind(size(w),ii,j+1):end)=dumvar4;
end; ii = fix(mrows+1);
end; j = fix(nsetb+1);
break;
end;
[w(i,i),w(i+1,i),sc,ss]=drotg(w(i,i),w(i+1,i),sc,ss);
w(i+1,i) = zero;
mdw_orig=mdw; [dumvar1,dumvar2,mdw,dumvar4,dumvar5,sc,ss]=drot(ncols-i+1,w(sub2ind(size(w),i,i+1):end),mdw,w(sub2ind(size(w),i+1,i+1):end),mdw,sc,ss); mdw(dumvar5~=mdw_orig)=dumvar5(dumvar5~=mdw_orig); w(sub2ind(size(w),i,i+1):end)=dumvar2; w(sub2ind(size(w),i+1,i+1):end)=dumvar4;
end;
%
% See if the remaining coefficients are feasible. They should be
% because of the way MIN(ALPHA,BETA) was chosen. Any that are not
% feasible will be set to their bounds and appropriately translated.
%
jdrop1 = 0;
jdrop2 = 0;
lgopr = 2;
gt(3)=1;
continue;
%
% Find a variable to become non-active.
%
end;
gt(5)=0;
if( found )
lgopr = 1;
gt(3)=1;
continue;
end;
%
% Rescale and translate variables.
%
igopr = 2;
break;
end;
[nsetb,x,dumvar3,rw]=dcopy(nsetb,x,1,rw,1);
[ncols,zero,dumvar3,x]=dcopy(ncols,zero,0,x,1);
for j = 1 : nsetb;
jcol = fix(abs(ibasis(j)));
x(jcol) = rw(j).*abs(scl(jcol));
end; j = fix(nsetb+1);
%
for j = 1 : ncols;
if( rem(ibb(j),2)==0 )
x(j) = bu(j) - x(j);
end;
end; j = fix(ncols+1);
%
for j = 1 : ncols;
jcol = fix(ibasis(j));
if( jcol<0 )
x(-jcol) = bl(-jcol) + x(-jcol);
end;
end; j = fix(ncols+1);
%
for j = 1 : ncols;
if( scl(j)<zero )
x(j) = -x(j);
end;
end; j = fix(ncols+1);
%
i = fix(max(nsetb,mval));
inext= @(idum) min(idum+1,mrows);
[rnorm ,dumvar2,w(sub2ind(size(w),inext(i),ncols+1):end)]=dnrm2(mrows-i,w(sub2ind(size(w),inext(i),ncols+1):end),1);
%
if( igopr==2 )
mode = fix(nsetb);
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;
ww_shape=zeros(ww_shape);ww_shape(:)=ww(1:numel(ww_shape));ww=ww_shape;
scl_shape=zeros(scl_shape);scl_shape(:)=scl(1:numel(scl_shape));scl=scl_shape;
ibasis_shape=zeros(ibasis_shape);ibasis_shape(:)=ibasis(1:numel(ibasis_shape));ibasis=ibasis_shape;
ibb_shape=zeros(ibb_shape);ibb_shape(:)=ibb(1:numel(ibb_shape));ibb=ibb_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;
end %subroutine dbolsm
%DECK DBSGQ8
|
|
