[a,mda,m,n,b,mdb,nb,mode,krank,rnorm,h,w,ir,ic]=du12us(a,mda,m,n,b,mdb,nb,mode,krank,rnorm,h,w,ir,ic);

function [a,mda,m,n,b,mdb,nb,mode,krank,rnorm,h,w,ir,ic]=du12us(a,mda,m,n,b,mdb,nb,mode,krank,rnorm,h,w,ir,ic);
%***BEGIN PROLOGUE  DU12US
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DULSIA
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (U12US-S, DU12US-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%        Given the Householder LQ factorization of A, this
%        subroutine solves the system AX=B. If the system
%        is of reduced rank, this routine returns a solution
%        according to the selected mode.
%
%       Note - If MODE.NE.2, W is never accessed.
%
%***SEE ALSO  DULSIA
%***ROUTINES CALLED  DAXPY, DDOT, DNRM2, DSWAP
%***REVISION HISTORY  (YYMMDD)
%   810801  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  DU12US
persistent bb i ij ip1 j jb k kp1 mmk tt ; 

if isempty(bb), bb=0; end;
if isempty(tt), tt=0; end;
if isempty(i), i=0; end;
if isempty(ij), ij=0; end;
if isempty(ip1), ip1=0; end;
if isempty(j), j=0; end;
if isempty(jb), jb=0; end;
if isempty(k), k=0; end;
if isempty(kp1), kp1=0; end;
if isempty(mmk), mmk=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([mda])).*prod([mda])-numel(a))],mda,[]);
b_shape=size(b);b=reshape([b(:).',zeros(1,ceil(numel(b)./prod([mdb])).*prod([mdb])-numel(b))],mdb,[]);
rnorm_shape=size(rnorm);rnorm=reshape(rnorm,1,[]);
h_shape=size(h);h=reshape(h,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
ic_shape=size(ic);ic=reshape(ic,1,[]);
ir_shape=size(ir);ir=reshape(ir,1,[]);
%***FIRST EXECUTABLE STATEMENT  DU12US
k = fix(krank);
kp1 = fix(k + 1);
%
%        RANK=0
%
if( k>0 )
%
%     REORDER B TO REFLECT ROW INTERCHANGES
%
i = 0;
while( true );
i = fix(i + 1);
if( i==m )
break;
end;
j = fix(ir(i));
if( j~=i )
if( j>=0 )
ir(i) = fix(-ir(i));
for jb = 1 : nb;
rnorm(jb) = b(i,jb);
end; jb = fix(nb+1);
ij = fix(i);
while( true );
for jb = 1 : nb;
b(ij,jb) = b(j,jb);
end; jb = fix(nb+1);
ij = fix(j);
j = fix(ir(ij));
ir(ij) = fix(-ir(ij));
if( j==i )
break;
end;
end;
for jb = 1 : nb;
b(ij,jb) = rnorm(jb);
end; jb = fix(nb+1);
end;
end;
end;
for i = 1 : m;
ir(i) = fix(abs(ir(i)));
end; i = fix(m+1);
%
%     IF A IS OF REDUCED RANK AND MODE=2,
%     APPLY HOUSEHOLDER TRANSFORMATIONS TO B
%
if( mode>=2 && k~=m )
mmk = fix(m - k);
for jb = 1 : nb;
for j = 1 : k;
i = fix(kp1 - j);
tt = -ddot(mmk,a(sub2ind(size(a),kp1,i):end),1,b(sub2ind(size(b),kp1,jb):end),1)./w(i);
tt = tt - b(i,jb);
[mmk,tt,a(sub2ind(size(a),kp1,i):end),dumvar4,b(sub2ind(size(b),kp1,jb):end)]=daxpy(mmk,tt,a(sub2ind(size(a),kp1,i):end),1,b(sub2ind(size(b),kp1,jb):end),1);
b(i,jb) = b(i,jb) + tt.*w(i);
end; j = fix(k+1);
end; jb = fix(nb+1);
end;
%
%     FIND NORMS OF RESIDUAL VECTOR(S)..(BEFORE OVERWRITE B)
%
for jb = 1 : nb;
[rnorm(jb) ,dumvar2,b(sub2ind(size(b),kp1,jb):end)]=dnrm2((m-k),b(sub2ind(size(b),kp1,jb):end),1);
end; jb = fix(nb+1);
%
%     BACK SOLVE LOWER TRIANGULAR L
%
for jb = 1 : nb;
for i = 1 : k;
b(i,jb) = b(i,jb)./a(i,i);
if( i==k )
break;
end;
ip1 = fix(i + 1);
[dumvar1,dumvar2,a(sub2ind(size(a),ip1,i):end),dumvar4,dumvar5]=daxpy(k-i,--b(i,jb),a(sub2ind(size(a),ip1,i):end),1,b(sub2ind(size(b),ip1,jb):end),1);      -b(i,jb)=dumvar2; b(sub2ind(size(b),ip1,jb):end)=dumvar5; 
end;
end;
%
%
%      TRUNCATED SOLUTION
%
if( k~=n )
for jb = 1 : nb;
for i = kp1 : n;
b(i,jb) = 0.0d0;
end; i = fix(n+1);
end; jb = fix(nb+1);
end;
%
%     APPLY HOUSEHOLDER TRANSFORMATIONS TO B
%
for i = 1 : k;
j = fix(kp1 - i);
tt = a(j,j);
a(j,j) = h(j);
for jb = 1 : nb;
bb = -ddot(n-j+1,a(sub2ind(size(a),j,j):end),mda,b(sub2ind(size(b),j,jb):end),1)./h(j);
[dumvar1,bb,a(sub2ind(size(a),j,j):end),mda,b(sub2ind(size(b),j,jb):end)]=daxpy(n-j+1,bb,a(sub2ind(size(a),j,j):end),mda,b(sub2ind(size(b),j,jb):end),1);
end; jb = fix(nb+1);
a(j,j) = tt;
end; i = fix(k+1);
%
%
%     REORDER B TO REFLECT COLUMN INTERCHANGES
%
i = 0;
while( true );
i = fix(i + 1);
if( i==n )
break;
end;
j = fix(ic(i));
if( j~=i )
if( j>=0 )
ic(i) = fix(-ic(i));
while( true );
mdb_orig=mdb;    [nb,dumvar2,mdb,dumvar4,dumvar5]=dswap(nb,b(sub2ind(size(b),j,1):end),mdb,b(sub2ind(size(b),i,1):end),mdb);    mdb(dumvar5~=mdb_orig)=dumvar5(dumvar5~=mdb_orig);      b(sub2ind(size(b),j,1):end)=dumvar2; b(sub2ind(size(b),i,1):end)=dumvar4; 
ij = fix(ic(j));
ic(j) = fix(-ic(j));
j = fix(ij);
if( j==i )
break;
end;
end;
end;
end;
end;
for i = 1 : n;
ic(i) = fix(abs(ic(i)));
end; i = fix(n+1);
else;
for jb = 1 : nb;
[rnorm(jb) ,m,b(sub2ind(size(b),1,jb):end)]=dnrm2(m,b(sub2ind(size(b),1,jb):end),1);
end; jb = fix(nb+1);
for jb = 1 : nb;
for i = 1 : n;
b(i,jb) = 0.0d0;
end; i = fix(n+1);
end; jb = fix(nb+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
ic_shape=zeros(ic_shape);ic_shape(:)=ic(1:numel(ic_shape));ic=ic_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
return;
end;
%
%        SOLUTION VECTORS ARE IN FIRST N ROWS OF B()
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
ic_shape=zeros(ic_shape);ic_shape(:)=ic(1:numel(ic_shape));ic=ic_shape;
ir_shape=zeros(ir_shape);ir_shape(:)=ir(1:numel(ir_shape));ir=ir_shape;
end %subroutine du12us
%DECK DUIVP