[a,lda,n,ipvt,info]=dgefa(a,lda,n,ipvt,info);

function [a,lda,n,ipvt,info]=dgefa(a,lda,n,ipvt,info);
%***BEGIN PROLOGUE  DGEFA
%***PURPOSE  Factor a matrix using Gaussian elimination.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2A1
%***TYPE      doubleprecision (SGEFA-S, DGEFA-D, CGEFA-C)
%***KEYWORDS  GENERAL MATRIX, LINEAR ALGEBRA, LINPACK,
%             MATRIX FACTORIZATION
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     DGEFA factors a doubleprecision matrix by Gaussian elimination.
%
%     DGEFA is usually called by DGECO, but it can be called
%     directly with a saving in time if  RCOND  is not needed.
%     (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) .
%
%     On Entry
%
%        A       doubleprecision(LDA, N)
%                the matrix to be factored.
%
%        LDA     INTEGER
%                the leading dimension of the array  A .
%
%        N       INTEGER
%                the order of the matrix  A .
%
%     On Return
%
%        A       an upper triangular matrix and the multipliers
%                which were used to obtain it.
%                The factorization can be written  A = L*U  where
%                L  is a product of permutation and unit lower
%                triangular matrices and  U  is upper triangular.
%
%        IPVT    INTEGER(N)
%                an integer vector of pivot indices.
%
%        INFO    INTEGER
%                = 0  normal value.
%                = K  if  U(K,K) .EQ. 0.0 .  This is not an error
%                     condition for this subroutine, but it does
%                     indicate that DGESL or DGEDI will divide by zero
%                     if called.  Use  RCOND  in DGECO for a reliable
%                     indication of singularity.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  DAXPY, DSCAL, IDAMAX
%***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  DGEFA
persistent j k kp1 l nm1 t ; 

ipvt_shape=size(ipvt);ipvt=reshape(ipvt,1,[]);
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(j), j=0; end;
if isempty(k), k=0; end;
if isempty(kp1), kp1=0; end;
if isempty(l), l=0; end;
if isempty(nm1), nm1=0; end;
%
%     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
%
%***FIRST EXECUTABLE STATEMENT  DGEFA
info = 0;
nm1 = fix(n - 1);
if( nm1>=1 )
for k = 1 : nm1;
kp1 = fix(k + 1);
%
%        FIND L = PIVOT INDEX
%
l = fix(idamax(n-k+1,a(sub2ind(size(a),k,k):end),1) + k - 1);
ipvt(k) = fix(l);
%
%        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
%
if( a(l,k)==0.0d0 )
info = fix(k);
else;
%
%           INTERCHANGE IF NECESSARY
%
if( l~=k )
t = a(l,k);
a(l,k) = a(k,k);
a(k,k) = t;
end;
%
%           COMPUTE MULTIPLIERS
%
t = -1.0d0./a(k,k);
[dumvar1,t,a(sub2ind(size(a),k+1,k):end)]=dscal(n-k,t,a(sub2ind(size(a),k+1,k):end),1);
%
%           ROW ELIMINATION WITH COLUMN INDEXING
%
for j = kp1 : n;
t = a(l,j);
if( l~=k )
a(l,j) = a(k,j);
a(k,j) = t;
end;
[dumvar1,t,dumvar3,dumvar4,dumvar5]=daxpy(n-k,t,a(sub2ind(size(a),k+1,k):end),1,a(sub2ind(size(a),k+1,j):end),1);      a(sub2ind(size(a),k+1,k):end)=dumvar3; a(sub2ind(size(a),k+1,j):end)=dumvar5; 
end; j = fix(n+1);
end;
end; k = fix(nm1+1);
end;
ipvt(n) = fix(n);
if( a(n,n)==0.0d0 )
info = fix(n);
end;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
end
%DECK DGEFS