| [a,mda,m,n,ub,db,mode,np,krank,ksure,h,w,eb,ir,ic]=du11us(a,mda,m,n,ub,db,mode,np,krank,ksure,h,w,eb,ir,ic); |
function [a,mda,m,n,ub,db,mode,np,krank,ksure,h,w,eb,ir,ic]=du11us(a,mda,m,n,ub,db,mode,np,krank,ksure,h,w,eb,ir,ic);
%***BEGIN PROLOGUE DU11US
%***SUBSIDIARY
%***PURPOSE Subsidiary to DULSIA
%***LIBRARY SLATEC
%***TYPE doubleprecision (U11US-S, DU11US-D)
%***AUTHOR (UNKNOWN)
%***DESCRIPTION
%
% This routine performs an LQ factorization of the
% matrix A using Householder transformations. Row
% and column pivots are chosen to reduce the growth
% of round-off and to help detect possible rank
% deficiency.
%
%***SEE ALSO DULSIA
%***ROUTINES CALLED DAXPY, DDOT, DNRM2, DSCAL, DSWAP, IDAMAX, ISWAP,
% XERMSG
%***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)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
% 900328 Added TYPE section. (WRB)
%***end PROLOGUE DU11US
persistent bb gt gt2 i ii im1 imin is j jm1 jmax jp1 kk km1 kmi kp1 kz l lm1 mmk nn r2 rmin summlv t temp tn tt ;
if isempty(bb), bb=0; end;
if isempty(r2), r2=0; end;
if isempty(rmin), rmin=0; end;
if isempty(summlv), summlv=0; end;
if isempty(t), t=0; end;
if isempty(temp), temp=0; end;
if isempty(tn), tn=0; end;
if isempty(tt), tt=0; end;
if isempty(i), i=0; end;
if isempty(ii), ii=0; end;
if isempty(im1), im1=0; end;
if isempty(imin), imin=0; end;
if isempty(is), is=0; end;
if isempty(j), j=0; end;
if isempty(jm1), jm1=0; end;
if isempty(jmax), jmax=0; end;
if isempty(jp1), jp1=0; end;
if isempty(kk), kk=0; end;
if isempty(km1), km1=0; end;
if isempty(kmi), kmi=0; end;
if isempty(kp1), kp1=0; end;
if isempty(kz), kz=0; end;
if isempty(l), l=0; end;
if isempty(lm1), lm1=0; end;
if isempty(mmk), mmk=0; end;
if isempty(nn), nn=0; end;
if isempty(gt), gt=0; end;
if isempty(gt2), gt2=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([mda])).*prod([mda])-numel(a))],mda,[]);
ub_shape=size(ub);ub=reshape(ub,1,[]);
db_shape=size(db);db=reshape(db,1,[]);
h_shape=size(h);h=reshape(h,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
eb_shape=size(eb);eb=reshape(eb,1,[]);
ic_shape=size(ic);ic=reshape(ic,1,[]);
ir_shape=size(ir);ir=reshape(ir,1,[]);
%
% INITIALIZATION
%
%***FIRST EXECUTABLE STATEMENT DU11US
j = 0;
krank = fix(m);
for i = 1 : n;
ic(i) = fix(i);
end; i = fix(n+1);
for i = 1 : m;
ir(i) = fix(i);
end; i = fix(m+1);
%
% DETERMINE REL AND ABS ERROR VECTORS
%
%
%
% CALCULATE ROW LENGTH
%
for i = 1 : m;
[h(i) ,n,a(sub2ind(size(a),i,1):end),mda]=dnrm2(n,a(sub2ind(size(a),i,1):end),mda);
w(i) = h(i);
end; i = fix(m+1);
%
% INITIALIZE ERROR BOUNDS
%
for i = 1 : m;
eb(i) = max(db(i),ub(i).*h(i));
ub(i) = eb(i);
db(i) = 0.0d0;
end; i = fix(m+1);
%
% DISCARD SELF DEPENDENT ROWS
%
i = 1;
while( true );
if( eb(i)>=h(i) )
%
% MATRIX REDUCTION
%
kk = fix(krank);
krank = fix(krank - 1);
if( mode==0 )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
ub_shape=zeros(ub_shape);ub_shape(:)=ub(1:numel(ub_shape));ub=ub_shape;
db_shape=zeros(db_shape);db_shape(:)=db(1:numel(db_shape));db=db_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;
eb_shape=zeros(eb_shape);eb_shape(:)=eb(1:numel(eb_shape));eb=eb_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;
if( i<=np )
xermsg('SLATEC','DU11US','FIRST NP ROWS ARE LINEARLY DEPENDENT',8,0);
krank = fix(i - 1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
ub_shape=zeros(ub_shape);ub_shape(:)=ub(1:numel(ub_shape));ub=ub_shape;
db_shape=zeros(db_shape);db_shape(:)=db(1:numel(db_shape));db=db_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;
eb_shape=zeros(eb_shape);eb_shape(:)=eb(1:numel(eb_shape));eb=eb_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;
else;
if( i>krank )
break;
end;
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,eb(sub2ind(size(eb),max(i,1)):end),1,eb(sub2ind(size(eb),max(kk,1)):end),1); dumvar2i=find((eb(sub2ind(size(eb),max(i,1)):end))~=(dumvar2));dumvar4i=find((eb(sub2ind(size(eb),max(kk,1)):end))~=(dumvar4)); eb(i-1+dumvar2i)=dumvar2(dumvar2i); eb(kk-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,ub(sub2ind(size(ub),max(i,1)):end),1,ub(sub2ind(size(ub),max(kk,1)):end),1); dumvar2i=find((ub(sub2ind(size(ub),max(i,1)):end))~=(dumvar2));dumvar4i=find((ub(sub2ind(size(ub),max(kk,1)):end))~=(dumvar4)); ub(i-1+dumvar2i)=dumvar2(dumvar2i); ub(kk-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,w(sub2ind(size(w),max(i,1)):end),1,w(sub2ind(size(w),max(kk,1)):end),1); dumvar2i=find((w(sub2ind(size(w),max(i,1)):end))~=(dumvar2));dumvar4i=find((w(sub2ind(size(w),max(kk,1)):end))~=(dumvar4)); w(i-1+dumvar2i)=dumvar2(dumvar2i); w(kk-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,h(sub2ind(size(h),max(i,1)):end),1,h(sub2ind(size(h),max(kk,1)):end),1); dumvar2i=find((h(sub2ind(size(h),max(i,1)):end))~=(dumvar2));dumvar4i=find((h(sub2ind(size(h),max(kk,1)):end))~=(dumvar4)); h(i-1+dumvar2i)=dumvar2(dumvar2i); h(kk-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=iswap(1,ir(sub2ind(size(ir),max(i,1)):end),1,ir(sub2ind(size(ir),max(kk,1)):end),1); dumvar2i=find((ir(sub2ind(size(ir),max(i,1)):end))~=(dumvar2));dumvar4i=find((ir(sub2ind(size(ir),max(kk,1)):end))~=(dumvar4)); ir(i-1+dumvar2i)=dumvar2(dumvar2i); ir(kk-1+dumvar4i)=dumvar4(dumvar4i);
mda_orig=mda; [n,dumvar2,mda,dumvar4,dumvar5]=dswap(n,a(sub2ind(size(a),i,1):end),mda,a(sub2ind(size(a),kk,1):end),mda); mda(dumvar5~=mda_orig)=dumvar5(dumvar5~=mda_orig); a(sub2ind(size(a),i,1):end)=dumvar2; a(sub2ind(size(a),kk,1):end)=dumvar4;
end;
elseif( i==krank ) ;
break;
else;
i = fix(i + 1);
end;
end;
%
% TEST FOR ZERO RANK
%
if( krank<=0 )
krank = 0;
ksure = 0;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
ub_shape=zeros(ub_shape);ub_shape(:)=ub(1:numel(ub_shape));ub=ub_shape;
db_shape=zeros(db_shape);db_shape(:)=db(1:numel(db_shape));db=db_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;
eb_shape=zeros(eb_shape);eb_shape(:)=eb(1:numel(eb_shape));eb=eb_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;
else;
%
% M A I N L O O P
%
gt2=0;
while( true );
j = fix(j + 1);
jp1 = fix(j + 1);
jm1 = fix(j - 1);
kz = fix(krank);
if( j<=np )
kz = fix(j);
end;
%
% EACH ROW HAS NN=N-J+1 COMPONENTS
%
nn = fix(n - j + 1);
%
% UB DETERMINES ROW PIVOT
%
while( true );
imin = fix(j);
gt=0;
if( h(j)~=0.0d0 )
rmin = ub(j)./h(j);
for i = j : kz;
if( ub(i)<h(i).*rmin )
rmin = ub(i)./h(i);
imin = fix(i);
end;
end; i = fix(kz+1);
%
% TEST FOR RANK DEFICIENCY
%
if( rmin>=1.0d0 )
tt =(eb(imin)+abs(db(imin)))./h(imin);
if( tt>=1.0d0 )
gt=1;
end;
% COMPUTE EXACT UB
if(gt==0)
for i = 1 : jm1;
w(i) = a(imin,i);
end; i = fix(jm1+1);
l = fix(jm1);
while( true );
w(l) = w(l)./a(l,l);
if( l==1 )
break;
end;
lm1 = fix(l - 1);
for i = l : jm1;
w(lm1) = w(lm1) - a(i,lm1).*w(i);
end; i = fix(jm1+1);
l = fix(lm1);
end;
tt = eb(imin);
for i = 1 : jm1;
tt = tt + abs(w(i)).*eb(i);
end; i = fix(jm1+1);
ub(imin) = tt;
if( ub(imin)./h(imin)>=1.0d0 )
gt=1;
end;
end;
end;
if(gt==0)
%
% ROW PIVOT
%
if( imin~=j )
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,h(sub2ind(size(h),max(j,1)):end),1,h(sub2ind(size(h),max(imin,1)):end),1); dumvar2i=find((h(sub2ind(size(h),max(j,1)):end))~=(dumvar2));dumvar4i=find((h(sub2ind(size(h),max(imin,1)):end))~=(dumvar4)); h(j-1+dumvar2i)=dumvar2(dumvar2i); h(imin-1+dumvar4i)=dumvar4(dumvar4i);
mda_orig=mda; [n,dumvar2,mda,dumvar4,dumvar5]=dswap(n,a(sub2ind(size(a),j,1):end),mda,a(sub2ind(size(a),imin,1):end),mda); mda(dumvar5~=mda_orig)=dumvar5(dumvar5~=mda_orig); a(sub2ind(size(a),j,1):end)=dumvar2; a(sub2ind(size(a),imin,1):end)=dumvar4;
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,eb(sub2ind(size(eb),max(j,1)):end),1,eb(sub2ind(size(eb),max(imin,1)):end),1); dumvar2i=find((eb(sub2ind(size(eb),max(j,1)):end))~=(dumvar2));dumvar4i=find((eb(sub2ind(size(eb),max(imin,1)):end))~=(dumvar4)); eb(j-1+dumvar2i)=dumvar2(dumvar2i); eb(imin-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,ub(sub2ind(size(ub),max(j,1)):end),1,ub(sub2ind(size(ub),max(imin,1)):end),1); dumvar2i=find((ub(sub2ind(size(ub),max(j,1)):end))~=(dumvar2));dumvar4i=find((ub(sub2ind(size(ub),max(imin,1)):end))~=(dumvar4)); ub(j-1+dumvar2i)=dumvar2(dumvar2i); ub(imin-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,db(sub2ind(size(db),max(j,1)):end),1,db(sub2ind(size(db),max(imin,1)):end),1); dumvar2i=find((db(sub2ind(size(db),max(j,1)):end))~=(dumvar2));dumvar4i=find((db(sub2ind(size(db),max(imin,1)):end))~=(dumvar4)); db(j-1+dumvar2i)=dumvar2(dumvar2i); db(imin-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,w(sub2ind(size(w),max(j,1)):end),1,w(sub2ind(size(w),max(imin,1)):end),1); dumvar2i=find((w(sub2ind(size(w),max(j,1)):end))~=(dumvar2));dumvar4i=find((w(sub2ind(size(w),max(imin,1)):end))~=(dumvar4)); w(j-1+dumvar2i)=dumvar2(dumvar2i); w(imin-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=iswap(1,ir(sub2ind(size(ir),max(j,1)):end),1,ir(sub2ind(size(ir),max(imin,1)):end),1); dumvar2i=find((ir(sub2ind(size(ir),max(j,1)):end))~=(dumvar2));dumvar4i=find((ir(sub2ind(size(ir),max(imin,1)):end))~=(dumvar4)); ir(j-1+dumvar2i)=dumvar2(dumvar2i); ir(imin-1+dumvar4i)=dumvar4(dumvar4i);
end;
%
% COLUMN PIVOT
%
[jmax ,nn,a(sub2ind(size(a),j,j):end),mda]=idamax(nn,a(sub2ind(size(a),j,j):end),mda);
jmax = fix(jmax + j - 1);
if( jmax~=j )
[m,dumvar2,dumvar3,dumvar4]=dswap(m,a(sub2ind(size(a),1,j):end),1,a(sub2ind(size(a),1,jmax):end),1); a(sub2ind(size(a),1,j):end)=dumvar2; a(sub2ind(size(a),1,jmax):end)=dumvar4;
[dumvar1,dumvar2,dumvar3,dumvar4]=iswap(1,ic(sub2ind(size(ic),max(j,1)):end),1,ic(sub2ind(size(ic),max(jmax,1)):end),1); dumvar2i=find((ic(sub2ind(size(ic),max(j,1)):end))~=(dumvar2));dumvar4i=find((ic(sub2ind(size(ic),max(jmax,1)):end))~=(dumvar4)); ic(j-1+dumvar2i)=dumvar2(dumvar2i); ic(jmax-1+dumvar4i)=dumvar4(dumvar4i);
end;
%
% APPLY HOUSEHOLDER TRANSFORMATION
%
[tn ,nn,a(sub2ind(size(a),j,j):end),mda]=dnrm2(nn,a(sub2ind(size(a),j,j):end),mda);
if( tn~=0.0d0 )
break;
end;
end;
end;
%
% MATRIX REDUCTION
%
kk = fix(krank);
krank = fix(krank - 1);
kz = fix(krank);
if( mode==0 )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
ub_shape=zeros(ub_shape);ub_shape(:)=ub(1:numel(ub_shape));ub=ub_shape;
db_shape=zeros(db_shape);db_shape(:)=db(1:numel(db_shape));db=db_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;
eb_shape=zeros(eb_shape);eb_shape(:)=eb(1:numel(eb_shape));eb=eb_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;
if( j<=np )
xermsg('SLATEC','DU11US','FIRST NP ROWS ARE LINEARLY DEPENDENT',8,0);
krank = fix(j - 1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
ub_shape=zeros(ub_shape);ub_shape(:)=ub(1:numel(ub_shape));ub=ub_shape;
db_shape=zeros(db_shape);db_shape(:)=db(1:numel(db_shape));db=db_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;
eb_shape=zeros(eb_shape);eb_shape(:)=eb(1:numel(eb_shape));eb=eb_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;
if( imin<=krank )
[dumvar1,dumvar2,dumvar3,dumvar4]=iswap(1,ir(sub2ind(size(ir),max(imin,1)):end),1,ir(sub2ind(size(ir),max(kk,1)):end),1); dumvar2i=find((ir(sub2ind(size(ir),max(imin,1)):end))~=(dumvar2));dumvar4i=find((ir(sub2ind(size(ir),max(kk,1)):end))~=(dumvar4)); ir(imin-1+dumvar2i)=dumvar2(dumvar2i); ir(kk-1+dumvar4i)=dumvar4(dumvar4i);
mda_orig=mda; [n,dumvar2,mda,dumvar4,dumvar5]=dswap(n,a(sub2ind(size(a),imin,1):end),mda,a(sub2ind(size(a),kk,1):end),mda); mda(dumvar5~=mda_orig)=dumvar5(dumvar5~=mda_orig); a(sub2ind(size(a),imin,1):end)=dumvar2; a(sub2ind(size(a),kk,1):end)=dumvar4;
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,eb(sub2ind(size(eb),max(imin,1)):end),1,eb(sub2ind(size(eb),max(kk,1)):end),1); dumvar2i=find((eb(sub2ind(size(eb),max(imin,1)):end))~=(dumvar2));dumvar4i=find((eb(sub2ind(size(eb),max(kk,1)):end))~=(dumvar4)); eb(imin-1+dumvar2i)=dumvar2(dumvar2i); eb(kk-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,ub(sub2ind(size(ub),max(imin,1)):end),1,ub(sub2ind(size(ub),max(kk,1)):end),1); dumvar2i=find((ub(sub2ind(size(ub),max(imin,1)):end))~=(dumvar2));dumvar4i=find((ub(sub2ind(size(ub),max(kk,1)):end))~=(dumvar4)); ub(imin-1+dumvar2i)=dumvar2(dumvar2i); ub(kk-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,db(sub2ind(size(db),max(imin,1)):end),1,db(sub2ind(size(db),max(kk,1)):end),1); dumvar2i=find((db(sub2ind(size(db),max(imin,1)):end))~=(dumvar2));dumvar4i=find((db(sub2ind(size(db),max(kk,1)):end))~=(dumvar4)); db(imin-1+dumvar2i)=dumvar2(dumvar2i); db(kk-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,w(sub2ind(size(w),max(imin,1)):end),1,w(sub2ind(size(w),max(kk,1)):end),1); dumvar2i=find((w(sub2ind(size(w),max(imin,1)):end))~=(dumvar2));dumvar4i=find((w(sub2ind(size(w),max(kk,1)):end))~=(dumvar4)); w(imin-1+dumvar2i)=dumvar2(dumvar2i); w(kk-1+dumvar4i)=dumvar4(dumvar4i);
[dumvar1,dumvar2,dumvar3,dumvar4]=dswap(1,h(sub2ind(size(h),max(imin,1)):end),1,h(sub2ind(size(h),max(kk,1)):end),1); dumvar2i=find((h(sub2ind(size(h),max(imin,1)):end))~=(dumvar2));dumvar4i=find((h(sub2ind(size(h),max(kk,1)):end))~=(dumvar4)); h(imin-1+dumvar2i)=dumvar2(dumvar2i); h(kk-1+dumvar4i)=dumvar4(dumvar4i);
end;
if( j>krank )
gt2=1;
break;
end;
end;
if(gt2==1)
break;
end;
if( a(j,j)~=0.0d0 )
tn = (abs(tn).*sign(a(j,j)));
end;
[nn,dumvar2,a(sub2ind(size(a),j,j):end),mda]=dscal(nn,1.0d0./tn,a(sub2ind(size(a),j,j):end),mda);
a(j,j) = a(j,j) + 1.0d0;
if( j~=m )
for i = jp1 : m;
bb = -ddot(nn,a(sub2ind(size(a),j,j):end),mda,a(sub2ind(size(a),i,j):end),mda)./a(j,j);
mda_orig=mda; [nn,bb,dumvar3,mda,dumvar5,dumvar6]=daxpy(nn,bb,a(sub2ind(size(a),j,j):end),mda,a(sub2ind(size(a),i,j):end),mda); mda(dumvar6~=mda_orig)=dumvar6(dumvar6~=mda_orig); a(sub2ind(size(a),j,j):end)=dumvar3; a(sub2ind(size(a),i,j):end)=dumvar5;
if( i>np )
if( h(i)~=0.0d0 )
tt = 1.0d0 -(abs(a(i,j))./h(i)).^2;
tt = max(tt,0.0d0);
t = tt;
tt = 1.0d0 + .05d0.*tt.*(h(i)./w(i)).^2;
if( tt==1.0d0 )
[h(i) ,dumvar2,a(sub2ind(size(a),i,j+1):end),mda]=dnrm2(n-j,a(sub2ind(size(a),i,j+1):end),mda);
w(i) = h(i);
else;
h(i) = h(i).*sqrt(t);
end;
end;
end;
end; i = fix(m+1);
end;
h(j) = a(j,j);
a(j,j) = -tn;
%
%
% UPDATE UB, DB
%
ub(j) = ub(j)./abs(a(j,j));
db(j) =((abs(eb(j)).*sign(db(j)))+db(j))./a(j,j);
if( j==krank )
gt2=1;
break;
end;
for i = jp1 : krank;
ub(i) = ub(i) + abs(a(i,j)).*ub(j);
db(i) = db(i) - a(i,j).*db(j);
end; i = fix(krank+1);
end;
if(gt2==0)
xermsg('SLATEC','DU11US','FIRST NP ROWS ARE LINEARLY DEPENDENT',8,0);
krank = fix(j - 1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
ub_shape=zeros(ub_shape);ub_shape(:)=ub(1:numel(ub_shape));ub=ub_shape;
db_shape=zeros(db_shape);db_shape(:)=db(1:numel(db_shape));db=db_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;
eb_shape=zeros(eb_shape);eb_shape(:)=eb(1:numel(eb_shape));eb=eb_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;
%
% E N D M A I N L O O P
%
%
% COMPUTE KSURE
%
km1 = fix(krank - 1);
for i = 1 : km1;
is = 0;
kmi = fix(krank - i);
for ii = 1 : kmi;
if( ub(ii)>ub(ii+1) )
is = 1;
temp = ub(ii);
ub(ii) = ub(ii+1);
ub(ii+1) = temp;
end;
end; ii = fix(kmi+1);
if( is==0 )
break;
end;
end;
ksure = 0;
summlv = 0.0d0;
for i = 1 : krank;
r2 = ub(i).*ub(i);
if( r2+summlv>=1.0d0 )
break;
end;
summlv = summlv + r2;
ksure = fix(ksure + 1);
end;
%
% IF SYSTEM IS OF REDUCED RANK AND MODE = 2
% COMPLETE THE DECOMPOSITION FOR SHORTEST LEAST SQUARES SOLUTION
%
if( krank~=m && mode>=2 )
mmk = fix(m - krank);
kp1 = fix(krank + 1);
i = fix(krank);
while( true );
tn = dnrm2(mmk,a(sub2ind(size(a),kp1,i):end),1)./a(i,i);
tn = a(i,i).*sqrt(1.0d0+tn.*tn);
[mmk,dumvar2,a(sub2ind(size(a),kp1,i):end)]=dscal(mmk,1.0d0./tn,a(sub2ind(size(a),kp1,i):end),1);
w(i) = a(i,i)./tn + 1.0d0;
a(i,i) = -tn;
if( i==1 )
break;
end;
im1 = fix(i - 1);
for ii = 1 : im1;
tt = -ddot(mmk,a(sub2ind(size(a),kp1,ii):end),1,a(sub2ind(size(a),kp1,i):end),1)./w(i);
tt = tt - a(i,ii);
[mmk,tt,dumvar3,dumvar4,dumvar5]=daxpy(mmk,tt,a(sub2ind(size(a),kp1,i):end),1,a(sub2ind(size(a),kp1,ii):end),1); a(sub2ind(size(a),kp1,i):end)=dumvar3; a(sub2ind(size(a),kp1,ii):end)=dumvar5;
a(i,ii) = a(i,ii) + tt.*w(i);
end; ii = fix(im1+1);
i = fix(i - 1);
end;
end;
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
ub_shape=zeros(ub_shape);ub_shape(:)=ub(1:numel(ub_shape));ub=ub_shape;
db_shape=zeros(db_shape);db_shape(:)=db(1:numel(db_shape));db=db_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;
eb_shape=zeros(eb_shape);eb_shape(:)=eb(1:numel(eb_shape));eb=eb_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 du11us
%DECK DU12LS
|
|
