[ap,n,info]=dppfa(ap,n,info);

function [ap,n,info]=dppfa(ap,n,info);
%***BEGIN PROLOGUE  DPPFA
%***PURPOSE  Factor a real symmetric positive definite matrix stored in
%            packed form.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B1B
%***TYPE      doubleprecision (SPPFA-S, DPPFA-D, CPPFA-C)
%***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION, PACKED,
%             POSITIVE DEFINITE
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     DPPFA factors a doubleprecision symmetric positive definite
%     matrix stored in packed form.
%
%     DPPFA is usually called by DPPCO, but it can be called
%     directly with a saving in time if  RCOND  is not needed.
%     (time for DPPCO) = (1 + 18/N)*(time for DPPFA) .
%
%     On Entry
%
%        AP      doubleprecision (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 .
%
%     On Return
%
%        AP      an upper triangular matrix  R , stored in packed
%                form, so that  A = TRANS(R)*R .
%
%        INFO    INTEGER
%                = 0  for normal return.
%                = K  if the leading minor of order  K  is not
%                     positive definite.
%
%
%     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  DDOT
%***REVISION HISTORY  (YYMMDD)
%   780814  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)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DPPFA
persistent j jj jm1 k kj kk s t ; 

ap_shape=size(ap);ap=reshape(ap,1,[]);
%
if isempty(t), t=0; end;
if isempty(s), s=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(jm1), jm1=0; end;
if isempty(k), k=0; end;
if isempty(kj), kj=0; end;
if isempty(kk), kk=0; end;
%***FIRST EXECUTABLE STATEMENT  DPPFA
jj = 0;
for j = 1 : n;
info = fix(j);
s = 0.0d0;
jm1 = fix(j - 1);
kj = fix(jj);
kk = 0;
if( jm1>=1 )
for k = 1 : jm1;
kj = fix(kj + 1);
t = ap(kj) - ddot(k-1,ap(sub2ind(size(ap),max(kk+1,1)):end),1,ap(sub2ind(size(ap),max(jj+1,1)):end),1);
kk = fix(kk + k);
t = t./ap(kk);
ap(kj) = t;
s = s + t.*t;
end; k = fix(jm1+1);
end;
jj = fix(jj + j);
s = ap(jj) - s;
if( s<=0.0d0 )
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
return;
end;
ap(jj) = sqrt(s);
end; j = fix(n+1);
info = 0;
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
end
%DECK DPPGQ8