[x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info]=sqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info);

function [x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info]=sqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info);
%***BEGIN PROLOGUE  SQRSL
%***PURPOSE  Apply the output of SQRDC to compute coordinate transfor-
%            mations, projections, and least squares solutions.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D9, D2A1
%***TYPE      SINGLE PRECISION (SQRSL-S, DQRSL-D, CQRSL-C)
%***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX, ORTHOGONAL TRIANGULAR,
%             SOLVE
%***AUTHOR  Stewart, G. W., (U. of Maryland)
%***DESCRIPTION
%
%     SQRSL applies the output of SQRDC to compute coordinate
%     transformations, projections, and least squares solutions.
%     For K .LE. MIN(N,P), let XK be the matrix
%
%            XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K)))
%
%     formed from columns JPVT(1), ... ,JPVT(K) of the original
%     N x P matrix X that was input to SQRDC (if no pivoting was
%     done, XK consists of the first K columns of X in their
%     original order).  SQRDC produces a factored orthogonal matrix Q
%     and an upper triangular matrix R such that
%
%              XK = Q * (R)
%                       (0)
%
%     This information is contained in coded form in the arrays
%     X and QRAUX.
%
%     On Entry
%
%        X      REAL(LDX,P)
%               X contains the output of SQRDC.
%
%        LDX    INTEGER
%               LDX is the leading dimension of the array X.
%
%        N      INTEGER
%               N is the number of rows of the matrix XK.  It must
%               have the same value as N in SQRDC.
%
%        K      INTEGER
%               K is the number of columns of the matrix XK.  K
%               must not be greater than MIN(N,P), where P is the
%               same as in the calling sequence to SQRDC.
%
%        QRAUX  REAL(P)
%               QRAUX contains the auxiliary output from SQRDC.
%
%        Y      REAL(N)
%               Y contains an N-vector that is to be manipulated
%               by SQRSL.
%
%        JOB    INTEGER
%               JOB specifies what is to be computed.  JOB has
%               the decimal expansion ABCDE, with the following
%               meaning.
%
%                    If A .NE. 0, compute QY.
%                    If B,C,D, or E .NE. 0, compute QTY.
%                    If C .NE. 0, compute B.
%                    If D .NE. 0, compute RSD.
%                    If E .NE. 0, compute XB.
%
%               Note that a request to compute B, RSD, or XB
%               automatically triggers the computation of QTY, for
%               which an array must be provided in the calling
%               sequence.
%
%     On Return
%
%        QY     REAL(N).
%               QY contains Q*Y, if its computation has been
%               requested.
%
%        QTY    REAL(N).
%               QTY contains TRANS(Q)*Y, if its computation has
%               been requested.  Here TRANS(Q) is the
%               transpose of the matrix Q.
%
%        B      REAL(K)
%               B contains the solution of the least squares problem
%
%                    minimize norm2(Y - XK*B),
%
%               if its computation has been requested.  (Note that
%               if pivoting was requested in SQRDC, the J-th
%               component of B will be associated with column JPVT(J)
%               of the original matrix X that was input into SQRDC.)
%
%        RSD    REAL(N).
%               RSD contains the least squares residual Y - XK*B,
%               if its computation has been requested.  RSD is
%               also the orthogonal projection of Y onto the
%               orthogonal complement of the column space of XK.
%
%        XB     REAL(N).
%               XB contains the least squares approximation XK*B,
%               if its computation has been requested.  XB is also
%               the orthogonal projection of Y onto the column space
%               of X.
%
%        INFO   INTEGER.
%               INFO is zero unless the computation of B has
%               been requested and R is exactly singular.  In
%               this case, INFO is the index of the first zero
%               diagonal element of R and B is left unaltered.
%
%     The parameters QY, QTY, B, RSD, and XB are not referenced
%     if their computation is not requested and in this case
%     can be replaced by dummy variables in the calling program.
%     To savemlv storage, the user may in some cases use the same
%     array for different parameters in the calling sequence.  A
%     frequently occurring example is when one wishes to compute
%     any of B, RSD, or XB and does not need Y or QTY.  In this
%     case one may identify Y, QTY, and one of B, RSD, or XB, while
%     providing separate arrays for anything else that is to be
%     computed.  Thus the calling sequence
%
%          CALL SQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO)
%
%     will result in the computation of B and RSD, with RSD
%     overwriting Y.  More generally, each item in the following
%     list contains groups of permissible identifications for
%     a single calling sequence.
%
%          1. (Y,QTY,B) (RSD) (XB) (QY)
%
%          2. (Y,QTY,RSD) (B) (XB) (QY)
%
%          3. (Y,QTY,XB) (B) (RSD) (QY)
%
%          4. (Y,QY) (QTY,B) (RSD) (XB)
%
%          5. (Y,QY) (QTY,RSD) (B) (XB)
%
%          6. (Y,QY) (QTY,XB) (B) (RSD)
%
%     In any group the value returned in the array allocated to
%     the group corresponds to the last member of the group.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  SAXPY, SCOPY, SDOT
%***REVISION HISTORY  (YYMMDD)
%   780814  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  SQRSL
persistent cb cqty cqy cr cxb i j jj ju kp1 t temp ; 

x_shape=size(x);x=reshape([x(:).',zeros(1,ceil(numel(x)./prod([ldx])).*prod([ldx])-numel(x))],ldx,[]);
qraux_shape=size(qraux);qraux=reshape(qraux,1,[]);
y_shape=size(y);y=reshape(y,1,[]);
qy_shape=size(qy);qy=reshape(qy,1,[]);
qty_shape=size(qty);qty=reshape(qty,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
rsd_shape=size(rsd);rsd=reshape(rsd,1,[]);
xb_shape=size(xb);xb=reshape(xb,1,[]);
%
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(ju), ju=0; end;
if isempty(kp1), kp1=0; end;
if isempty(t), t=0; end;
if isempty(temp), temp=0; end;
if isempty(cb), cb=false; end;
if isempty(cqy), cqy=false; end;
if isempty(cqty), cqty=false; end;
if isempty(cr), cr=false; end;
if isempty(cxb), cxb=false; end;
%***FIRST EXECUTABLE STATEMENT  SQRSL
%
%     SET INFO FLAG.
%
info = 0;
%
%     DETERMINE WHAT IS TO BE COMPUTED.
%
cqy = fix(job./10000)~=0;
cqty = rem(job,10000)~=0;
cb = rem(job,1000)./100~=0;
cr = rem(job,100)./10~=0;
cxb = rem(job,10)~=0;
ju = fix(min(k,n-1));
%
%     SPECIAL ACTION WHEN N=1.
%
if( ju~=0 )
%
%        SET UP TO COMPUTE QY OR QTY.
%
if( cqy )
[n,y,dumvar3,qy]=scopy(n,y,1,qy,1);
end;
if( cqty )
[n,y,dumvar3,qty]=scopy(n,y,1,qty,1);
end;
if( cqy )
%
%           COMPUTE QY.
%
for jj = 1 : ju;
j = fix(ju - jj + 1);
if( qraux(j)~=0.0e0 )
temp = x(j,j);
x(j,j) = qraux(j);
t = -sdot(n-j+1,x(sub2ind(size(x),j,j):end),1,qy(sub2ind(size(qy),max(j,1)):end),1)./x(j,j);
[dumvar1,t,x(sub2ind(size(x),j,j):end),dumvar4,qy(sub2ind(size(qy),max(j,1)):end)]=saxpy(n-j+1,t,x(sub2ind(size(x),j,j):end),1,qy(sub2ind(size(qy),max(j,1)):end),1);
x(j,j) = temp;
end;
end; jj = fix(ju+1);
end;
if( cqty )
%
%           COMPUTE TRANS(Q)*Y.
%
for j = 1 : ju;
if( qraux(j)~=0.0e0 )
temp = x(j,j);
x(j,j) = qraux(j);
t = -sdot(n-j+1,x(sub2ind(size(x),j,j):end),1,qty(sub2ind(size(qty),max(j,1)):end),1)./x(j,j);
[dumvar1,t,x(sub2ind(size(x),j,j):end),dumvar4,qty(sub2ind(size(qty),max(j,1)):end)]=saxpy(n-j+1,t,x(sub2ind(size(x),j,j):end),1,qty(sub2ind(size(qty),max(j,1)):end),1);
x(j,j) = temp;
end;
end; j = fix(ju+1);
end;
%
%        SET UP TO COMPUTE B, RSD, OR XB.
%
if( cb )
[k,qty,dumvar3,b]=scopy(k,qty,1,b,1);
end;
kp1 = fix(k + 1);
if( cxb )
[k,qty,dumvar3,xb]=scopy(k,qty,1,xb,1);
end;
if( cr && k<n )
[dumvar1,qty(sub2ind(size(qty),max(kp1,1)):end),dumvar3,rsd(sub2ind(size(rsd),max(kp1,1)):end)]=scopy(n-k,qty(sub2ind(size(qty),max(kp1,1)):end),1,rsd(sub2ind(size(rsd),max(kp1,1)):end),1);
end;
if( ~(~cxb || kp1>n) )
for i = kp1 : n;
xb(i) = 0.0e0;
end; i = fix(n+1);
end;
if( cr )
for i = 1 : k;
rsd(i) = 0.0e0;
end; i = fix(k+1);
end;
if( cb )
%
%           COMPUTE B.
%
for jj = 1 : k;
j = fix(k - jj + 1);
if( x(j,j)~=0.0e0 )
b(j) = b(j)./x(j,j);
if( j~=1 )
t = -b(j);
[dumvar1,t,x(sub2ind(size(x),1,j):end),dumvar4,b]=saxpy(j-1,t,x(sub2ind(size(x),1,j):end),1,b,1);
end;
else;
info = fix(j);
break;
end;
end;
end;
if( ~(~cr && ~cxb) )
%
%           COMPUTE RSD OR XB AS REQUIRED.
%
for jj = 1 : ju;
j = fix(ju - jj + 1);
if( qraux(j)~=0.0e0 )
temp = x(j,j);
x(j,j) = qraux(j);
if( cr )
t = -sdot(n-j+1,x(sub2ind(size(x),j,j):end),1,rsd(sub2ind(size(rsd),max(j,1)):end),1)./x(j,j);
[dumvar1,t,x(sub2ind(size(x),j,j):end),dumvar4,rsd(sub2ind(size(rsd),max(j,1)):end)]=saxpy(n-j+1,t,x(sub2ind(size(x),j,j):end),1,rsd(sub2ind(size(rsd),max(j,1)):end),1);
end;
if( cxb )
t = -sdot(n-j+1,x(sub2ind(size(x),j,j):end),1,xb(sub2ind(size(xb),max(j,1)):end),1)./x(j,j);
[dumvar1,t,x(sub2ind(size(x),j,j):end),dumvar4,xb(sub2ind(size(xb),max(j,1)):end)]=saxpy(n-j+1,t,x(sub2ind(size(x),j,j):end),1,xb(sub2ind(size(xb),max(j,1)):end),1);
end;
x(j,j) = temp;
end;
end; jj = fix(ju+1);
end;
else;
if( cqy )
qy(1) = y(1);
end;
if( cqty )
qty(1) = y(1);
end;
if( cxb )
xb(1) = y(1);
end;
if( cb )
if( x(1,1)~=0.0e0 )
b(1) = y(1)./x(1,1);
else;
info = 1;
end;
end;
if( cr )
rsd(1) = 0.0e0;
end;
end;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
qraux_shape=zeros(qraux_shape);qraux_shape(:)=qraux(1:numel(qraux_shape));qraux=qraux_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
qy_shape=zeros(qy_shape);qy_shape(:)=qy(1:numel(qy_shape));qy=qy_shape;
qty_shape=zeros(qty_shape);qty_shape(:)=qty(1:numel(qty_shape));qty=qty_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
rsd_shape=zeros(rsd_shape);rsd_shape(:)=rsd(1:numel(rsd_shape));rsd=rsd_shape;
xb_shape=zeros(xb_shape);xb_shape(:)=xb(1:numel(xb_shape));xb=xb_shape;
end
%DECK SREADP