Hi Eduardo,
Today I had a closer look to your past weeks work. Your code looks good, but it definitely needs at least a little documentation. Sure, this is very annoying, especially if you are still in an very early stage, but I wasn't able to verify your ILUTP function, I'm also very interested in. I tried to use your ILUTP function this way directly from MATLAB:
options.type = 'ilutp';
options.thresh = 0.1;
options.droptol = 1e-10;
% Generate some diagonal dominant matrix A
A = sprand (10,10,0.1);
A = A + diag (rand (10,1) + 1);
A = sparse (A)
% Call MATLAB ILUTP
[L,U,P] = ilu (A,options)
% Call new ILUTP
[A2,P] = ILU_pc (A,1e-10,0.1)
Did I use your function correctly? I think you are trying to optimize your function doing some in-place calculation in A. This is not a very common programming pattern, as the user might want to keep the original A even after doing the factorization step.
Maybe you can post a tiny usage example with small comments in your blog?
Best,
Kai