[a,lda,n,kpvt,info]=dsifa(a,lda,n,kpvt,info);

function [a,lda,n,kpvt,info]=dsifa(a,lda,n,kpvt,info);
%***BEGIN PROLOGUE  DSIFA
%***PURPOSE  Factor a real symmetric matrix by elimination with
%            symmetric pivoting.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B1A
%***TYPE      doubleprecision (SSIFA-S, DSIFA-D, CHIFA-C, CSIFA-C)
%***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION, SYMMETRIC
%***AUTHOR  Bunch, J., (UCSD)
%***DESCRIPTION
%
%     DSIFA factors a doubleprecision symmetric matrix by elimination
%     with symmetric pivoting.
%
%     To solve  A*X = B , follow DSIFA by DSISL.
%     To compute  INVERSE(A)*C , follow DSIFA by DSISL.
%     To compute  DETERMINANT(A) , follow DSIFA by DSIDI.
%     To compute  INERTIA(A) , follow DSIFA by DSIDI.
%     To compute  INVERSE(A) , follow DSIFA by DSIDI.
%
%     On Entry
%
%        A       doubleprecision(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       a block diagonal matrix and the multipliers which
%                were used to obtain it.
%                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 DSISL or DSIDI may
%                     divide by zero if called.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  DAXPY, DSWAP, IDAMAX
%***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  DSIFA
persistent absakk ak akm1 alpha bk bkm1 colmax denom imax imaxp1 j jj jmax k km1 km2 kstep mulk mulkm1 rowmax swap t ; 

kpvt_shape=size(kpvt);kpvt=reshape(kpvt,1,[]);
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
%
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(imax), imax=0; end;
if isempty(imaxp1), imaxp1=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(jmax), jmax=0; end;
if isempty(k), k=0; end;
if isempty(km1), km1=0; end;
if isempty(km2), km2=0; end;
if isempty(kstep), kstep=0; end;
if isempty(swap), swap=false; end;
%***FIRST EXECUTABLE STATEMENT  DSIFA
%
%     INITIALIZE
%
%     ALPHA IS USED IN CHOOSING PIVOT BLOCK SIZE.
%
alpha =(1.0d0+sqrt(17.0d0))./8.0d0;
%
info = 0;
%
%     MAIN LOOP ON K, WHICH GOES FROM N TO 1.
%
k = fix(n);
%
%        LEAVE THE LOOP IF K=0 OR K=1.
%
while( k~=0 );
if( k<=1 )
kpvt(1) = 1;
if( a(1,1)==0.0d0 )
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);
absakk = abs(a(k,k));
%
%        DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
%        COLUMN K.
%
[imax ,dumvar2,a(sub2ind(size(a),1,k):end)]=idamax(k-1,a(sub2ind(size(a),1,k):end),1);
colmax = abs(a(imax,k));
if( absakk<alpha.*colmax )
%
%           DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
%           ROW IMAX.
%
rowmax = 0.0d0;
imaxp1 = fix(imax + 1);
for j = imaxp1 : k;
rowmax = max(rowmax,abs(a(imax,j)));
end; j = fix(k+1);
if( imax~=1 )
[jmax ,dumvar2,a(sub2ind(size(a),1,imax):end)]=idamax(imax-1,a(sub2ind(size(a),1,imax):end),1);
rowmax = max(rowmax,abs(a(jmax,imax)));
end;
if( abs(a(imax,imax))>=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.0d0 )
%
%           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.
%
if( swap )
%
%              PERFORM AN INTERCHANGE.
%
[imax,dumvar2,dumvar3,dumvar4]=dswap(imax,a(sub2ind(size(a),1,imax):end),1,a(sub2ind(size(a),1,k-1):end),1);      a(sub2ind(size(a),1,imax):end)=dumvar2; a(sub2ind(size(a),1,k-1):end)=dumvar4; 
for jj = imax : km1;
j = fix(km1 + imax - jj);
t = a(j,k-1);
a(j,k-1) = a(imax,j);
a(imax,j) = t;
end; jj = fix(km1+1);
t = a(k-1,k);
a(k-1,k) = a(imax,k);
a(imax,k) = t;
end;
%
%           PERFORM THE ELIMINATION.
%
km2 = fix(k - 2);
if( km2~=0 )
ak = a(k,k)./a(k-1,k);
akm1 = a(k-1,k-1)./a(k-1,k);
denom = 1.0d0 - ak.*akm1;
for jj = 1 : km2;
j = fix(km1 - jj);
bk = a(j,k)./a(k-1,k);
bkm1 = a(j,k-1)./a(k-1,k);
mulk =(akm1.*bk-bkm1)./denom;
mulkm1 =(ak.*bkm1-bk)./denom;
t = mulk;
[j,t,dumvar3,dumvar4,dumvar5]=daxpy(j,t,a(sub2ind(size(a),1,k):end),1,a(sub2ind(size(a),1,j):end),1);      a(sub2ind(size(a),1,k):end)=dumvar3; a(sub2ind(size(a),1,j):end)=dumvar5; 
t = mulkm1;
[j,t,dumvar3,dumvar4,dumvar5]=daxpy(j,t,a(sub2ind(size(a),1,k-1):end),1,a(sub2ind(size(a),1,j):end),1);      a(sub2ind(size(a),1,k-1):end)=dumvar3; a(sub2ind(size(a),1,j):end)=dumvar5; 
a(j,k) = mulk;
a(j,k-1) = mulkm1;
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]=dswap(imax,a(sub2ind(size(a),1,imax):end),1,a(sub2ind(size(a),1,k):end),1);      a(sub2ind(size(a),1,imax):end)=dumvar2; a(sub2ind(size(a),1,k):end)=dumvar4; 
for jj = imax : k;
j = fix(k + imax - jj);
t = a(j,k);
a(j,k) = a(imax,j);
a(imax,j) = t;
end; jj = fix(k+1);
end;
%
%           PERFORM THE ELIMINATION.
%
for jj = 1 : km1;
j = fix(k - jj);
mulk = -a(j,k)./a(k,k);
t = mulk;
[j,t,dumvar3,dumvar4,dumvar5]=daxpy(j,t,a(sub2ind(size(a),1,k):end),1,a(sub2ind(size(a),1,j):end),1);      a(sub2ind(size(a),1,k):end)=dumvar3; a(sub2ind(size(a),1,j):end)=dumvar5; 
a(j,k) = mulk;
end; jj = fix(km1+1);
%
%           SET THE PIVOT ARRAY.
%
kpvt(k) = fix(k);
if( swap )
kpvt(k) = fix(imax);
end;
end;
k = fix(k - kstep);
end;
end;
kpvt_shape=zeros(kpvt_shape);kpvt_shape(:)=kpvt(1:numel(kpvt_shape));kpvt=kpvt_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
end %subroutine dsifa
%DECK DSILUR