[ap,n,b]=sppsl(ap,n,b);

function [ap,n,b]=sppsl(ap,n,b);
%***BEGIN PROLOGUE  SPPSL
%***PURPOSE  Solve the real symmetric positive definite system using
%            the factors computed by SPPCO or SPPFA.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B1B
%***TYPE      SINGLE PRECISION (SPPSL-S, DPPSL-D, CPPSL-C)
%***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX, PACKED,
%             POSITIVE DEFINITE, SOLVE
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     SPPSL solves the real symmetric positive definite system
%     A * X = B
%     using the factors computed by SPPCO or SPPFA.
%
%     On Entry
%
%        AP      REAL (N*(N+1)/2)
%                the output from SPPCO or SPPFA.
%
%        N       INTEGER
%                the order of the matrix  A .
%
%        B       REAL(N)
%                the right hand side vector.
%
%     On Return
%
%        B       the solution vector  X .
%
%     Error Condition
%
%        A division by zero will occur if the input factor contains
%        a zero on the diagonal.  Technically, this indicates
%        singularity, but it is usually caused by improper subroutine
%        arguments.  It will not occur if the subroutines are called
%        correctly and  INFO .EQ. 0 .
%
%     To compute  INVERSE(A) * C  where  C  is a matrix
%     with  P  columns
%           CALL SPPCO(AP,N,RCOND,Z,INFO)
%           IF (RCOND is too small .OR. INFO .NE. 0) GO TO ...
%           DO 10 J = 1, P
%              CALL SPPSL(AP,N,C(1,J))
%        10 CONTINUE
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  SAXPY, SDOT
%***REVISION HISTORY  (YYMMDD)
%   780814  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  SPPSL
persistent k kb kk t ; 

ap_shape=size(ap);ap=reshape(ap,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
%
if isempty(t), t=0; end;
if isempty(k), k=0; end;
if isempty(kb), kb=0; end;
if isempty(kk), kk=0; end;
%***FIRST EXECUTABLE STATEMENT  SPPSL
kk = 0;
for k = 1 : n;
[t ,dumvar2,ap(sub2ind(size(ap),max(kk+1,1)):end),dumvar4,b(sub2ind(size(b),max(1,1)):end)]=sdot(k-1,ap(sub2ind(size(ap),max(kk+1,1)):end),1,b(sub2ind(size(b),max(1,1)):end),1);
kk = fix(kk + k);
b(k) =(b(k)-t)./ap(kk);
end; k = fix(n+1);
for kb = 1 : n;
k = fix(n + 1 - kb);
b(k) = b(k)./ap(kk);
kk = fix(kk - k);
t = -b(k);
[dumvar1,t,ap(sub2ind(size(ap),max(kk+1,1)):end),dumvar4,b(sub2ind(size(b),max(1,1)):end)]=saxpy(k-1,t,ap(sub2ind(size(ap),max(kk+1,1)):end),1,b(sub2ind(size(b),max(1,1)):end),1);
end; kb = fix(n+1);
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
end
%DECK SPSORT