Code covered by the BSD License

### Highlights fromslatec

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

[ap,n,kpvt,info]=sspfa(ap,n,kpvt,info);
```function [ap,n,kpvt,info]=sspfa(ap,n,kpvt,info);
%***BEGIN PROLOGUE  SSPFA
%***PURPOSE  Factor a real symmetric matrix stored in packed form by
%            elimination with symmetric pivoting.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B1A
%***TYPE      SINGLE PRECISION (SSPFA-S, DSPFA-D, CHPFA-C, CSPFA-C)
%***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION, PACKED,
%             SYMMETRIC
%***AUTHOR  Bunch, J., (UCSD)
%***DESCRIPTION
%
%     SSPFA factors a real symmetric matrix stored in
%     packed form by elimination with symmetric pivoting.
%
%     To solve  A*X = B , follow SSPFA by SSPSL.
%     To compute  INVERSE(A)*C , follow SSPFA by SSPSL.
%     To compute  DETERMINANT(A) , follow SSPFA by SSPDI.
%     To compute  INERTIA(A) , follow SSPFA by SSPDI.
%     To compute  INVERSE(A) , follow SSPFA by SSPDI.
%
%     On Entry
%
%        AP      REAL (N*(N+1)/2)
%                the packed form of a symmetric matrix  A .  The
%                columns of the upper triangle are stored sequentially
%                in a one-dimensional array of length  N*(N+1)/2 .
%                See comments below for details.
%
%        N       INTEGER
%                the order of the matrix  A .
%
%     Output
%
%        AP      a block diagonal matrix and the multipliers which
%                were used to obtain it stored in packed form.
%                The factorization can be written  A = U*D*TRANS(U)
%                where  U  is a product of permutation and unit
%                upper triangular matrices , TRANS(U) is the
%                transpose of  U , and  D  is block diagonal
%                with 1 by 1 and 2 by 2 blocks.
%
%        KPVT    INTEGER(N)
%                an integer vector of pivot indices.
%
%        INFO    INTEGER
%                = 0  normal value.
%                = K  if the K-th pivot block is singular.  This is
%                     not an error condition for this subroutine,
%                     but it does indicate that SSPSL or SSPDI may
%                     divide by zero if called.
%
%     Packed Storage
%
%          The following program segment will pack the upper
%          triangle of a symmetric matrix.
%
%                K = 0
%                DO 20 J = 1, N
%                   DO 10 I = 1, J
%                      K = K + 1
%                      AP(K)  = A(I,J)
%             10    CONTINUE
%             20 CONTINUE
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  ISAMAX, SAXPY, SSWAP
%***REVISION HISTORY  (YYMMDD)
%   780814  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891107  Modified routine equivalence list.  (WRB)
%   891107  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  SSPFA
persistent absakk ak akm1 alpha bk bkm1 colmax denom ij ik ikm1 im imax imaxp1 imim imj imk j jj jk jkm1 jmax jmim k kk km1 km1k km1km1 km2 kstep mulk mulkm1 rowmax swap t ;

kpvt_shape=size(kpvt);kpvt=reshape(kpvt,1,[]);
ap_shape=size(ap);ap=reshape(ap,1,[]);
%
if isempty(ak), ak=0; end;
if isempty(akm1), akm1=0; end;
if isempty(bk), bk=0; end;
if isempty(bkm1), bkm1=0; end;
if isempty(denom), denom=0; end;
if isempty(mulk), mulk=0; end;
if isempty(mulkm1), mulkm1=0; end;
if isempty(t), t=0; end;
if isempty(absakk), absakk=0; end;
if isempty(alpha), alpha=0; end;
if isempty(colmax), colmax=0; end;
if isempty(rowmax), rowmax=0; end;
if isempty(ij), ij=0; end;
if isempty(ik), ik=0; end;
if isempty(ikm1), ikm1=0; end;
if isempty(im), im=0; end;
if isempty(imax), imax=0; end;
if isempty(imaxp1), imaxp1=0; end;
if isempty(imim), imim=0; end;
if isempty(imj), imj=0; end;
if isempty(imk), imk=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(jk), jk=0; end;
if isempty(jkm1), jkm1=0; end;
if isempty(jmax), jmax=0; end;
if isempty(jmim), jmim=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(km1), km1=0; end;
if isempty(km1k), km1k=0; end;
if isempty(km1km1), km1km1=0; end;
if isempty(km2), km2=0; end;
if isempty(kstep), kstep=0; end;
if isempty(swap), swap=false; end;
%***FIRST EXECUTABLE STATEMENT  SSPFA
%
%     INITIALIZE
%
%     ALPHA IS USED IN CHOOSING PIVOT BLOCK SIZE.
%
alpha =(1.0e0+sqrt(17.0e0))./8.0e0;
%
info = 0;
%
%     MAIN LOOP ON K, WHICH GOES FROM N TO 1.
%
k = fix(n);
ik =fix(fix((n.*(n-1))./2));
%
%        LEAVE THE LOOP IF K=0 OR K=1.
%
while( k~=0 );
if( k<=1 )
kpvt(1) = 1;
if( ap(1)==0.0e0 )
info = 1;
end;
break;
else;
%
%        THIS SECTION OF CODE DETERMINES THE KIND OF
%        ELIMINATION TO BE PERFORMED.  WHEN IT IS COMPLETED,
%        KSTEP WILL BE SET TO THE SIZE OF THE PIVOT BLOCK, AND
%        SWAP WILL BE SET TO true IF AN INTERCHANGE IS
%        REQUIRED.
%
km1 = fix(k - 1);
kk = fix(ik + k);
absakk = abs(ap(kk));
%
%        DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
%        COLUMN K.
%
[imax ,dumvar2,ap(sub2ind(size(ap),max(ik+1,1)):end)]=isamax(k-1,ap(sub2ind(size(ap),max(ik+1,1)):end),1);
imk = fix(ik + imax);
colmax = abs(ap(imk));
if( absakk<alpha.*colmax )
%
%           DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
%           ROW IMAX.
%
rowmax = 0.0e0;
imaxp1 = fix(imax + 1);
im = fix(fix((imax.*(imax-1))./2));
imj = fix(im + 2.*imax);
for j = imaxp1 : k;
rowmax = max(rowmax,abs(ap(imj)));
imj = fix(imj + j);
end; j = fix(k+1);
if( imax~=1 )
[jmax ,dumvar2,ap(sub2ind(size(ap),max(im+1,1)):end)]=isamax(imax-1,ap(sub2ind(size(ap),max(im+1,1)):end),1);
jmim = fix(jmax + im);
rowmax = max(rowmax,abs(ap(jmim)));
end;
imim = fix(imax + im);
if( abs(ap(imim))>=alpha.*rowmax )
kstep = 1;
swap = true;
elseif( absakk<alpha.*colmax.*(colmax./rowmax) ) ;
kstep = 2;
swap = imax~=km1;
else;
kstep = 1;
swap = false;
end;
else;
kstep = 1;
swap = false;
end;
if( max(absakk,colmax)==0.0e0 )
%
%           COLUMN K IS ZERO.  SET INFO AND ITERATE THE LOOP.
%
kpvt(k) = fix(k);
info = fix(k);
elseif( kstep==2 ) ;
%
%           2 X 2 PIVOT BLOCK.
%
km1k = fix(ik + k - 1);
ikm1 = fix(ik -(k-1));
if( swap )
%
%              PERFORM AN INTERCHANGE.
%
[imax,dumvar2,dumvar3,dumvar4]=sswap(imax,ap(sub2ind(size(ap),max(im+1,1)):end),1,ap(sub2ind(size(ap),max(ikm1+1,1)):end),1);   dumvar2i=find((ap(sub2ind(size(ap),max(im+1,1)):end))~=(dumvar2));dumvar4i=find((ap(sub2ind(size(ap),max(ikm1+1,1)):end))~=(dumvar4));   ap(im+1-1+dumvar2i)=dumvar2(dumvar2i); ap(ikm1+1-1+dumvar4i)=dumvar4(dumvar4i);
imj = fix(ikm1 + imax);
for jj = imax : km1;
j = fix(km1 + imax - jj);
jkm1 = fix(ikm1 + j);
t = ap(jkm1);
ap(jkm1) = ap(imj);
ap(imj) = t;
imj = fix(imj -(j-1));
end; jj = fix(km1+1);
t = ap(km1k);
ap(km1k) = ap(imk);
ap(imk) = t;
end;
%
%           PERFORM THE ELIMINATION.
%
km2 = fix(k - 2);
if( km2~=0 )
ak = ap(kk)./ap(km1k);
km1km1 = fix(ikm1 + k - 1);
akm1 = ap(km1km1)./ap(km1k);
denom = 1.0e0 - ak.*akm1;
ij = fix(ik -(k-1) -(k-2));
for jj = 1 : km2;
j = fix(km1 - jj);
jk = fix(ik + j);
bk = ap(jk)./ap(km1k);
jkm1 = fix(ikm1 + j);
bkm1 = ap(jkm1)./ap(km1k);
mulk =(akm1.*bk-bkm1)./denom;
mulkm1 =(ak.*bkm1-bk)./denom;
t = mulk;
[j,t,dumvar3,dumvar4,dumvar5]=saxpy(j,t,ap(sub2ind(size(ap),max(ik+1,1)):end),1,ap(sub2ind(size(ap),max(ij+1,1)):end),1);   dumvar3i=find((ap(sub2ind(size(ap),max(ik+1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(ij+1,1)):end))~=(dumvar5));   ap(ik+1-1+dumvar3i)=dumvar3(dumvar3i); ap(ij+1-1+dumvar5i)=dumvar5(dumvar5i);
t = mulkm1;
[j,t,dumvar3,dumvar4,dumvar5]=saxpy(j,t,ap(sub2ind(size(ap),max(ikm1+1,1)):end),1,ap(sub2ind(size(ap),max(ij+1,1)):end),1);   dumvar3i=find((ap(sub2ind(size(ap),max(ikm1+1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(ij+1,1)):end))~=(dumvar5));   ap(ikm1+1-1+dumvar3i)=dumvar3(dumvar3i); ap(ij+1-1+dumvar5i)=dumvar5(dumvar5i);
ap(jk) = mulk;
ap(jkm1) = mulkm1;
ij = fix(ij -(j-1));
end; jj = fix(km2+1);
end;
%
%           SET THE PIVOT ARRAY.
%
kpvt(k) = fix(1 - k);
if( swap )
kpvt(k) = fix(-imax);
end;
kpvt(k-1) = fix(kpvt(k));
else;
%
%           1 X 1 PIVOT BLOCK.
%
if( swap )
%
%              PERFORM AN INTERCHANGE.
%
[imax,dumvar2,dumvar3,dumvar4]=sswap(imax,ap(sub2ind(size(ap),max(im+1,1)):end),1,ap(sub2ind(size(ap),max(ik+1,1)):end),1);   dumvar2i=find((ap(sub2ind(size(ap),max(im+1,1)):end))~=(dumvar2));dumvar4i=find((ap(sub2ind(size(ap),max(ik+1,1)):end))~=(dumvar4));   ap(im+1-1+dumvar2i)=dumvar2(dumvar2i); ap(ik+1-1+dumvar4i)=dumvar4(dumvar4i);
imj = fix(ik + imax);
for jj = imax : k;
j = fix(k + imax - jj);
jk = fix(ik + j);
t = ap(jk);
ap(jk) = ap(imj);
ap(imj) = t;
imj = fix(imj -(j-1));
end; jj = fix(k+1);
end;
%
%           PERFORM THE ELIMINATION.
%
ij = fix(ik -(k-1));
for jj = 1 : km1;
j = fix(k - jj);
jk = fix(ik + j);
mulk = -ap(jk)./ap(kk);
t = mulk;
[j,t,dumvar3,dumvar4,dumvar5]=saxpy(j,t,ap(sub2ind(size(ap),max(ik+1,1)):end),1,ap(sub2ind(size(ap),max(ij+1,1)):end),1);   dumvar3i=find((ap(sub2ind(size(ap),max(ik+1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(ij+1,1)):end))~=(dumvar5));   ap(ik+1-1+dumvar3i)=dumvar3(dumvar3i); ap(ij+1-1+dumvar5i)=dumvar5(dumvar5i);
ap(jk) = mulk;
ij = fix(ij -(j-1));
end; jj = fix(km1+1);
%
%           SET THE PIVOT ARRAY.
%
kpvt(k) = fix(k);
if( swap )
kpvt(k) = fix(imax);
end;
end;
ik = fix(ik -(k-1));
if( kstep==2 )
ik = fix(ik -(k-2));
end;
k = fix(k - kstep);
end;
end;
kpvt_shape=zeros(kpvt_shape);kpvt_shape(:)=kpvt(1:numel(kpvt_shape));kpvt=kpvt_shape;
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
end %subroutine sspfa
%DECK SSPMV
```