| [a,mda,m,n,b,mdb,nb,re,ae,key,mode,np,krank,ksure,rnorm,w,lw,iwork,liw,info]=ulsia(a,mda,m,n,b,mdb,nb,re,ae,key,mode,np,krank,ksure,rnorm,w,lw,iwork,liw,info); |
function [a,mda,m,n,b,mdb,nb,re,ae,key,mode,np,krank,ksure,rnorm,w,lw,iwork,liw,info]=ulsia(a,mda,m,n,b,mdb,nb,re,ae,key,mode,np,krank,ksure,rnorm,w,lw,iwork,liw,info);
persistent eps i it m1 m2 m3 m4 m5 ;
if isempty(eps), eps=0; end;
if isempty(i), i=0; end;
if isempty(it), it=0; end;
if isempty(m1), m1=0; end;
if isempty(m2), m2=0; end;
if isempty(m3), m3=0; end;
if isempty(m4), m4=0; end;
if isempty(m5), m5=0; end;
%***BEGIN PROLOGUE ULSIA
%***PURPOSE Solve an underdetermined linear system of equations by
% performing an LQ factorization of the matrix using
% Householder transformations. Emphasis is put on detecting
% possible rank deficiency.
%***LIBRARY SLATEC
%***CATEGORY D9
%***TYPE SINGLE PRECISION (ULSIA-S, DULSIA-D)
%***KEYWORDS LINEAR LEAST SQUARES, LQ FACTORIZATION,
% UNDERDETERMINED LINEAR SYSTEM
%***AUTHOR Manteuffel, T. A., (LANL)
%***DESCRIPTION
%
% ULSIA computes the minimal length solution(s) to the problem AX=B
% where A is an M by N matrix with M.LE.N and B is the M by NB
% matrix of right hand sides. User input bounds on the uncertainty
% in the elements of A are used to detect numerical rank deficiency.
% The algorithm employs a row and column pivot strategy to
% minimize the growth of uncertainty and round-off errors.
%
% ULSIA requires (MDA+1)*N + (MDB+1)*NB + 6*M dimensioned space
%
% ******************************************************************
% * *
% * WARNING - All input arrays are changed on exit. *
% * *
% ******************************************************************
%
% Input..
%
% A() Linear coefficient matrix of AX=B, with MDA the
% MDA,M,N actual first dimension of A in the calling program.
% M is the row dimension (no. of EQUATIONS of the
% problem) and N the col dimension (no. of UNKNOWNS).
% Must have MDA.GE.M and M.LE.N.
%
% B() Right hand side(s), with MDB the actual first
% MDB,NB dimension of B in the calling program. NB is the
% number of M by 1 right hand sides. Since the
% solution is returned in B, must have MDB.GE.N. If
% NB = 0, B is never accessed.
%
% ******************************************************************
% * *
% * Note - Use of RE and AE are what make this *
% * code significantly different from *
% * other linear least squares solvers. *
% * However, the inexperienced user is *
% * advised to set RE=0.,AE=0.,KEY=0. *
% * *
% ******************************************************************
%
% RE(),AE(),KEY
% RE() RE() is a vector of length N such that RE(I) is
% the maximum relative uncertainty in row I of
% the matrix A. The values of RE() must be between
% 0 and 1. A minimum of 10*machine precision will
% be enforced.
%
% AE() AE() is a vector of length N such that AE(I) is
% the maximum absolute uncertainty in row I of
% the matrix A. The values of AE() must be greater
% than or equal to 0.
%
% KEY For ease of use, RE and AE may be input as either
% vectors or scalars. If a scalar is input, the algo-
% rithm will use that value for each column of A.
% The parameter KEY indicates whether scalars or
% vectors are being input.
% KEY=0 RE scalar AE scalar
% KEY=1 RE vector AE scalar
% KEY=2 RE scalar AE vector
% KEY=3 RE vector AE vector
%
%
% MODE The integer MODE indicates how the routine
% is to react if rank deficiency is detected.
% If MODE = 0 return immediately, no solution
% 1 compute truncated solution
% 2 compute minimal length least squares sol
% The inexperienced user is advised to set MODE=0
%
% NP The first NP rows of A will not be interchanged
% with other rows even though the pivot strategy
% would suggest otherwise.
% The inexperienced user is advised to set NP=0.
%
% WORK() A real work array dimensioned 5*M. However, if
% RE or AE have been specified as vectors, dimension
% WORK 4*M. If both RE and AE have been specified
% as vectors, dimension WORK 3*M.
%
% LW Actual dimension of WORK
%
% IWORK() Integer work array dimensioned at least N+M.
%
% LIW Actual dimension of IWORK.
%
%
% INFO Is a flag which provides for the efficient
% solution of subsequent problems involving the
% same A but different B.
% If INFO = 0 original call
% INFO = 1 subsequent calls
% On subsequent calls, the user must supply A, KRANK,
% LW, IWORK, LIW, and the first 2*M locations of WORK
% as output by the original call to ULSIA. MODE must
% be equal to the value of MODE in the original call.
% If MODE.LT.2, only the first N locations of WORK
% are accessed. AE, RE, KEY, and NP are not accessed.
%
%
%
%
% Output..
%
% A() Contains the lower triangular part of the reduced
% matrix and the transformation information. It togeth
% with the first M elements of WORK (see below)
% completely specify the LQ factorization of A.
%
% B() Contains the N by NB solution matrix for X.
%
% KRANK,KSURE The numerical rank of A, based upon the relative
% and absolute bounds on uncertainty, is bounded
% above by KRANK and below by KSURE. The algorithm
% returns a solution based on KRANK. KSURE provides
% an indication of the precision of the rank.
%
% RNORM() Contains the Euclidean length of the NB residual
% vectors B(I)-AX(I), I=1,NB. If the matrix A is of
% full rank, then RNORM=0.0.
%
% WORK() The first M locations of WORK contain values
% necessary to reproduce the Householder
% transformation.
%
% IWORK() The first N locations contain the order in
% which the columns of A were used. The next
% M locations contain the order in which the
% rows of A were used.
%
% INFO Flag to indicate status of computation on completion
% -1 Parameter error(s)
% 0 - Rank deficient, no solution
% 1 - Rank deficient, truncated solution
% 2 - Rank deficient, minimal length least squares sol
% 3 - Numerical rank 0, zero solution
% 4 - Rank .LT. NP
% 5 - Full rank
%
%***REFERENCES T. Manteuffel, An interval analysis approach to rank
% determination in linear least squares problems,
% Report SAND80-0655, Sandia Laboratories, June 1980.
%***ROUTINES CALLED R1MACH, U11US, U12US, XERMSG
%***REVISION HISTORY (YYMMDD)
% 810801 DATE WRITTEN
% 890831 Modified array declarations. (WRB)
% 891009 Removed unreferenced variable. (WRB)
% 891009 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 Fixed an error message. (RWC)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE ULSIA
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,[]);
re_shape=size(re);re=reshape(re,1,[]);
ae_shape=size(ae);ae=reshape(ae,1,[]);
rnorm_shape=size(rnorm);rnorm=reshape(rnorm,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
%
%***FIRST EXECUTABLE STATEMENT ULSIA
if( info<0 || info>1 )
xermsg('SLATEC','ULSIA','INFO OUT OF RANGE',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
it = fix(info);
info = -1;
if( nb==0 && it==1 )
%
% ERROR MESSAGES
%
xermsg('SLATEC','ULSIA','SOLUTION ONLY (INFO=1) BUT NO RIGHT HAND SIDE (NB=0)',1,0);
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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif( m<1 ) ;
xermsg('SLATEC','ULSIA','M.LT.1',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
if( n<1 )
xermsg('SLATEC','ULSIA','N.LT.1',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
if( n<m )
xermsg('SLATEC','ULSIA','N.LT.M',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
if( mda<m )
xermsg('SLATEC','ULSIA','MDA.LT.M',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
if( liw<m+n )
xermsg('SLATEC','ULSIA','LIW.LT.M+N',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
while (1);
if( mode<0 || mode>3 )
xermsg('SLATEC','ULSIA','MODE OUT OF RANGE',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
if( nb~=0 )
if( nb<0 )
xermsg('SLATEC','ULSIA','NB.LT.0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif( mdb<n ) ;
xermsg('SLATEC','ULSIA','MDB.LT.N',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif( it~=0 ) ;
break;
end;
end;
if( key<0 || key>3 )
xermsg('SLATEC','ULSIA','KEY OUT OF RANGE',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
if( key==0 && lw<5.*m )
xermsg('SLATEC','ULSIA','INSUFFICIENT WORK SPACE',8,1);
info = -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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( key==1 && lw<4.*m )
xermsg('SLATEC','ULSIA','INSUFFICIENT WORK SPACE',8,1);
info = -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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( key==2 && lw<4.*m )
xermsg('SLATEC','ULSIA','INSUFFICIENT WORK SPACE',8,1);
info = -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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( key==3 && lw<3.*m )
xermsg('SLATEC','ULSIA','INSUFFICIENT WORK SPACE',8,1);
info = -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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( np<0 || np>m )
xermsg('SLATEC','ULSIA','NP OUT OF RANGE',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
%
eps = 10..*r1mach(4);
m1 = 1;
m2 = fix(m1 + m);
m3 = fix(m2 + m);
m4 = fix(m3 + m);
m5 = fix(m4 + m);
%
if( key==1 )
%
if( ae(1)<0.0 )
xermsg('SLATEC','ULSIA','AE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
for i = 1 : m;
if( re(i)<0.0 )
xermsg('SLATEC','ULSIA','RE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( re(i)>1.0 )
xermsg('SLATEC','ULSIA','RE(I) .GT. 1',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( re(i)<eps )
re(i) = eps;
end;
w(m4-1+i) = ae(1);
end; i = fix(m+1);
[a,mda,m,n,re,dumvar6,mode,np,krank,ksure,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15]=u11us(a,mda,m,n,re,w(sub2ind(size(w),max(m4,1)):end),mode,np,krank,ksure,w(sub2ind(size(w),max(m1,1)):end),w(sub2ind(size(w),max(m2,1)):end),w(sub2ind(size(w),max(m3,1)):end),iwork(sub2ind(size(iwork),max(m1,1)):end),iwork(sub2ind(size(iwork),max(m2,1)):end)); dumvar6i=find((w(sub2ind(size(w),max(m4,1)):end))~=(dumvar6));dumvar11i=find((w(sub2ind(size(w),max(m1,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(m2,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(m3,1)):end))~=(dumvar13));dumvar14i=find((iwork(sub2ind(size(iwork),max(m1,1)):end))~=(dumvar14));dumvar15i=find((iwork(sub2ind(size(iwork),max(m2,1)):end))~=(dumvar15)); w(m4-1+dumvar6i)=dumvar6(dumvar6i); w(m1-1+dumvar11i)=dumvar11(dumvar11i); w(m2-1+dumvar12i)=dumvar12(dumvar12i); w(m3-1+dumvar13i)=dumvar13(dumvar13i); iwork(m1-1+dumvar14i)=dumvar14(dumvar14i); iwork(m2-1+dumvar15i)=dumvar15(dumvar15i);
elseif( key==2 ) ;
%
if( re(1)<0.0 )
xermsg('SLATEC','ULSIA','RE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( re(1)>1.0 )
xermsg('SLATEC','ULSIA','RE(I) .GT. 1',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( re(1)<eps )
re(1) = eps;
end;
for i = 1 : m;
w(m4-1+i) = re(1);
if( ae(i)<0.0 )
xermsg('SLATEC','ULSIA','AE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
end; i = fix(m+1);
[a,mda,m,n,dumvar5,ae,mode,np,krank,ksure,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15]=u11us(a,mda,m,n,w(sub2ind(size(w),max(m4,1)):end),ae,mode,np,krank,ksure,w(sub2ind(size(w),max(m1,1)):end),w(sub2ind(size(w),max(m2,1)):end),w(sub2ind(size(w),max(m3,1)):end),iwork(sub2ind(size(iwork),max(m1,1)):end),iwork(sub2ind(size(iwork),max(m2,1)):end)); dumvar5i=find((w(sub2ind(size(w),max(m4,1)):end))~=(dumvar5));dumvar11i=find((w(sub2ind(size(w),max(m1,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(m2,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(m3,1)):end))~=(dumvar13));dumvar14i=find((iwork(sub2ind(size(iwork),max(m1,1)):end))~=(dumvar14));dumvar15i=find((iwork(sub2ind(size(iwork),max(m2,1)):end))~=(dumvar15)); w(m4-1+dumvar5i)=dumvar5(dumvar5i); w(m1-1+dumvar11i)=dumvar11(dumvar11i); w(m2-1+dumvar12i)=dumvar12(dumvar12i); w(m3-1+dumvar13i)=dumvar13(dumvar13i); iwork(m1-1+dumvar14i)=dumvar14(dumvar14i); iwork(m2-1+dumvar15i)=dumvar15(dumvar15i);
elseif( key==3 ) ;
%
for i = 1 : m;
if( re(i)<0.0 )
xermsg('SLATEC','ULSIA','RE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( re(i)>1.0 )
xermsg('SLATEC','ULSIA','RE(I) .GT. 1',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( re(i)<eps )
re(i) = eps;
end;
if( ae(i)<0.0 )
xermsg('SLATEC','ULSIA','AE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
end; i = fix(m+1);
[a,mda,m,n,re,ae,mode,np,krank,ksure,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15]=u11us(a,mda,m,n,re,ae,mode,np,krank,ksure,w(sub2ind(size(w),max(m1,1)):end),w(sub2ind(size(w),max(m2,1)):end),w(sub2ind(size(w),max(m3,1)):end),iwork(sub2ind(size(iwork),max(m1,1)):end),iwork(sub2ind(size(iwork),max(m2,1)):end)); dumvar11i=find((w(sub2ind(size(w),max(m1,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(m2,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(m3,1)):end))~=(dumvar13));dumvar14i=find((iwork(sub2ind(size(iwork),max(m1,1)):end))~=(dumvar14));dumvar15i=find((iwork(sub2ind(size(iwork),max(m2,1)):end))~=(dumvar15)); w(m1-1+dumvar11i)=dumvar11(dumvar11i); w(m2-1+dumvar12i)=dumvar12(dumvar12i); w(m3-1+dumvar13i)=dumvar13(dumvar13i); iwork(m1-1+dumvar14i)=dumvar14(dumvar14i); iwork(m2-1+dumvar15i)=dumvar15(dumvar15i);
else;
%
if( re(1)<0.0 )
xermsg('SLATEC','ULSIA','RE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( re(1)>1.0 )
xermsg('SLATEC','ULSIA','RE(I) .GT. 1',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( re(1)<eps )
re(1) = eps;
end;
if( ae(1)<0.0 )
xermsg('SLATEC','ULSIA','AE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
for i = 1 : m;
w(m4-1+i) = re(1);
w(m5-1+i) = ae(1);
end; i = fix(m+1);
[a,mda,m,n,dumvar5,dumvar6,mode,np,krank,ksure,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15]=u11us(a,mda,m,n,w(sub2ind(size(w),max(m4,1)):end),w(sub2ind(size(w),max(m5,1)):end),mode,np,krank,ksure,w(sub2ind(size(w),max(m1,1)):end),w(sub2ind(size(w),max(m2,1)):end),w(sub2ind(size(w),max(m3,1)):end),iwork(sub2ind(size(iwork),max(m1,1)):end),iwork(sub2ind(size(iwork),max(m2,1)):end)); dumvar5i=find((w(sub2ind(size(w),max(m4,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(m5,1)):end))~=(dumvar6));dumvar11i=find((w(sub2ind(size(w),max(m1,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(m2,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(m3,1)):end))~=(dumvar13));dumvar14i=find((iwork(sub2ind(size(iwork),max(m1,1)):end))~=(dumvar14));dumvar15i=find((iwork(sub2ind(size(iwork),max(m2,1)):end))~=(dumvar15)); w(m4-1+dumvar5i)=dumvar5(dumvar5i); w(m5-1+dumvar6i)=dumvar6(dumvar6i); w(m1-1+dumvar11i)=dumvar11(dumvar11i); w(m2-1+dumvar12i)=dumvar12(dumvar12i); w(m3-1+dumvar13i)=dumvar13(dumvar13i); iwork(m1-1+dumvar14i)=dumvar14(dumvar14i); iwork(m2-1+dumvar15i)=dumvar15(dumvar15i);
end;
end;
end;
end;
break;
end;
%
% DETERMINE INFO
%
if( krank==m )
info = 5;
elseif( krank==0 ) ;
info = 3;
elseif( krank>=np ) ;
info = fix(mode);
if( mode==0 )
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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
else;
info = 4;
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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( nb==0 )
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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
%
%
% SOLUTION PHASE
%
m1 = 1;
m2 = fix(m1 + m);
m3 = fix(m2 + m);
if( info==2 )
%
if( lw>=m3-1 )
[a,mda,m,n,b,mdb,nb,mode,krank,rnorm,dumvar11,dumvar12,dumvar13,dumvar14]=u12us(a,mda,m,n,b,mdb,nb,mode,krank,rnorm,w(sub2ind(size(w),max(m1,1)):end),w(sub2ind(size(w),max(m2,1)):end),iwork(sub2ind(size(iwork),max(m1,1)):end),iwork(sub2ind(size(iwork),max(m2,1)):end)); dumvar11i=find((w(sub2ind(size(w),max(m1,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(m2,1)):end))~=(dumvar12));dumvar13i=find((iwork(sub2ind(size(iwork),max(m1,1)):end))~=(dumvar13));dumvar14i=find((iwork(sub2ind(size(iwork),max(m2,1)):end))~=(dumvar14)); w(m1-1+dumvar11i)=dumvar11(dumvar11i); w(m2-1+dumvar12i)=dumvar12(dumvar12i); iwork(m1-1+dumvar13i)=dumvar13(dumvar13i); iwork(m2-1+dumvar14i)=dumvar14(dumvar14i);
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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
elseif( lw>=m2-1 ) ;
[a,mda,m,n,b,mdb,nb,mode,krank,rnorm,dumvar11,dumvar11,dumvar13,dumvar14]=u12us(a,mda,m,n,b,mdb,nb,mode,krank,rnorm,w(sub2ind(size(w),max(m1,1)):end),w(sub2ind(size(w),max(m1,1)):end),iwork(sub2ind(size(iwork),max(m1,1)):end),iwork(sub2ind(size(iwork),max(m2,1)):end)); dumvar11i=find((w(sub2ind(size(w),max(m1,1)):end))~=(dumvar11));dumvar13i=find((iwork(sub2ind(size(iwork),max(m1,1)):end))~=(dumvar13));dumvar14i=find((iwork(sub2ind(size(iwork),max(m2,1)):end))~=(dumvar14)); w(m1-1+dumvar11i)=dumvar11(dumvar11i); iwork(m1-1+dumvar13i)=dumvar13(dumvar13i); iwork(m2-1+dumvar14i)=dumvar14(dumvar14i);
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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
end;
xermsg('SLATEC','ULSIA','INSUFFICIENT WORK SPACE',8,1);
info = -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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
xermsg('SLATEC','ULSIA','RE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
xermsg('SLATEC','ULSIA','RE(I) .GT. 1',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
xermsg('SLATEC','ULSIA','AE(I) .LT. 0',2,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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
end;
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;
re_shape=zeros(re_shape);re_shape(:)=re(1:numel(re_shape));re=re_shape;
ae_shape=zeros(ae_shape);ae_shape(:)=ae(1:numel(ae_shape));ae=ae_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
end %subroutine ulsia
%DECK USRMAT
|
|
