Weekly post: Matlab ilutp behaviour disclosed and implemented in a m-file.

Here I am again.

This week has been more about researching than coding. I have finally been able to reproduce the output from the ilutp(ilu with threshold and pivoting) Matlab's algorithm with an m-script (named ILU_pc.m in my project's directory). The fact is that Matlab does not implement the algorithm as is described in  Yousef Saad's book in a few ways. Because of that I had to do  reverse engineering, testing many cases and matrices. That is the function, ugly as hell, but is just for testing purposes.



function [A, P] = ILU_pc(A, tau, thresh)

  B = A;
  n = length(A);
  P = speye(n);
  for i = 1:n
    for k = i:n
      A(k:n,k) *= thresh;
      A(k,k) /= thresh;
      [m,mi] = max(abs(A(k:n,k)))
      A(k,k) *= thresh;
      A(k:n,k) /= thresh;
      mi = mi + k -1;
      tmp = A(mi,:);
      A(mi,:) = A(k,:);
      A(k,:) = tmp;
      e = speye(n);
      e(mi,mi) = 0; e(k,mi) = 1;
      e(k,k) = 0; e(mi,k) = 1;
      P = e*P;
    endfor
    for k = 1:i-1
         if ( (A(i,k) == 0) || (abs(A(i,k)) < (tau*norm(B(:,k)))))
            A(i,k) = 0;
            continue
         endif
         A(i,k) = A(i,k) / A(k,k);
         A(i,k+1:n) = A(i,k+1:n) - A(i,k) * A(k,k+1:n);
    endfor
  endfor

  for i = 1:n
    for j = i+1:n
      if (abs(A(i,j)) < (tau*norm(B(:,j))))
        A(i,j) = 0;
      end
    end
  end
end

• The next goal to achieve is obviously to implement the function as .oct file translating this algorithm into a sparse one using Octave's internal data types. 

•  All the testing I did was at college using their Matlab license. That delayed me because I couldn't do almost nothing in the weekend. Now I have a function that reproduce the behavior of Matlab's version I can test against it my c++ code.

 See you next week!

The starting line.

As code period is starting today, I want to write a brief timeline for the first period of the GSOC here:

FIRST PERIOD

• 19 May-20 June: Implement ilu related functions (ilu0.cc, iluc.cc, ilutp.cc)  and merge them together with ilu.m script 

• 20-25 June: Automated test writing and documentation. Integration to mainstream octave code should be achieved here. 

• 27 June: (Millstone 1) ilu.m is fully functional and integrated with Octave core.

Taking the idea from Kai's last year blog, I will keep track of what is already done with the following figure.

Regarding repository setup, Kai helped me to configure a subrepository using bitbucket service. At present, it only contains an outdated Octave development version just to make sure things work.  For cloning:

hg clone https://edu159@bitbucket.org/edu159/octave-subrepo

However, I would not need to use this subrepo until the final integration of my code into Octave. For development purposes I have set another repository for daily work, as I am working with .oct files that compile standalone. Here is the repo you should check for my updated work.

hg clone https://edu159@bitbucket.org/edu159/gsoc2014-edu159


See you next week!