Code covered by the BSD License

### Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

```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 .
%             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;
%
ind = -3;
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;
%
%
%        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
```