[ap,n,det,job]=dppdi(ap,n,det,job);

function [ap,n,det,job]=dppdi(ap,n,det,job);
%***BEGIN PROLOGUE  DPPDI
%***PURPOSE  Compute the determinant and inverse of a real symmetric
%            positive definite matrix using factors from DPPCO or DPPFA.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B1B, D3B1B
%***TYPE      doubleprecision (SPPDI-S, DPPDI-D, CPPDI-C)
%***KEYWORDS  DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK, MATRIX,
%             PACKED, POSITIVE DEFINITE
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     DPPDI computes the determinant and inverse
%     of a doubleprecision symmetric positive definite matrix
%     using the factors computed by DPPCO or DPPFA .
%
%     On Entry
%
%        AP      doubleprecision (N*(N+1)/2)
%                the output from DPPCO or DPPFA.
%
%        N       INTEGER
%                the order of the matrix  A .
%
%        JOB     INTEGER
%                = 11   both determinant and inverse.
%                = 01   inverse only.
%                = 10   determinant only.
%
%     On Return
%
%        AP      the upper triangular half of the inverse .
%                The strict lower triangle is unaltered.
%
%        DET     doubleprecision(2)
%                determinant of original matrix if requested.
%                Otherwise not referenced.
%                DETERMINANT = DET(1) * 10.0**DET(2)
%                with  1.0 .LE. DET(1) .LT. 10.0
%                or  DET(1) .EQ. 0.0 .
%
%     Error Condition
%
%        A division by zero will occur if the input factor contains
%        a zero on the diagonal and the inverse is requested.
%        It will not occur if the subroutines are called correctly
%        and if DPOCO or DPOFA has set INFO .EQ. 0 .
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  DAXPY, DSCAL
%***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  DPPDI
persistent i ii j j1 jj jm1 k k1 kj kk kp1 s t ; 

ap_shape=size(ap);ap=reshape(ap,1,[]);
%
if isempty(t), t=0; end;
if isempty(s), s=0; end;
if isempty(i), i=0; end;
if isempty(ii), ii=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(jm1), jm1=0; end;
if isempty(j1), j1=0; end;
if isempty(k), k=0; end;
if isempty(kj), kj=0; end;
if isempty(kk), kk=0; end;
if isempty(kp1), kp1=0; end;
if isempty(k1), k1=0; end;
%***FIRST EXECUTABLE STATEMENT  DPPDI
%
%     COMPUTE DETERMINANT
%
if( fix(job./10)~=0 )
det(1) = 1.0d0;
det(2) = 0.0d0;
s = 10.0d0;
ii = 0;
for i = 1 : n;
ii = fix(ii + i);
det(1) = ap(ii).^2.*det(1);
if( det(1)==0.0d0 )
break;
end;
while( det(1)<1.0d0 );
det(1) = s.*det(1);
det(2) = det(2) - 1.0d0;
end;
while( det(1)>=s );
det(1) = det(1)./s;
det(2) = det(2) + 1.0d0;
end;
end;
end;
%
%     COMPUTE INVERSE(R)
%
if( rem(job,10)~=0 )
kk = 0;
for k = 1 : n;
k1 = fix(kk + 1);
kk = fix(kk + k);
ap(kk) = 1.0d0./ap(kk);
t = -ap(kk);
[dumvar1,t,ap(sub2ind(size(ap),max(k1,1)):end)]=dscal(k-1,t,ap(sub2ind(size(ap),max(k1,1)):end),1);
kp1 = fix(k + 1);
j1 = fix(kk + 1);
kj = fix(kk + k);
if( n>=kp1 )
for j = kp1 : n;
t = ap(kj);
ap(kj) = 0.0d0;
[k,t,dumvar3,dumvar4,dumvar5]=daxpy(k,t,ap(sub2ind(size(ap),max(k1,1)):end),1,ap(sub2ind(size(ap),max(j1,1)):end),1);   dumvar3i=find((ap(sub2ind(size(ap),max(k1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(j1,1)):end))~=(dumvar5));   ap(k1-1+dumvar3i)=dumvar3(dumvar3i); ap(j1-1+dumvar5i)=dumvar5(dumvar5i); 
j1 = fix(j1 + j);
kj = fix(kj + j);
end; j = fix(n+1);
end;
end; k = fix(n+1);
%
%        FORM  INVERSE(R) * TRANS(INVERSE(R))
%
jj = 0;
for j = 1 : n;
j1 = fix(jj + 1);
jj = fix(jj + j);
jm1 = fix(j - 1);
k1 = 1;
kj = fix(j1);
if( jm1>=1 )
for k = 1 : jm1;
t = ap(kj);
[k,t,dumvar3,dumvar4,dumvar5]=daxpy(k,t,ap(sub2ind(size(ap),max(j1,1)):end),1,ap(sub2ind(size(ap),max(k1,1)):end),1);   dumvar3i=find((ap(sub2ind(size(ap),max(j1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(k1,1)):end))~=(dumvar5));   ap(j1-1+dumvar3i)=dumvar3(dumvar3i); ap(k1-1+dumvar5i)=dumvar5(dumvar5i); 
k1 = fix(k1 + k);
kj = fix(kj + 1);
end; k = fix(jm1+1);
end;
t = ap(jj);
[j,t,ap(sub2ind(size(ap),max(j1,1)):end)]=dscal(j,t,ap(sub2ind(size(ap),max(j1,1)):end),1);
end; j = fix(n+1);
end;
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
end %subroutine dppdi
%DECK DPPERM