function [L,U]=mluinc_matlab(A,degree);
%General modified incomplete factorization
%degree of accuracy of the factorization specified
%be aware of modification parameter for 2D and 3D PDE problems
n=length(A);
si=10/n; %2D Poisson problem
%si=5/n^(2/3); %3D Poisson problem
R=sparse(zeros(n,n));
for i=1:degree
  B=A;
  if i==degree B=A+si*sparse(diag(diag(A))); end
  for k=1:n-1,
      B(k+1:n,k)=B(k+1:n,k)/B(k,k);
      C=B(k+1:n,k+1:n);C=(C~=0|R(k+1:n,k+1:n)~=0);D=B(k+1:n,k)*B(k,k+1:n);
      B(k+1:n,k+1:n)=B(k+1:n,k+1:n)-C.*D;
      if i==degree 
        E=(~C)&(D~=0);
        modify=(E.*D)*ones(n-k,1);
        B(k+1:n,k+1:n)=B(k+1:n,k+1:n)-diag(modify);
      end
  end
  U=triu(B);
  L=tril(B,-1)+speye(n);
  R=R+L*U-A;
end
