Firm pencils down.

During this week I have been reorganizing all the code, docs and tests in a better way for integrating into Octave. As Rik kindly suggested, I decided to organize things this way:

• Inside libinterp/dldfcn directory I have created two files, __ichol__.cc and __ilu__.cc

•  Within those files there are the dld functions that implements the each of the algorithms. They are ment to be built-in functions and follows the __foo__.cc naming convention.

              * __ilu__.cc: contains __ilu0__() , __iluc__() and __ilutp__()

              * __ichol__.cc: contains __ichol0__() and __icholt__().

• I have moved all the tests from .cc files to .m scripts so no tests are performed for built-in functions.

 The code is ready to be pulled from my repo to be reviewed :

 https://edu159@bitbucket.org/edu159/octave-edu159

Practical examples: preconditioned conjugate gradient

It is interesting to show how preconditioning techniques can improve the convergency of some iterative solvers. In this case I am running a Matlab example using the Poisson matrix (that is positive definite) obtained with gallery() function. The scritp:

1. In this first case the convergency of pcg using ICHOL(0) algorithm, modified ICHOL(0) algorithm and without preconditioning are compared.

N = 100;

A = gallery ('Poisson', N);

b = ones (size (A, 1), 1);

tol = 1e-6; maxit = 100;

[x0, fl0, rr0, it0, rv0] = pcg (A, b, tol, maxit);

L1 = ichol (A);

[x1, fl1, rr1, it1, rv1] = pcg(A, b, tol, maxit, L1, L1');

opts.type = 'nofill'; opts.michol = 'on';

L2 = ichol (A, opts);

e = ones (size (A, 2), 1);

norm(A * e - L2 * (L2' * e))

[x2, fl2, rr2, it2, rv2] = pcg (A, b, tol, maxit, L2, L2');

semilogy (0:maxit, rv0 ./ norm (b), 'b.');

hold on;

semilogy (0:it1, rv1 ./ norm(b), 'r.');

semilogy (0:it2, rv2 ./ norm(b), 'k.');

xlabel ('iterations');

ylabel ('error');

legend ('No Preconditioner', 'IC(0)', 'MIC(0)'); 

Octave

Matlab

2. In this second part of the script what is compared is the      convergency of ICHOLT algorithm with different values of droptol.

 

L3 = ichol(A, struct('type', 'ict', 'droptol', 1e-1));

[x3, fl3, rr3, it3, rv3] = pcg (A, b, tol, maxit, L3, L3');

L4 = ichol (A, struct ('type', 'ict', 'droptol', 1e-2));

[x4, fl4, rr4, it4, rv4] = pcg (A, b, tol, maxit, L4, L4');

L5 = ichol (A, struct ('type', 'ict', 'droptol', 1e-3));

[x5, fl5, rr5, it5, rv5] = pcg (A, b, tol, maxit, L5, L5');

figure; semilogy (0:maxit, rv0 ./ norm (b), 'b-', 'linewidth', 2);

hold on;

semilogy (0:it3, rv3 ./ norm(b), 'b-.', 'linewidth', 2);

semilogy (0:it4, rv4 ./ norm(b), 'b--', 'linewidth', 2);

semilogy (0:it5, rv5 ./ norm(b), 'b:', 'linewidth', 2);

ylabel ('error');

xlabel ('iterations');

legend ('No Preconditioner', 'ICT(1e-1)', 'ICT(1e-2)', ...

   'ICT(1e-3)', 'Location', 'SouthEast');

Octave

Matlab

As it can be seen Octave plots are the same as Matlab's ones.  Both lead to a decrease in the number of steps upt to convergence of the pcg method. ILU algorithms could also have been used here, but due to the simetry of the problem matrix ICHOL is faster.

Regards,

Eduardo