| [w,mdw,me,ma,mg,n,prgopt,x,rnorme,rnorml,mode,ws,ip]=lsei(w,mdw,me,ma,mg,n,prgopt,x,rnorme,rnorml,mode,ws,ip); |
function [w,mdw,me,ma,mg,n,prgopt,x,rnorme,rnorml,mode,ws,ip]=lsei(w,mdw,me,ma,mg,n,prgopt,x,rnorme,rnorml,mode,ws,ip);
%***BEGIN PROLOGUE LSEI
%***PURPOSE Solve a linearly constrained least squares problem with
% equality and inequality constraints, and optionally compute
% a covariance matrix.
%***LIBRARY SLATEC
%***CATEGORY K1A2A, D9
%***TYPE SINGLE PRECISION (LSEI-S, DLSEI-D)
%***KEYWORDS CONSTRAINED LEAST SQUARES, CURVE FITTING, DATA FITTING,
% EQUALITY CONSTRAINTS, INEQUALITY CONSTRAINTS,
% QUADRATIC PROGRAMMING
%***AUTHOR Hanson, R. J., (SNLA)
% Haskell, K. H., (SNLA)
%***DESCRIPTION
%
% Abstract
%
% This subprogram solves a linearly constrained least squares
% problem with both equality and inequality constraints, and, if the
% user requests, obtains a covariance matrix of the solution
% parameters.
%
% Suppose there are given matrices E, A and G of respective
% dimensions ME by N, MA by N and MG by N, and vectors F, B and H of
% respective lengths ME, MA and MG. This subroutine solves the
% linearly constrained least squares problem
%
% EX = F, (E ME by N) (equations to be exactly
% satisfied)
% AX = B, (A MA by N) (equations to be
% approximately satisfied,
% least squares sense)
% GX .GE. H,(G MG by N) (inequality constraints)
%
% The inequalities GX .GE. H mean that every component of the
% product GX must be .GE. the corresponding component of H.
%
% In case the equality constraints cannot be satisfied, a
% generalized inverse solution residual vector length is obtained
% for F-EX. This is the minimal length possible for F-EX.
%
% Any values ME .GE. 0, MA .GE. 0, or MG .GE. 0 are permitted. The
% rank of the matrix E is estimated during the computation. We call
% this value KRANKE. It is an output parameter in IP(1) defined
% below. Using a generalized inverse solution of EX=F, a reduced
% least squares problem with inequality constraints is obtained.
% The tolerances used in these tests for determining the rank
% of E and the rank of the reduced least squares problem are
% given in Sandia Tech. Rept. SAND-78-1290. They can be
% modified by the user if new values are provided in
% the option list of the array PRGOPT(*).
%
% The user must dimension all arrays appearing in the call list..
% W(MDW,N+1),PRGOPT(*),X(N),WS(2*(ME+N)+K+(MG+2)*(N+7)),IP(MG+2*N+2)
% where K=MAX(MA+MG,N). This allows for a solution of a range of
% problems in the given working space. The dimension of WS(*)
% given is a necessary overestimate. Once a particular problem
% has been run, the output parameter IP(3) gives the actual
% dimension required for that problem.
%
% The parameters for LSEI( ) are
%
% Input..
%
% W(*,*),MDW, The array W(*,*) is doubly subscripted with
% ME,MA,MG,N first dimensioning parameter equal to MDW.
% For this discussion let us call M = ME+MA+MG. Then
% MDW must satisfy MDW .GE. M. The condition
% MDW .LT. M is an error.
%
% The array W(*,*) contains the matrices and vectors
%
% (E F)
% (A B)
% (G H)
%
% in rows and columns 1,...,M and 1,...,N+1
% respectively.
%
% The integers ME, MA, and MG are the
% respective matrix row dimensions
% of E, A and G. Each matrix has N columns.
%
% PRGOPT(*) This real-valued array is the option vector.
% If the user is satisfied with the nominal
% subprogram features set
%
% PRGOPT(1)=1 (or PRGOPT(1)=1.0)
%
% Otherwise PRGOPT(*) is a linked list consisting of
% groups of data of the following form
%
% LINK
% KEY
% DATA SET
%
% The parameters LINK and KEY are each one word.
% The DATA SET can be comprised of several words.
% The number of items depends on the value of KEY.
% The value of LINK points to the first
% entry of the next group of data within
% PRGOPT(*). The exception is when there are
% no more options to change. In that
% case, LINK=1 and the values KEY and DATA SET
% are not referenced. The general layout of
% PRGOPT(*) is as follows.
%
% ...PRGOPT(1) = LINK1 (link to first entry of next group)
% . PRGOPT(2) = KEY1 (key to the option change)
% . PRGOPT(3) = data value (data value for this change)
% . .
% . .
% . .
% ...PRGOPT(LINK1) = LINK2 (link to the first entry of
% . next group)
% . PRGOPT(LINK1+1) = KEY2 (key to the option change)
% . PRGOPT(LINK1+2) = data value
% ... .
% . .
% . .
% ...PRGOPT(LINK) = 1 (no more options to change)
%
% Values of LINK that are nonpositive are errors.
% A value of LINK .GT. NLINK=100000 is also an error.
% This helps prevent using invalid but positive
% values of LINK that will probably extend
% beyond the program limits of PRGOPT(*).
% Unrecognized values of KEY are ignored. The
% order of the options is arbitrary and any number
% of options can be changed with the following
% restriction. To prevent cycling in the
% processing of the option array, a count of the
% number of options changed is maintained.
% Whenever this count exceeds NOPT=1000, an error
% message is printed and the subprogram returns.
%
% Options..
%
% KEY=1
% Compute in W(*,*) the N by N
% covariance matrix of the solution variables
% as an output parameter. Nominally the
% covariance matrix will not be computed.
% (This requires no user input.)
% The data set for this option is a single value.
% It must be nonzero when the covariance matrix
% is desired. If it is zero, the covariance
% matrix is not computed. When the covariance matrix
% is computed, the first dimensioning parameter
% of the array W(*,*) must satisfy MDW .GE. MAX(M,N).
%
% KEY=10
% Suppress scaling of the inverse of the
% normal matrix by the scale factor RNORM**2/
% MAX(1, no. of degrees of freedom). This option
% only applies when the option for computing the
% covariance matrix (KEY=1) is used. With KEY=1 and
% KEY=10 used as options the unscaled inverse of the
% normal matrix is returned in W(*,*).
% The data set for this option is a single value.
% When it is nonzero no scaling is done. When it is
% zero scaling is done. The nominal case is to do
% scaling so if option (KEY=1) is used alone, the
% matrix will be scaled on output.
%
% KEY=2
% Scale the nonzero columns of the
% entire data matrix.
% (E)
% (A)
% (G)
%
% to have length one. The data set for this
% option is a single value. It must be
% nonzero if unit length column scaling
% is desired.
%
% KEY=3
% Scale columns of the entire data matrix
% (E)
% (A)
% (G)
%
% with a user-provided diagonal matrix.
% The data set for this option consists
% of the N diagonal scaling factors, one for
% each matrix column.
%
% KEY=4
% Change the rank determination tolerance for
% the equality constraint equations from
% the nominal value of SQRT(SRELPR). This quantity can
% be no smaller than SRELPR, the arithmetic-
% storage precision. The quantity SRELPR is the
% largest positive number such that T=1.+SRELPR
% satisfies T .EQ. 1. The quantity used
% here is internally restricted to be at
% least SRELPR. The data set for this option
% is the new tolerance.
%
% KEY=5
% Change the rank determination tolerance for
% the reduced least squares equations from
% the nominal value of SQRT(SRELPR). This quantity can
% be no smaller than SRELPR, the arithmetic-
% storage precision. The quantity used
% here is internally restricted to be at
% least SRELPR. The data set for this option
% is the new tolerance.
%
% For example, suppose we want to change
% the tolerance for the reduced least squares
% problem, compute the covariance matrix of
% the solution parameters, and provide
% column scaling for the data matrix. For
% these options the dimension of PRGOPT(*)
% must be at least N+9. The Fortran statements
% defining these options would be as follows:
%
% PRGOPT(1)=4 (link to entry 4 in PRGOPT(*))
% PRGOPT(2)=1 (covariance matrix key)
% PRGOPT(3)=1 (covariance matrix wanted)
%
% PRGOPT(4)=7 (link to entry 7 in PRGOPT(*))
% PRGOPT(5)=5 (least squares equas. tolerance key)
% PRGOPT(6)=... (new value of the tolerance)
%
% PRGOPT(7)=N+9 (link to entry N+9 in PRGOPT(*))
% PRGOPT(8)=3 (user-provided column scaling key)
%
% CALL SCOPY (N, D, 1, PRGOPT(9), 1) (Copy the N
% scaling factors from the user array D(*)
% to PRGOPT(9)-PRGOPT(N+8))
%
% PRGOPT(N+9)=1 (no more options to change)
%
% The contents of PRGOPT(*) are not modified
% by the subprogram.
% The options for WNNLS( ) can also be included
% in this array. The values of KEY recognized
% by WNNLS( ) are 6, 7 and 8. Their functions
% are documented in the usage instructions for
% subroutine WNNLS( ). Normally these options
% do not need to be modified when using LSEI( ).
%
% IP(1), The amounts of working storage actually
% IP(2) allocated for the working arrays WS(*) and
% IP(*), respectively. These quantities are
% compared with the actual amounts of storage
% needed by LSEI( ). Insufficient storage
% allocated for either WS(*) or IP(*) is an
% error. This feature was included in LSEI( )
% because miscalculating the storage formulas
% for WS(*) and IP(*) might very well lead to
% subtle and hard-to-find execution errors.
%
% The length of WS(*) must be at least
%
% LW = 2*(ME+N)+K+(MG+2)*(N+7)
%
% where K = max(MA+MG,N)
% This test will not be made if IP(1).LE.0.
%
% The length of IP(*) must be at least
%
% LIP = MG+2*N+2
% This test will not be made if IP(2).LE.0.
%
% Output..
%
% X(*),RNORME, The array X(*) contains the solution parameters
% RNORML if the integer output flag MODE = 0 or 1.
% The definition of MODE is given directly below.
% When MODE = 0 or 1, RNORME and RNORML
% respectively contain the residual vector
% Euclidean lengths of F - EX and B - AX. When
% MODE=1 the equality constraint equations EX=F
% are contradictory, so RNORME .NE. 0. The residual
% vector F-EX has minimal Euclidean length. For
% MODE .GE. 2, none of these parameters is defined.
%
% MODE Integer flag that indicates the subprogram
% status after completion. If MODE .GE. 2, no
% solution has been computed.
%
% MODE =
%
% 0 Both equality and inequality constraints
% are compatible and have been satisfied.
%
% 1 Equality constraints are contradictory.
% A generalized inverse solution of EX=F was used
% to minimize the residual vector length F-EX.
% In this sense, the solution is still meaningful.
%
% 2 Inequality constraints are contradictory.
%
% 3 Both equality and inequality constraints
% are contradictory.
%
% The following interpretation of
% MODE=1,2 or 3 must be made. The
% sets consisting of all solutions
% of the equality constraints EX=F
% and all vectors satisfying GX .GE. H
% have no points in common. (In
% particular this does not say that
% each individual set has no points
% at all, although this could be the
% case.)
%
% 4 Usage error occurred. The value
% of MDW is .LT. ME+MA+MG, MDW is
% .LT. N and a covariance matrix is
% requested, or the option vector
% PRGOPT(*) is not properly defined,
% or the lengths of the working arrays
% WS(*) and IP(*), when specified in
% IP(1) and IP(2) respectively, are not
% long enough.
%
% W(*,*) The array W(*,*) contains the N by N symmetric
% covariance matrix of the solution parameters,
% provided this was requested on input with
% the option vector PRGOPT(*) and the output
% flag is returned with MODE = 0 or 1.
%
% IP(*) The integer working array has three entries
% that provide rank and working array length
% information after completion.
%
% IP(1) = rank of equality constraint
% matrix. Define this quantity
% as KRANKE.
%
% IP(2) = rank of reduced least squares
% problem.
%
% IP(3) = the amount of storage in the
% working array WS(*) that was
% actually used by the subprogram.
% The formula given above for the length
% of WS(*) is a necessary overestimate.
% If exactly the same problem matrices
% are used in subsequent executions,
% the declared dimension of WS(*) can
% be reduced to this output value.
% User Designated
% Working Arrays..
%
% WS(*),IP(*) These are respectively type real
% and type integer working arrays.
% Their required minimal lengths are
% given above.
%
%***REFERENCES K. H. Haskell and R. J. Hanson, An algorithm for
% linear least squares problems with equality and
% nonnegativity constraints, Report SAND77-0552, Sandia
% Laboratories, June 1978.
% K. H. Haskell and R. J. Hanson, Selected algorithms for
% the linearly constrained least squares problem - a
% users guide, Report SAND78-1290, Sandia Laboratories,
% August 1979.
% K. H. Haskell and R. J. Hanson, An algorithm for
% linear least squares problems with equality and
% nonnegativity constraints, Mathematical Programming
% 21 (1981), pp. 98-118.
% R. J. Hanson and K. H. Haskell, Two algorithms for the
% linearly constrained least squares problem, ACM
% Transactions on Mathematical Software, September 1982.
%***ROUTINES CALLED H12, LSI, R1MACH, SASUM, SAXPY, SCOPY, SDOT, SNRM2,
% SSCAL, SSWAP, XERMSG
%***REVISION HISTORY (YYMMDD)
% 790701 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890618 Completely restructured and extensively revised (WRBRWC)
% 890831 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
% 900510 Convert XERRWV calls to XERMSG calls. (RWC)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE LSEI
persistent cov enorm first firstCall fnorm gam gt i imax j jp1 k key kranke last lchk link m mapke1 mdeqc mend mep1 n1 n2 next nlink nopt np1 ntimes rb rn rnmax sizemlv sn snmax srelpr t tau uj up vj xern1 xern2 xern3 xern4 xnorm xnrme ; if isempty(firstCall),firstCall=1;end;
if isempty(gt), gt=0; end;
prgopt_shape=size(prgopt);prgopt=reshape(prgopt,1,[]);
w_shape=size(w);w=reshape([w(:).',zeros(1,ceil(numel(w)./prod([mdw])).*prod([mdw])-numel(w))],mdw,[]);
ws_shape=size(ws);ws=reshape(ws,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
%
%
if isempty(enorm), enorm=0; end;
if isempty(fnorm), fnorm=0; end;
if isempty(gam), gam=0; end;
if isempty(rb), rb=0; end;
if isempty(rn), rn=0; end;
if isempty(rnmax), rnmax=0; end;
if isempty(sizemlv), sizemlv=0; end;
if isempty(sn), sn=0; end;
if isempty(snmax), snmax=0; end;
if isempty(srelpr), srelpr=0; end;
if isempty(t), t=0; end;
if isempty(tau), tau=0; end;
if isempty(uj), uj=0; end;
if isempty(up), up=0; end;
if isempty(vj), vj=0; end;
if isempty(xnorm), xnorm=0; end;
if isempty(xnrme), xnrme=0; end;
if isempty(i), i=0; end;
if isempty(imax), imax=0; end;
if isempty(j), j=0; end;
if isempty(jp1), jp1=0; end;
if isempty(k), k=0; end;
if isempty(key), key=0; end;
if isempty(kranke), kranke=0; end;
if isempty(last), last=0; end;
if isempty(lchk), lchk=0; end;
if isempty(link), link=0; end;
if isempty(m), m=0; end;
if isempty(mapke1), mapke1=0; end;
if isempty(mdeqc), mdeqc=0; end;
if isempty(mend), mend=0; end;
if isempty(mep1), mep1=0; end;
if isempty(n1), n1=0; end;
if isempty(n2), n2=0; end;
if isempty(next), next=0; end;
if isempty(nlink), nlink=0; end;
if isempty(nopt), nopt=0; end;
if isempty(np1), np1=0; end;
if isempty(ntimes), ntimes=0; end;
if isempty(cov), cov=false; end;
if isempty(first), first=false; end;
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern2), xern2=repmat(' ',1,8); end;
if isempty(xern3), xern3=repmat(' ',1,8); end;
if isempty(xern4), xern4=repmat(' ',1,8); end;
%
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT LSEI
%
% Set the nominal tolerance used in the code for the equality
% constraint equations.
%
if( first )
[ srelpr ]=r1mach(4);
end;
first = false;
tau = sqrt(srelpr);
%
% Check that enough storage was allocated in WS(*) and IP(*).
%
mode = 4;
if( min([n,me,ma,mg])<0 )
xern1=sprintf(['%8i'], n);
xern2=sprintf(['%8i'], me);
xern3=sprintf(['%8i'], ma);
xern4=sprintf(['%8i'], mg);
xermsg('SLATEC','LSEI',['ALL OF THE VARIABLES N, ME,',[' MA, MG MUST BE .GE. 0$$ENTERED ROUTINE WITH',['$$N = ',[xern1,['$$ME = ',[xern2,['$$MA = ',[xern3,['$$MG = ',xern4]]]]]]]]],2,1);
else;
%
if( ip(1)>0 )
lchk = fix(2.*(me+n) + max(ma+mg,n) +(mg+2).*(n+7));
if( ip(1)<lchk )
xern1=sprintf(['%8i'], lchk);
xermsg('SLATEC','LSEI',['INSUFFICIENT STORAGE ',['ALLOCATED FOR WS(*), NEED LW = ',xern1]],2,1);
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
ws_shape=zeros(ws_shape);ws_shape(:)=ws(1:numel(ws_shape));ws=ws_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
return;
end;
end;
%
if( ip(2)>0 )
lchk = fix(mg + 2.*n + 2);
if( ip(2)<lchk )
xern1=sprintf(['%8i'], lchk);
xermsg('SLATEC','LSEI',['INSUFFICIENT STORAGE ',['ALLOCATED FOR IP(*), NEED LIP = ',xern1]],2,1);
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
ws_shape=zeros(ws_shape);ws_shape(:)=ws(1:numel(ws_shape));ws=ws_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
return;
end;
end;
%
% Compute number of possible right multiplying Householder
% transformations.
%
m = fix(me + ma + mg);
if( n<=0 || m<=0 )
mode = 0;
rnorme = 0;
rnorml = 0;
%
elseif( mdw<m ) ;
xermsg('SLATEC','LSEI','MDW.LT.ME+MA+MG IS AN ERROR',2,1);
else;
%
np1 = fix(n + 1);
kranke = fix(min(me,n));
n1 = fix(2.*kranke + 1);
n2 = fix(n1 + n);
%
% Set nominal values.
%
% The nominal column scaling used in the code is
% the identity scaling.
%
[n,dumvar2,dumvar3,ws(sub2ind(size(ws),max(n1,1)):end)]=scopy(n,1.0e0,0,ws(sub2ind(size(ws),max(n1,1)):end),1);
%
% No covariance matrix is nominally computed.
%
cov = false;
%
% Process option vector.
% Define bound for number of options to change.
%
nopt = 1000;
ntimes = 0;
%
% Define bound for positive values of LINK.
%
nlink = 100000;
last = 1;
link = fix(prgopt(1));
if( link==0 || link>nlink )
xermsg('SLATEC','LSEI','THE OPTION VECTOR IS UNDEFINED',2,1);
else;
%
while( link>1 );
ntimes = fix(ntimes + 1);
if( ntimes>nopt )
xermsg('SLATEC','LSEI','THE LINKS IN THE OPTION VECTOR ARE CYCLING.',2,1);
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
ws_shape=zeros(ws_shape);ws_shape(:)=ws(1:numel(ws_shape));ws=ws_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
return;
else;
%
key = fix(prgopt(last+1));
if( key==1 )
cov = prgopt(last+2)~=0.0e0;
elseif( key==2 && prgopt(last+2)~=0.0e0 ) ;
for j = 1 : n;
[t ,m,w(sub2ind(size(w),1,j):end)]=snrm2(m,w(sub2ind(size(w),1,j):end),1);
if( t~=0.0e0 )
t = 1.0e0./t;
end;
ws(j+n1-1) = t;
end; j = fix(n+1);
elseif( key==3 ) ;
[n,prgopt(sub2ind(size(prgopt),max(last+2,1)):end),dumvar3,ws(sub2ind(size(ws),max(n1,1)):end)]=scopy(n,prgopt(sub2ind(size(prgopt),max(last+2,1)):end),1,ws(sub2ind(size(ws),max(n1,1)):end),1);
elseif( key==4 ) ;
tau = max(srelpr,prgopt(last+2));
end;
%
next = fix(prgopt(link));
if( next<=0 || next>nlink )
xermsg('SLATEC','LSEI','THE OPTION VECTOR IS UNDEFINED',2,1);
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
ws_shape=zeros(ws_shape);ws_shape(:)=ws(1:numel(ws_shape));ws=ws_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
return;
else;
%
last = fix(link);
link = fix(next);
end;
end;
end;
%
for j = 1 : n;
[m,ws(n1+j-1),w(sub2ind(size(w),1,j):end)]=sscal(m,ws(n1+j-1),w(sub2ind(size(w),1,j):end),1);
end; j = fix(n+1);
%
if( cov && mdw<n )
xermsg('SLATEC','LSEI','MDW .LT. N WHEN COV MATRIX NEEDED, IS AN ERROR',2,1);
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
ws_shape=zeros(ws_shape);ws_shape(:)=ws(1:numel(ws_shape));ws=ws_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
return;
else;
%
% Problem definition and option vector OK.
%
mode = 0;
%
% Compute norm of equality constraint matrix and right side.
%
enorm = 0.0e0;
for j = 1 : n;
enorm = max(enorm,sasum(me,w(sub2ind(size(w),1,j):end),1));
end; j = fix(n+1);
%
[fnorm ,me,w(sub2ind(size(w),1,np1):end)]=sasum(me,w(sub2ind(size(w),1,np1):end),1);
snmax = 0.0e0;
rnmax = 0.0e0;
for i = 1 : kranke;
%
% Compute maximum ratio of vector lengths. Partition is at
% column I.
%
for k = i : me;
mdw_orig=mdw; [sn ,dumvar2,dumvar2,mdw,dumvar2,dumvar6]=sdot(n-i+1,w(sub2ind(size(w),k,i):end),mdw,w(sub2ind(size(w),k,i):end),mdw); mdw(dumvar6~=mdw_orig)=dumvar6(dumvar6~=mdw_orig); w(sub2ind(size(w),k,i):end)=dumvar2;
mdw_orig=mdw; [rn ,dumvar2,dumvar2,mdw,dumvar2,dumvar6]=sdot(i-1,w(sub2ind(size(w),k,1):end),mdw,w(sub2ind(size(w),k,1):end),mdw); mdw(dumvar6~=mdw_orig)=dumvar6(dumvar6~=mdw_orig); w(sub2ind(size(w),k,1):end)=dumvar2;
if( rn==0.0e0 && sn>snmax )
snmax = sn;
imax = fix(k);
elseif( k==i || sn.*rnmax>rn.*snmax ) ;
snmax = sn;
rnmax = rn;
imax = fix(k);
end;
end; k = fix(me+1);
%
% Interchange rows if necessary.
%
if( i~=imax )
mdw_orig=mdw; [np1,dumvar2,mdw,dumvar4,dumvar5]=sswap(np1,w(sub2ind(size(w),i,1):end),mdw,w(sub2ind(size(w),imax,1):end),mdw); mdw(dumvar5~=mdw_orig)=dumvar5(dumvar5~=mdw_orig); w(sub2ind(size(w),i,1):end)=dumvar2; w(sub2ind(size(w),imax,1):end)=dumvar4;
end;
if( snmax<=rnmax.*tau.^2 )
kranke = fix(i - 1);
break;
else;
%
% Eliminate elements I+1,...,N in row I.
%
mdw_orig=mdw; [dumvar1,i,dumvar3,n,dumvar5,mdw,ws(i),dumvar8,dumvar9]=h12(1,i,i+1,n,w(sub2ind(size(w),i,1):end),mdw,ws(i),w(sub2ind(size(w),i+1,1):end),mdw,1,m-i); mdw(dumvar9~=mdw_orig)=dumvar9(dumvar9~=mdw_orig); w(sub2ind(size(w),i,1):end)=dumvar5; w(sub2ind(size(w),i+1,1):end)=dumvar8;
end;
end;
%
% Save diagonal terms of lower trapezoidal matrix.
%
[kranke,w,dumvar3,ws(sub2ind(size(ws),max(kranke+1,1)):end)]=scopy(kranke,w,mdw+1,ws(sub2ind(size(ws),max(kranke+1,1)):end),1);
%
% use Householder transformation from left to achieve
% KRANKE by KRANKE upper triangular form.
%
if( kranke<me )
for k = kranke : -1: 1 ;
%
% Apply transformation to matrix cols. 1,...,K-1.
%
[dumvar1,k,dumvar3,me,dumvar5,dumvar6,up,w,dumvar9,mdw]=h12(1,k,kranke+1,me,w(sub2ind(size(w),1,k):end),1,up,w,1,mdw,k-1); w(sub2ind(size(w),1,k):end)=dumvar5;
%
% Apply to rt side vector.
%
[dumvar1,k,dumvar3,me,dumvar5,dumvar6,up,dumvar8]=h12(2,k,kranke+1,me,w(sub2ind(size(w),1,k):end),1,up,w(sub2ind(size(w),1,np1):end),1,1,1); w(sub2ind(size(w),1,k):end)=dumvar5; w(sub2ind(size(w),1,np1):end)=dumvar8;
end; k = fix(1 -1);
end;
%
% Solve for variables 1,...,KRANKE in new coordinates.
%
[kranke,w(sub2ind(size(w),1,np1):end),dumvar3,x]=scopy(kranke,w(sub2ind(size(w),1,np1):end),1,x,1);
for i = 1 : kranke;
x(i) =(x(i)-sdot(i-1,w(sub2ind(size(w),i,1):end),mdw,x,1))./w(i,i);
end; i = fix(kranke+1);
%
% Compute residuals for reduced problem.
%
mep1 = fix(me + 1);
rnorml = 0.0e0;
for i = mep1 : m;
w(i,np1) = w(i,np1) - sdot(kranke,w(sub2ind(size(w),i,1):end),mdw,x,1);
mdw_orig=mdw; [sn ,kranke,dumvar2,mdw,dumvar2,dumvar6]=sdot(kranke,w(sub2ind(size(w),i,1):end),mdw,w(sub2ind(size(w),i,1):end),mdw); mdw(dumvar6~=mdw_orig)=dumvar6(dumvar6~=mdw_orig); w(sub2ind(size(w),i,1):end)=dumvar2;
mdw_orig=mdw; [rn ,dumvar2,dumvar2,mdw,dumvar2,dumvar6]=sdot(n-kranke,w(sub2ind(size(w),i,kranke+1):end),mdw,w(sub2ind(size(w),i,kranke+1):end),mdw); mdw(dumvar6~=mdw_orig)=dumvar6(dumvar6~=mdw_orig); w(sub2ind(size(w),i,kranke+1):end)=dumvar2;
if( rn<=sn.*tau.^2 && kranke<n )
[dumvar1,dumvar2,dumvar3,w(sub2ind(size(w),i,kranke+1):end),mdw]=scopy(n-kranke,0.0e0,0,w(sub2ind(size(w),i,kranke+1):end),mdw);
end;
end; i = fix(m+1);
%
% Compute equality constraint equations residual length.
%
[rnorme ,dumvar2,w(sub2ind(size(w),kranke+1,np1):end)]=snrm2(me-kranke,w(sub2ind(size(w),kranke+1,np1):end),1);
%
% Move reduced problem data upward if KRANKE.LT.ME.
%
if( kranke<me )
for j = 1 : np1;
[dumvar1,dumvar2,dumvar3,dumvar4]=scopy(m-me,w(sub2ind(size(w),me+1,j):end),1,w(sub2ind(size(w),kranke+1,j):end),1); w(sub2ind(size(w),me+1,j):end)=dumvar2; w(sub2ind(size(w),kranke+1,j):end)=dumvar4;
end; j = fix(np1+1);
end;
%
% Compute solution of reduced problem.
%
[w(sub2ind(size(w),kranke+1,kranke+1):end),mdw,ma,mg,dumvar5,prgopt,x(sub2ind(size(x),max(kranke+1,1)):end),rnorml,mode,ws(sub2ind(size(ws),max(n2,1)):end),ip(sub2ind(size(ip),max(2,1)):end)]=lsi(w(sub2ind(size(w),kranke+1,kranke+1):end),mdw,ma,mg,n-kranke,prgopt,x(sub2ind(size(x),max(kranke+1,1)):end),rnorml,mode,ws(sub2ind(size(ws),max(n2,1)):end),ip(sub2ind(size(ip),max(2,1)):end));
%
% Test for consistency of equality constraints.
%
gt=0;
if( me>0 )
mdeqc = 0;
[xnrme ,kranke,w(sub2ind(size(w),1,np1):end)]=sasum(kranke,w(sub2ind(size(w),1,np1):end),1);
if( rnorme>tau.*(enorm.*xnrme+fnorm) )
mdeqc = 1;
end;
mode = fix(mode + mdeqc);
%
% Check if solution to equality constraints satisfies inequality
% constraints when there are no degrees of freedom left.
%
if( kranke==n && mg>0 )
[xnorm ,n,x]=sasum(n,x,1);
mapke1 = fix(ma + kranke + 1);
mend = fix(ma + kranke + mg);
for i = mapke1 : mend;
sizemlv = sasum(n,w(sub2ind(size(w),i,1):end),mdw).*xnorm + abs(w(i,np1));
if( w(i,np1)>tau.*sizemlv )
mode = fix(mode + 2);
gt=1;
break;
end;
end;
end;
end;
%
% Replace diagonal terms of lower trapezoidal matrix.
%
if(gt==0)
if( kranke>0 )
[kranke,ws(sub2ind(size(ws),max(kranke+1,1)):end),dumvar3,w]=scopy(kranke,ws(sub2ind(size(ws),max(kranke+1,1)):end),1,w,mdw+1);
%
% Reapply transformation to put solution in original coordinates.
%
for i = kranke : -1: 1 ;
[dumvar1,i,dumvar3,n,w(sub2ind(size(w),i,1):end),mdw,ws(i),x]=h12(2,i,i+1,n,w(sub2ind(size(w),i,1):end),mdw,ws(i),x,1,1,1);
end; i = fix(1 -1);
%
% Compute covariance matrix of equality constrained problem.
%
if( cov )
for j = min(kranke,n-1) : -1: 1 ;
rb = ws(j).*w(j,j);
if( rb~=0.0e0 )
rb = 1.0e0./rb;
end;
jp1 = fix(j + 1);
for i = jp1 : n;
w(i,j) = rb.*sdot(n-j,w(sub2ind(size(w),i,jp1):end),mdw,w(sub2ind(size(w),j,jp1):end),mdw);
end; i = fix(n+1);
%
gam = 0.5e0.*rb.*sdot(n-j,w(sub2ind(size(w),jp1,j):end),1,w(sub2ind(size(w),j,jp1):end),mdw);
[dumvar1,gam,dumvar3,mdw,dumvar5]=saxpy(n-j,gam,w(sub2ind(size(w),j,jp1):end),mdw,w(sub2ind(size(w),jp1,j):end),1); w(sub2ind(size(w),j,jp1):end)=dumvar3; w(sub2ind(size(w),jp1,j):end)=dumvar5;
for i = jp1 : n;
for k = i : n;
w(i,k) = w(i,k) + w(j,i).*w(k,j) + w(i,j).*w(j,k);
w(k,i) = w(i,k);
end; k = fix(n+1);
end; i = fix(n+1);
uj = ws(j);
vj = gam.*uj;
w(j,j) = uj.*vj + uj.*vj;
for i = jp1 : n;
w(j,i) = uj.*w(i,j) + vj.*w(j,i);
end; i = fix(n+1);
[dumvar1,dumvar2,mdw,dumvar4]=scopy(n-j,w(sub2ind(size(w),j,jp1):end),mdw,w(sub2ind(size(w),jp1,j):end),1); w(sub2ind(size(w),j,jp1):end)=dumvar2; w(sub2ind(size(w),jp1,j):end)=dumvar4;
end; j = fix(1 -1);
end;
end;
%
% Apply the scaling to the covariance matrix.
%
if( cov )
for i = 1 : n;
[n,ws(i+n1-1),w(sub2ind(size(w),i,1):end),mdw]=sscal(n,ws(i+n1-1),w(sub2ind(size(w),i,1):end),mdw);
[n,ws(i+n1-1),w(sub2ind(size(w),1,i):end)]=sscal(n,ws(i+n1-1),w(sub2ind(size(w),1,i):end),1);
end; i = fix(n+1);
end;
end;
%
% Rescale solution vector.
%
if( mode<=1 )
for j = 1 : n;
x(j) = x(j).*ws(n1+j-1);
end; j = fix(n+1);
end;
%
ip(1) = fix(kranke);
ip(3) = fix(ip(3) + 2.*kranke + n);
end;
end;
end;
end;
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
ws_shape=zeros(ws_shape);ws_shape(:)=ws(1:numel(ws_shape));ws=ws_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
end %subroutine lsei
%DECK LSI
|
|
