[a,lda,n,info]=spofa(a,lda,n,info);

function [a,lda,n,info]=spofa(a,lda,n,info);
%***BEGIN PROLOGUE  SPOFA
%***PURPOSE  Factor a real symmetric positive definite matrix.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B1B
%***TYPE      SINGLE PRECISION (SPOFA-S, DPOFA-D, CPOFA-C)
%***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION,
%             POSITIVE DEFINITE
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     SPOFA factors a real symmetric positive definite matrix.
%
%     SPOFA is usually called by SPOCO, but it can be called
%     directly with a saving in time if  RCOND  is not needed.
%     (Time for SPOCO) = (1 + 18/N)*(Time for SPOFA) .
%
%     On Entry
%
%        A       REAL(LDA, N)
%                the symmetric matrix to be factored.  Only the
%                diagonal and upper triangle are used.
%
%        LDA     INTEGER
%                the leading dimension of the array  A .
%
%        N       INTEGER
%                the order of the matrix  A .
%
%     On Return
%
%        A       an upper triangular matrix  R  so that  A = TRANS(R)*R
%                where  TRANS(R)  is the transpose.
%                The strict lower triangle is unaltered.
%                If  INFO .NE. 0 , the factorization is not complete.
%
%        INFO    INTEGER
%                = 0  for normal return.
%                = K  signals an error condition.  The leading minor
%                     of order  K  is not positive definite.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  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  SPOFA
persistent j jm1 k s t ; 

a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
%
if isempty(t), t=0; end;
if isempty(s), s=0; end;
if isempty(j), j=0; end;
if isempty(jm1), jm1=0; end;
if isempty(k), k=0; end;
%***FIRST EXECUTABLE STATEMENT  SPOFA
for j = 1 : n;
info = fix(j);
s = 0.0e0;
jm1 = fix(j - 1);
if( jm1>=1 )
for k = 1 : jm1;
t = a(k,j) - sdot(k-1,a(sub2ind(size(a),1,k):end),1,a(sub2ind(size(a),1,j):end),1);
t = t./a(k,k);
a(k,j) = t;
s = s + t.*t;
end; k = fix(jm1+1);
end;
s = a(j,j) - s;
if( s<=0.0e0 )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
return;
end;
a(j,j) = sqrt(s);
end; j = fix(n+1);
info = 0;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
end
%DECK SPOFS