| [a,x,b,n,m,nrda,u,nrdu,iflag,mlso,irank,iscale,q,diag,kpivot,s,div,td,isflg,scales]=lssuds(a,x,b,n,m,nrda,u,nrdu,iflag,mlso,irank,iscale,q,diag,kpivot,s,div,td,isflg,scales); |
function [a,x,b,n,m,nrda,u,nrdu,iflag,mlso,irank,iscale,q,diag,kpivot,s,div,td,isflg,scales]=lssuds(a,x,b,n,m,nrda,u,nrdu,iflag,mlso,irank,iscale,q,diag,kpivot,s,div,td,isflg,scales);
persistent gam gamma i irp j jr k kp l maxmes mj nfat nfatal nmir nu res ss uro ;
if isempty(gam), gam=0; end;
if isempty(gamma), gamma=0; end;
if isempty(res), res=0; end;
if isempty(ss), ss=0; end;
if isempty(uro), uro=0; end;
if isempty(i), i=0; end;
if isempty(irp), irp=0; end;
if isempty(j), j=0; end;
if isempty(jr), jr=0; end;
if isempty(k), k=0; end;
if isempty(kp), kp=0; end;
if isempty(l), l=0; end;
if isempty(maxmes), maxmes=0; end;
if isempty(mj), mj=0; end;
if isempty(nfat), nfat=0; end;
if isempty(nfatal), nfatal=0; end;
if isempty(nmir), nmir=0; end;
if isempty(nu), nu=0; end;
%***BEGIN PROLOGUE LSSUDS
%***SUBSIDIARY
%***PURPOSE Subsidiary to BVSUP
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (LSSUDS-S, DLSSUD-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% LSSUDS solves the underdetermined system of equations A Z = B,
% where A is N by M and N .LE. M. In particular, if rank A equals
% IRA, a vector X and a matrix U are determined such that X is the
% UNIQUE solution of smallest length, satisfying A X = B, and the
% columns of U form an orthonormal basis for the null space of A,
% satisfying A U = 0 . Then all solutions Z are given by
% Z = X + C(1)*U(1) + ..... + C(M-IRA)*U(M-IRA)
% where U(J) represents the J-th column of U and the C(J) are
% arbitrary constants.
% If the system of equations are not compatible, only the least
% squares solution of minimal length is computed.
%
% *********************************************************************
% INPUT
% *********************************************************************
%
% A -- Contains the matrix of N equations in M unknowns, A remains
% unchanged, must be dimensioned NRDA by M.
% X -- Solution array of length at least M.
% B -- Given constant vector of length N, B remains unchanged.
% N -- Number of equations, N greater or equal to 1.
% M -- Number of unknowns, M greater or equal to N.
% NRDA -- Row dimension of A, NRDA greater or equal to N.
% U -- Matrix used for solution, must be dimensioned NRDU by
% (M - rank of A).
% (storage for U may be ignored when only the minimal length
% solution X is desired)
% NRDU -- Row dimension of U, NRDU greater or equal to M.
% (if only the minimal length solution is wanted,
% NRDU=0 is acceptable)
% IFLAG -- Status indicator
% =0 for the first call (and for each new problem defined by
% a new matrix A) when the matrix data is treated as exact
% =-K for the first call (and for each new problem defined by
% a new matrix A) when the matrix data is assumed to be
% accurate to about K digits.
% =1 for subsequent calls whenever the matrix A has already
% been decomposed (problems with new vectors B but
% same matrix A can be handled efficiently).
% MLSO -- =0 if only the minimal length solution is wanted.
% =1 if the complete solution is wanted, includes the
% linear space defined by the matrix U.
% IRANK -- Variable used for the rank of A, set by the code.
% ISCALE -- Scaling indicator
% =-1 if the matrix A is to be pre-scaled by
% columns when appropriate.
% If the scaling indicator is not equal to -1
% no scaling will be attempted.
% For most problems scaling will probably not be necessary.
% Q -- Matrix used for the transformation, must be dimensioned
% NRDA by M.
% DIAG,KPIVOT,S, -- Arrays of length at least N used for internal
% DIV,TD,SCALES storage (except for SCALES which is M).
% ISFLG -- Storage for an internal variable.
%
% *********************************************************************
% OUTPUT
% *********************************************************************
%
% IFLAG -- Status indicator
% =1 if solution was obtained.
% =2 if improper input is detected.
% =3 if rank of matrix is less than N.
% To continue, simply reset IFLAG=1 and call LSSUDS again.
% =4 if the system of equations appears to be inconsistent.
% However, the least squares solution of minimal length
% was obtained.
% X -- Minimal length least squares solution of A Z = B
% IRANK -- Numerically determined rank of A, must not be altered
% on succeeding calls with input values of IFLAG=1.
% U -- Matrix whose M-IRANK columns are mutually orthogonal unit
% vectors which span the null space of A. This is to be ignored
% when MLSO was set to zero or IFLAG=4 on output.
% Q -- Contains the strictly upper triangular part of the reduced
% matrix and transformation information.
% DIAG -- Contains the diagonal elements of the triangular reduced
% matrix.
% KPIVOT -- Contains the pivotal information. The row interchanges
% performed on the original matrix are recorded here.
% S -- Contains the solution of the lower triangular system.
% DIV,TD -- Contains transformation information for rank
% deficient problems.
% SCALES -- Contains the column scaling parameters.
%
% *********************************************************************
%
%***SEE ALSO BVSUP
%***REFERENCES H. A. Watts, Solving linear least squares problems
% using SODS/SUDS/CODS, Sandia Report SAND77-0683,
% Sandia Laboratories, 1977.
%***ROUTINES CALLED J4SAVE, OHTROL, ORTHOR, R1MACH, SDOT, XERMAX,
% XERMSG, XGETF, XSETF
%***REVISION HISTORY (YYMMDD)
% 750601 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)
% 900510 Fixed an error message. (RWC)
% 910408 Updated the AUTHOR and REFERENCES sections. (WRB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE LSSUDS
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nrda])).*prod([nrda])-numel(a))],nrda,[]);
x_shape=size(x);x=reshape(x,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
u_shape=size(u);u=reshape([u(:).',zeros(1,ceil(numel(u)./prod([nrdu])).*prod([nrdu])-numel(u))],nrdu,[]);
q_shape=size(q);q=reshape([q(:).',zeros(1,ceil(numel(q)./prod([nrda])).*prod([nrda])-numel(q))],nrda,[]);
diag_shape=size(diag);diag=reshape(diag,1,[]);
kpivot_shape=size(kpivot);kpivot=reshape(kpivot,1,[]);
s_shape=size(s);s=reshape(s,1,[]);
div_shape=size(div);div=reshape(div,1,[]);
td_shape=size(td);td=reshape(td,1,[]);
scales_shape=size(scales);scales=reshape(scales,1,[]);
%
% **********************************************************************
%
% MACHINE PRECISION (COMPUTER UNIT ROUNDOFF VALUE) IS DEFINED
% BY THE FUNCTION R1MACH.
%
%***FIRST EXECUTABLE STATEMENT LSSUDS
[uro ]=r1mach(4);
%
% **********************************************************************
%
if( n>=1 && m>=n && nrda>=n )
if( nrdu==0 || nrdu>=m )
if( iflag<=0 )
%
[nfatal]=xgetf(nfatal);
[maxmes ]=j4save(4,0,false);
isflg = -15;
if( iflag~=0 )
isflg = fix(iflag);
nfat = -1;
if( nfatal==0 )
nfat = 0;
end;
[nfat]=xsetf(nfat);
xermax(1);
end;
%
% COPY MATRIX A INTO MATRIX Q
%
for k = 1 : m;
for j = 1 : n;
q(j,k) = a(j,k);
end; j = fix(n+1);
end; k = fix(m+1);
%
% use ORTHOGONAL TRANSFORMATIONS TO REDUCE Q TO LOWER
% TRIANGULAR FORM
%
[q,n,m,nrda,iflag,irank,iscale,diag,kpivot,scales,div,td]=orthor(q,n,m,nrda,iflag,irank,iscale,diag,kpivot,scales,div,td);
%
[nfatal]=xsetf(nfatal);
[maxmes]=xermax(maxmes);
if( irank==n )
%
% STORE DIVISORS FOR THE TRIANGULAR SOLUTION
%
for k = 1 : n;
div(k) = diag(k);
end; k = fix(n+1);
else;
%
% FOR RANK DEFICIENT PROBLEMS USE ADDITIONAL ORTHOGONAL
% TRANSFORMATIONS TO FURTHER REDUCE Q
%
if( irank~=0 )
[q,n,nrda,diag,irank,div,td]=ohtrol(q,n,nrda,diag,irank,div,td);
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
kpivot_shape=zeros(kpivot_shape);kpivot_shape(:)=kpivot(1:numel(kpivot_shape));kpivot=kpivot_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
div_shape=zeros(div_shape);div_shape(:)=div(1:numel(div_shape));div=div_shape;
td_shape=zeros(td_shape);td_shape(:)=td(1:numel(td_shape));td=td_shape;
scales_shape=zeros(scales_shape);scales_shape(:)=scales(1:numel(scales_shape));scales=scales_shape;
return;
end;
elseif( iflag~=1 ) ;
iflag = 2;
xermsg('SLATEC','LSSUDS','INVALID INPUT PARAMETERS.',2,1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
kpivot_shape=zeros(kpivot_shape);kpivot_shape(:)=kpivot(1:numel(kpivot_shape));kpivot=kpivot_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
div_shape=zeros(div_shape);div_shape(:)=div(1:numel(div_shape));div=div_shape;
td_shape=zeros(td_shape);td_shape(:)=td(1:numel(td_shape));td=td_shape;
scales_shape=zeros(scales_shape);scales_shape(:)=scales(1:numel(scales_shape));scales=scales_shape;
return;
end;
%
%
if( irank>0 )
%
% COPY CONSTANT VECTOR INTO S AFTER FIRST INTERCHANGING
% THE ELEMENTS ACCORDING TO THE PIVOTAL SEQUENCE
%
for k = 1 : n;
kp = fix(kpivot(k));
x(k) = b(kp);
end; k = fix(n+1);
for k = 1 : n;
s(k) = x(k);
end; k = fix(n+1);
%
irp = fix(irank + 1);
nu = 1;
if( mlso==0 )
nu = 0;
end;
if( irank~=n )
%
% FOR RANK DEFICIENT PROBLEMS WE MUST APPLY THE
% ORTHOGONAL TRANSFORMATION TO S
% WE ALSO CHECK TO SEE IF THE SYSTEM APPEARS TO BE INCONSISTENT
%
nmir = fix(n - irank);
[ss ,n,dumvar2,dumvar4,dumvar2]=sdot(n,s(sub2ind(size(s),max(1,1)):end),1,s(sub2ind(size(s),max(1,1)):end),1); dumvar2i=find((s(sub2ind(size(s),max(1,1)):end))~=(dumvar2)); s(1-1+dumvar2i)=dumvar2(dumvar2i);
for l = 1 : irank;
k = fix(irp - l);
gam =((td(k).*s(k))+sdot(nmir,q(sub2ind(size(q),irp,k):end),1,s(sub2ind(size(s),max(irp,1)):end),1))./(td(k).*div(k));
s(k) = s(k) + gam.*td(k);
for j = irp : n;
s(j) = s(j) + gam.*q(j,k);
end; j = fix(n+1);
end; l = fix(irank+1);
[res ,nmir,dumvar2,dumvar4,dumvar2]=sdot(nmir,s(sub2ind(size(s),max(irp,1)):end),1,s(sub2ind(size(s),max(irp,1)):end),1); dumvar2i=find((s(sub2ind(size(s),max(irp,1)):end))~=(dumvar2)); s(irp-1+dumvar2i)=dumvar2(dumvar2i);
if( res>ss.*(10..*max(10..^isflg,10..*uro)).^2 )
%
% INCONSISTENT SYSTEM
iflag = 4;
nu = 0;
end;
end;
%
% APPLY FORWARD SUBSTITUTION TO SOLVE LOWER TRIANGULAR SYSTEM
%
s(1) = s(1)./div(1);
if( irank~=1 )
for k = 2 : irank;
s(k) =(s(k)-sdot(k-1,q(sub2ind(size(q),k,1):end),nrda,s(sub2ind(size(s),max(1,1)):end),1))./div(k);
end; k = fix(irank+1);
end;
%
% INITIALIZE X VECTOR AND THEN APPLY ORTHOGONAL TRANSFORMATION
%
for k = 1 : m;
x(k) = 0.;
if( k<=irank )
x(k) = s(k);
end;
end; k = fix(m+1);
%
for jr = 1 : irank;
j = fix(irp - jr);
mj = fix(m - j + 1);
gamma = sdot(mj,q(sub2ind(size(q),j,j):end),nrda,x(sub2ind(size(x),max(j,1)):end),1)./(diag(j).*q(j,j));
for k = j : m;
x(k) = x(k) + gamma.*q(j,k);
end; k = fix(m+1);
end; jr = fix(irank+1);
%
% RESCALE ANSWERS AS DICTATED
%
for k = 1 : m;
x(k) = x(k).*scales(k);
end; k = fix(m+1);
%
if((nu==0) ||(m==irank) )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
kpivot_shape=zeros(kpivot_shape);kpivot_shape(:)=kpivot(1:numel(kpivot_shape));kpivot=kpivot_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
div_shape=zeros(div_shape);div_shape(:)=div(1:numel(div_shape));div=div_shape;
td_shape=zeros(td_shape);td_shape(:)=td(1:numel(td_shape));td=td_shape;
scales_shape=zeros(scales_shape);scales_shape(:)=scales(1:numel(scales_shape));scales=scales_shape;
return;
end;
%
% INITIALIZE U MATRIX AND THEN APPLY ORTHOGONAL TRANSFORMATION
%
l = fix(m - irank);
for k = 1 : l;
for i = 1 : m;
u(i,k) = 0.;
if( i==irank+k )
u(i,k) = 1.;
end;
end; i = fix(m+1);
%
for jr = 1 : irank;
j = fix(irp - jr);
mj = fix(m - j + 1);
gamma = sdot(mj,q(sub2ind(size(q),j,j):end),nrda,u(sub2ind(size(u),j,k):end),1)./(diag(j).*q(j,j));
for i = j : m;
u(i,k) = u(i,k) + gamma.*q(j,i);
end; i = fix(m+1);
end; jr = fix(irank+1);
end; k = fix(l+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
kpivot_shape=zeros(kpivot_shape);kpivot_shape(:)=kpivot(1:numel(kpivot_shape));kpivot=kpivot_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
div_shape=zeros(div_shape);div_shape(:)=div(1:numel(div_shape));div=div_shape;
td_shape=zeros(td_shape);td_shape(:)=td(1:numel(td_shape));td=td_shape;
scales_shape=zeros(scales_shape);scales_shape(:)=scales(1:numel(scales_shape));scales=scales_shape;
return;
else;
%
% SPECIAL CASE FOR THE NULL MATRIX
for k = 1 : m;
x(k) = 0.;
if( mlso~=0 )
u(k,k) = 1.;
for j = 1 : m;
if( j~=k )
u(j,k) = 0.;
end;
end; j = fix(m+1);
end;
end; k = fix(m+1);
for k = 1 : n;
if( b(k)>0. )
iflag = 4;
end;
end; k = fix(n+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
kpivot_shape=zeros(kpivot_shape);kpivot_shape(:)=kpivot(1:numel(kpivot_shape));kpivot=kpivot_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
div_shape=zeros(div_shape);div_shape(:)=div(1:numel(div_shape));div=div_shape;
td_shape=zeros(td_shape);td_shape(:)=td(1:numel(td_shape));td=td_shape;
scales_shape=zeros(scales_shape);scales_shape(:)=scales(1:numel(scales_shape));scales=scales_shape;
return;
end;
end;
end;
%
% INVALID INPUT FOR LSSUDS
iflag = 2;
xermsg('SLATEC','LSSUDS','INVALID INPUT PARAMETERS.',2,1);
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
kpivot_shape=zeros(kpivot_shape);kpivot_shape(:)=kpivot(1:numel(kpivot_shape));kpivot=kpivot_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
div_shape=zeros(div_shape);div_shape(:)=div(1:numel(div_shape));div=div_shape;
td_shape=zeros(td_shape);td_shape(:)=td(1:numel(td_shape));td=td_shape;
scales_shape=zeros(scales_shape);scales_shape(:)=scales(1:numel(scales_shape));scales=scales_shape;
end %subroutine lssuds
%DECK MACON
|
|
