Code covered by the BSD License

### Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[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,
%***AUTHOR  Hanson, R. J., (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
%
%                  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(2) = KEY1 (key to the option change)
%               .  PRGOPT(3) = data value (data value for this change)
%               .       .
%               .       .
%               .       .
%               .                       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.
%                  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
%
%                   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(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(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.
%
last = 1;
xermsg('SLATEC','LSEI','THE OPTION VECTOR IS UNDEFINED',2,1);
else;
%
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;
%
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;
%
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
```