function f=Determinant(C)
% DETERMINANT computes the determinant of a matrix
%   f=Determinant(C) computes the the determinant of C via
%   LU-decomposition with partial pivoting

n=length(C);
f=1; 
for i=1:n
  [cmax,kmax]=max(abs(C(i:n,i)));
   if cmax==0                                        % Matrix singular
     error('Matrix is singular')
   end
   kmax=kmax+i-1;
   if kmax~=i
     h=C(i,:); C(i,:)=C(kmax,:); C(kmax,:)=h;
     f=-f; 
   end
   f=f*C(i,i);                                       % elimination step
   C(i+1:n,i)=C(i+1:n,i)/C(i,i);     
   C(i+1:n,i+1:n)=C(i+1:n,i+1:n)-C(i+1:n,i)*C(i,i+1:n);
end