| [a,lda,n,v,itask,ind,work]=spoir(a,lda,n,v,itask,ind,work); |
function [a,lda,n,v,itask,ind,work]=spoir(a,lda,n,v,itask,ind,work);
%***BEGIN PROLOGUE SPOIR
%***PURPOSE Solve a positive definite symmetric system of linear
% equations. Iterative refinement is used to obtain an error
% estimate.
%***LIBRARY SLATEC
%***CATEGORY D2B1B
%***TYPE SINGLE PRECISION (SPOIR-S, CPOIR-C)
%***KEYWORDS HERMITIAN, LINEAR EQUATIONS, POSITIVE DEFINITE, SYMMETRIC
%***AUTHOR Voorhees, E. A., (LANL)
%***DESCRIPTION
%
% subroutine SPOIR solves a real positive definite symmetric
% NxN system of single precision linear equations using LINPACK
% subroutines SPOFA and SPOSL. One pass of iterative refine-
% ment is used only to obtain an estimate of the accuracy. That
% is, if A is an NxN real positive definite symmetric matrix
% and if X and B are real N-vectors, then SPOIR solves the
% equation
%
% A*X=B.
%
% The matrix A is first factored into upper and lower
% triangular matrices R and R-TRANSPOSE. These
% factors are used to calculate the solution, X.
% Then the residual vector is found and used
% to calculate an estimate of the relative error, IND.
% IND estimates the accuracy of the solution only when the
% input matrix and the right hand side are represented
% exactly in the computer and does not take into account
% any errors in the input data.
%
% If the equation A*X=B is to be solved for more than one vector
% B, the factoring of A does not need to be performed again and
% the option to only solve (ITASK .GT. 1) will be faster for
% the succeeding solutions. In this case, the contents of A,
% LDA, N, and WORK must not have been altered by the user
% following factorization (ITASK=1). IND will not be changed
% by SPOIR in this case.
%
% Argument Description ***
% A REAL(LDA,N)
% the doubly subscripted array with dimension (LDA,N)
% which contains the coefficient matrix. Only the
% upper triangle, including the diagonal, of the
% coefficient matrix need be entered. A is not
% altered by the routine.
% LDA INTEGER
% the leading dimension of the array A. LDA must be great-
% er than or equal to N. (Terminal error message IND=-1)
% N INTEGER
% the order of the matrix A. N must be greater than
% or equal to one. (Terminal error message IND=-2)
% V REAL(N)
% on entry, the singly subscripted array(vector) of di-
% mension N which contains the right hand side B of a
% system of simultaneous linear equations A*X=B.
% on return, V contains the solution vector, X .
% ITASK INTEGER
% If ITASK = 1, the matrix A is factored and then the
% linear equation is solved.
% If ITASK .GT. 1, the equation is solved using the existing
% factored matrix A (stored in WORK).
% If ITASK .LT. 1, then terminal terminal error IND=-3 is
% printed.
% IND INTEGER
% GT. 0 IND is a rough estimate of the number of digits
% of accuracy in the solution, X. IND=75 means
% that the solution vector X is zero.
% LT. 0 See error message corresponding to IND below.
% WORK REAL(N*(N+1))
% a singly subscripted array of dimension at least N*(N+1).
%
% Error Messages Printed ***
%
% IND=-1 terminal N is greater than LDA.
% IND=-2 terminal N is less than one.
% IND=-3 terminal ITASK is less than one.
% IND=-4 Terminal The matrix A is computationally singular
% or is not positive definite.
% A solution has not been computed.
% IND=-10 warning The solution has no apparent significance.
% The solution may be inaccurate or the matrix
% A may be poorly scaled.
%
% Note- The above terminal(*fatal*) error messages are
% designed to be handled by XERMSG in which
% LEVEL=1 (recoverable) and IFLAG=2 . LEVEL=0
% for warning error messages from XERMSG. Unless
% the user provides otherwise, an error message
% will be printed followed by an abort.
%
%***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
% Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED DSDOT, R1MACH, SASUM, SCOPY, SPOFA, SPOSL, XERMSG
%***REVISION HISTORY (YYMMDD)
% 800528 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)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
% 900510 Convert XERRWV calls to XERMSG calls. (RWC)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE SPOIR
%
persistent dnorm info j xern1 xern2 xnorm ;
if isempty(info), info=0; end;
if isempty(j), j=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
v_shape=size(v);v=reshape(v,1,[]);
work_shape=size(work);work=reshape([work(:).',zeros(1,ceil(numel(work)./prod([n])).*prod([n])-numel(work))],n,[]);
if isempty(xnorm), xnorm=0; end;
if isempty(dnorm), dnorm=0; end;
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern2), xern2=repmat(' ',1,8); end;
%***FIRST EXECUTABLE STATEMENT SPOIR
if( lda<n )
ind = -1;
xern1=sprintf(['%8i'], lda);
xern2=sprintf(['%8i'], n);
xermsg('SLATEC','SPOIR',['LDA = ',[xern1,[' IS LESS THAN N = ',xern2]]],-1,1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
if( n<=0 )
ind = -2;
xern1=sprintf(['%8i'], n);
xermsg('SLATEC','SPOIR',['N = ',[xern1,' IS LESS THAN 1']],-2,1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
if( itask<1 )
ind = -3;
xern1=sprintf(['%8i'], itask);
xermsg('SLATEC','SPOIR',['ITASK = ',[xern1,' IS LESS THAN 1']],-3,1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
if( itask==1 )
%
% MOVE MATRIX A TO WORK
%
for j = 1 : n;
[n,a(sub2ind(size(a),1,j):end),dumvar3,work(sub2ind(size(work),1,j):end)]=scopy(n,a(sub2ind(size(a),1,j):end),1,work(sub2ind(size(work),1,j):end),1);
end; j = fix(n+1);
%
% FACTOR MATRIX A INTO R
n_orig=n; [work,n,dumvar3,info]=spofa(work,n,n,info); n(dumvar3~=n_orig)=dumvar3(dumvar3~=n_orig);
%
% CHECK FOR SINGULAR OR NOT POS.DEF. MATRIX
if( info~=0 )
ind = -4;
xermsg('SLATEC','SPOIR','SINGULAR OR NOT POSITIVE DEFINITE - NO SOLUTION',-4,1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end;
%
% SOLVE AFTER FACTORING
% MOVE VECTOR B TO WORK
%
[n,v(sub2ind(size(v),max(1,1)):end),dumvar3,work(sub2ind(size(work),1,n+1):end)]=scopy(n,v(sub2ind(size(v),max(1,1)):end),1,work(sub2ind(size(work),1,n+1):end),1);
n_orig=n; [work,n,dumvar3,v]=sposl(work,n,n,v); n(dumvar3~=n_orig)=dumvar3(dumvar3~=n_orig);
%
% FORM NORM OF X0
%
[xnorm ,n,v(sub2ind(size(v),max(1,1)):end)]=sasum(n,v(sub2ind(size(v),max(1,1)):end),1);
if( xnorm==0.0 )
ind = 75;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
% COMPUTE RESIDUAL
%
for j = 1 : n;
work(j,n+1) = -work(j,n+1) + dsdot(j-1,a(sub2ind(size(a),1,j):end),1,v(sub2ind(size(v),max(1,1)):end),1)+ dsdot(n-j+1,a(sub2ind(size(a),j,j):end),lda,v(sub2ind(size(v),max(j,1)):end),1);
end; j = fix(n+1);
%
% SOLVE A*DELTA=R
%
n_orig=n; [work,n,dumvar3,dumvar4]=sposl(work,n,n,work(sub2ind(size(work),1,n+1):end)); n(dumvar3~=n_orig)=dumvar3(dumvar3~=n_orig); work(sub2ind(size(work),1,n+1):end)=dumvar4;
%
% FORM NORM OF DELTA
%
[dnorm ,n,work(sub2ind(size(work),1,n+1):end)]=sasum(n,work(sub2ind(size(work),1,n+1):end),1);
%
% COMPUTE IND (ESTIMATE OF NO. OF SIGNIFICANT DIGITS)
% AND CHECK FOR IND GREATER THAN ZERO
%
ind = fix(-log10(max(r1mach(4),dnorm./xnorm)));
if( ind<=0 )
ind = -10;
xermsg('SLATEC','SPOIR','SOLUTION MAY HAVE NO SIGNIFICANCE',-10,0);
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK SPOPT
|
|
