[ap,n,kpvt,det,inert,work,job]=dspdi(ap,n,kpvt,det,inert,work,job);

function [ap,n,kpvt,det,inert,work,job]=dspdi(ap,n,kpvt,det,inert,work,job);
%***BEGIN PROLOGUE  DSPDI
%***PURPOSE  Compute the determinant, inertia, inverse of a real
%            symmetric matrix stored in packed form using the factors
%            from DSPFA.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B1A, D3B1A
%***TYPE      doubleprecision (SSPDI-S, DSPDI-D, CHPDI-C, CSPDI-C)
%***KEYWORDS  DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK, MATRIX,
%             PACKED, SYMMETRIC
%***AUTHOR  Bunch, J., (UCSD)
%***DESCRIPTION
%
%     DSPDI computes the determinant, inertia and inverse
%     of a doubleprecision symmetric matrix using the factors from
%     DSPFA, where the matrix is stored in packed form.
%
%     On Entry
%
%        AP      doubleprecision (N*(N+1)/2)
%                the output from DSPFA.
%
%        N       INTEGER
%                the order of the matrix A.
%
%        KPVT    INTEGER(N)
%                the pivot vector from DSPFA.
%
%        WORK    doubleprecision(N)
%                work vector.  Contents ignored.
%
%        JOB     INTEGER
%                JOB has the decimal expansion  ABC  where
%                   if  C .NE. 0, the inverse is computed,
%                   if  B .NE. 0, the determinant is computed,
%                   if  A .NE. 0, the inertia is computed.
%
%                For example, JOB = 111  gives all three.
%
%     On Return
%
%        Variables not requested by JOB are not used.
%
%        AP     contains the upper triangle of the inverse of
%               the original matrix, stored in packed form.
%               The columns of the upper triangle are stored
%               sequentially in a one-dimensional array.
%
%        DET    doubleprecision(2)
%               determinant of original matrix.
%               DETERMINANT = DET(1) * 10.0**DET(2)
%               with 1.0 .LE. ABS(DET(1)) .LT. 10.0
%               or DET(1) = 0.0.
%
%        INERT  INTEGER(3)
%               the inertia of the original matrix.
%               INERT(1)  =  number of positive eigenvalues.
%               INERT(2)  =  number of negative eigenvalues.
%               INERT(3)  =  number of zero eigenvalues.
%
%     Error Condition
%
%        A division by zero will occur if the inverse is requested
%        and  DSPCO  has set RCOND .EQ. 0.0
%        or  DSPFA  has set  INFO .NE. 0 .
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  DAXPY, DCOPY, DDOT, DSWAP
%***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  DSPDI
persistent ak akkp1 akp1 d ij ik ikp1 iks j jb jk jkp1 k kk kkp1 km1 ks ksj kskp1 kstep nodet noert noinv t temp ten ; 

ap_shape=size(ap);ap=reshape(ap,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
kpvt_shape=size(kpvt);kpvt=reshape(kpvt,1,[]);
%
if isempty(akkp1), akkp1=0; end;
if isempty(temp), temp=0; end;
if isempty(ten), ten=0; end;
if isempty(d), d=0; end;
if isempty(t), t=0; end;
if isempty(ak), ak=0; end;
if isempty(akp1), akp1=0; end;
if isempty(ij), ij=0; end;
if isempty(ik), ik=0; end;
if isempty(ikp1), ikp1=0; end;
if isempty(iks), iks=0; end;
if isempty(j), j=0; end;
if isempty(jb), jb=0; end;
if isempty(jk), jk=0; end;
if isempty(jkp1), jkp1=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(kkp1), kkp1=0; end;
if isempty(km1), km1=0; end;
if isempty(ks), ks=0; end;
if isempty(ksj), ksj=0; end;
if isempty(kskp1), kskp1=0; end;
if isempty(kstep), kstep=0; end;
if isempty(noinv), noinv=false; end;
if isempty(nodet), nodet=false; end;
if isempty(noert), noert=false; end;
%***FIRST EXECUTABLE STATEMENT  DSPDI
noinv = rem(job,10)==0;
nodet = rem(job,100)./10==0;
noert = rem(job,1000)./100==0;
%
if( ~(nodet && noert) )
if( ~(noert) )
inert(1) = 0;
inert(2) = 0;
inert(3) = 0;
end;
if( ~(nodet) )
det(1) = 1.0d0;
det(2) = 0.0d0;
ten = 10.0d0;
end;
t = 0.0d0;
ik = 0;
for k = 1 : n;
kk = fix(ik + k);
d = ap(kk);
%
%           CHECK IF 1 BY 1
%
if( kpvt(k)<=0 )
%
%              2 BY 2 BLOCK
%              use DET (D  S)  =  (D/T * C - T) * T  ,  T = ABS(S)
%                      (S  C)
%              TO AVOID UNDERFLOW/OVERFLOW TROUBLES.
%              TAKE TWO PASSES THROUGH SCALING.  USE  T  FOR FLAG.
%
if( t==0.0d0 )
ikp1 = fix(ik + k);
kkp1 = fix(ikp1 + k);
t = abs(ap(kkp1));
d =(d./t).*ap(kkp1+1) - t;
else;
d = t;
t = 0.0d0;
end;
end;
%
if( ~(noert) )
if( d>0.0d0 )
inert(1) = fix(inert(1) + 1);
end;
if( d<0.0d0 )
inert(2) = fix(inert(2) + 1);
end;
if( d==0.0d0 )
inert(3) = fix(inert(3) + 1);
end;
end;
%
if( ~(nodet) )
det(1) = d.*det(1);
if( det(1)~=0.0d0 )
while( abs(det(1))<1.0d0 );
det(1) = ten.*det(1);
det(2) = det(2) - 1.0d0;
end;
while( abs(det(1))>=ten );
det(1) = det(1)./ten;
det(2) = det(2) + 1.0d0;
end;
end;
end;
ik = fix(ik + k);
end; k = fix(n+1);
end;
%
%     COMPUTE INVERSE(A)
%
if( ~(noinv) )
k = 1;
ik = 0;
while( k<=n );
km1 = fix(k - 1);
kk = fix(ik + k);
ikp1 = fix(ik + k);
kkp1 = fix(ikp1 + k);
if( kpvt(k)<0 )
%
%              2 BY 2
%
t = abs(ap(kkp1));
ak = ap(kk)./t;
akp1 = ap(kkp1+1)./t;
akkp1 = ap(kkp1)./t;
d = t.*(ak.*akp1-1.0d0);
ap(kk) = akp1./d;
ap(kkp1+1) = ak./d;
ap(kkp1) = -akkp1./d;
if( km1>=1 )
[km1,ap(sub2ind(size(ap),max(ikp1+1,1)):end),dumvar3,work]=dcopy(km1,ap(sub2ind(size(ap),max(ikp1+1,1)):end),1,work,1);
ij = 0;
for j = 1 : km1;
jkp1 = fix(ikp1 + j);
[ap(jkp1) ,j,ap(sub2ind(size(ap),max(ij+1,1)):end),dumvar4,work]=ddot(j,ap(sub2ind(size(ap),max(ij+1,1)):end),1,work,1);
[dumvar1,work(j),dumvar3,dumvar4,dumvar5]=daxpy(j-1,work(j),ap(sub2ind(size(ap),max(ij+1,1)):end),1,ap(sub2ind(size(ap),max(ikp1+1,1)):end),1);   dumvar3i=find((ap(sub2ind(size(ap),max(ij+1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(ikp1+1,1)):end))~=(dumvar5));   ap(ij+1-1+dumvar3i)=dumvar3(dumvar3i); ap(ikp1+1-1+dumvar5i)=dumvar5(dumvar5i); 
ij = fix(ij + j);
end; j = fix(km1+1);
ap(kkp1+1) = ap(kkp1+1) + ddot(km1,work,1,ap(sub2ind(size(ap),max(ikp1+1,1)):end),1);
ap(kkp1) = ap(kkp1) + ddot(km1,ap(sub2ind(size(ap),max(ik+1,1)):end),1,ap(sub2ind(size(ap),max(ikp1+1,1)):end),1);
[km1,ap(sub2ind(size(ap),max(ik+1,1)):end),dumvar3,work]=dcopy(km1,ap(sub2ind(size(ap),max(ik+1,1)):end),1,work,1);
ij = 0;
for j = 1 : km1;
jk = fix(ik + j);
[ap(jk) ,j,ap(sub2ind(size(ap),max(ij+1,1)):end),dumvar4,work]=ddot(j,ap(sub2ind(size(ap),max(ij+1,1)):end),1,work,1);
[dumvar1,work(j),dumvar3,dumvar4,dumvar5]=daxpy(j-1,work(j),ap(sub2ind(size(ap),max(ij+1,1)):end),1,ap(sub2ind(size(ap),max(ik+1,1)):end),1);   dumvar3i=find((ap(sub2ind(size(ap),max(ij+1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(ik+1,1)):end))~=(dumvar5));   ap(ij+1-1+dumvar3i)=dumvar3(dumvar3i); ap(ik+1-1+dumvar5i)=dumvar5(dumvar5i); 
ij = fix(ij + j);
end; j = fix(km1+1);
ap(kk) = ap(kk) + ddot(km1,work,1,ap(sub2ind(size(ap),max(ik+1,1)):end),1);
end;
kstep = 2;
else;
%
%              1 BY 1
%
ap(kk) = 1.0d0./ap(kk);
if( km1>=1 )
[km1,ap(sub2ind(size(ap),max(ik+1,1)):end),dumvar3,work]=dcopy(km1,ap(sub2ind(size(ap),max(ik+1,1)):end),1,work,1);
ij = 0;
for j = 1 : km1;
jk = fix(ik + j);
[ap(jk) ,j,ap(sub2ind(size(ap),max(ij+1,1)):end),dumvar4,work]=ddot(j,ap(sub2ind(size(ap),max(ij+1,1)):end),1,work,1);
[dumvar1,work(j),dumvar3,dumvar4,dumvar5]=daxpy(j-1,work(j),ap(sub2ind(size(ap),max(ij+1,1)):end),1,ap(sub2ind(size(ap),max(ik+1,1)):end),1);   dumvar3i=find((ap(sub2ind(size(ap),max(ij+1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(ik+1,1)):end))~=(dumvar5));   ap(ij+1-1+dumvar3i)=dumvar3(dumvar3i); ap(ik+1-1+dumvar5i)=dumvar5(dumvar5i); 
ij = fix(ij + j);
end; j = fix(km1+1);
ap(kk) = ap(kk) + ddot(km1,work,1,ap(sub2ind(size(ap),max(ik+1,1)):end),1);
end;
kstep = 1;
end;
%
%           SWAP
%
ks = fix(abs(kpvt(k)));
if( ks~=k )
iks =fix(fix((ks.*(ks-1))./2));
[ks,dumvar2,dumvar3,dumvar4]=dswap(ks,ap(sub2ind(size(ap),max(iks+1,1)):end),1,ap(sub2ind(size(ap),max(ik+1,1)):end),1);   dumvar2i=find((ap(sub2ind(size(ap),max(iks+1,1)):end))~=(dumvar2));dumvar4i=find((ap(sub2ind(size(ap),max(ik+1,1)):end))~=(dumvar4));   ap(iks+1-1+dumvar2i)=dumvar2(dumvar2i); ap(ik+1-1+dumvar4i)=dumvar4(dumvar4i); 
ksj = fix(ik + ks);
for jb = ks : k;
j = fix(k + ks - jb);
jk = fix(ik + j);
temp = ap(jk);
ap(jk) = ap(ksj);
ap(ksj) = temp;
ksj = fix(ksj -(j-1));
end; jb = fix(k+1);
if( kstep~=1 )
kskp1 = fix(ikp1 + ks);
temp = ap(kskp1);
ap(kskp1) = ap(kkp1);
ap(kkp1) = temp;
end;
end;
ik = fix(ik + k);
if( kstep==2 )
ik = fix(ik + k + 1);
end;
k = fix(k + kstep);
end;
end;
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
kpvt_shape=zeros(kpvt_shape);kpvt_shape(:)=kpvt(1:numel(kpvt_shape));kpvt=kpvt_shape;
end %subroutine dspdi
%DECK DSPENC