Re: [eigen] initial sparse PCG solver

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi,

excellent!

I don't see any obvious way to optimize the implementation w/r to the
performance.

So if we put aside pure API discussion, I only have remarks you
probably already know:
- It should accept for the rhs and the result any objects
(sparse/dense, matrix/vector).  That means replacing the dots by
matrix products.
- The EPSILON must be a user parameter with a default value depending
on the scalar type.
- Check it works with complexes.
- The two divisions should be checked for 0.
- For preconditioners you can also try our SparseLLT with a relatively
large epsilon such that it is really fast and stay sparse.
- The preconditionner should be an arbitrary type such that we could
use "virtual" matrix, i.e., we would wrap SparseLLT into a tiny struct
having operator* overloaded doing:

XXXX operator* (const T& x) const {
  return m_llt.solve(x);
}

but this first requires we port the solve API we now have in Dense to
Sparse, i.e., returning proxy objects.... Not difficult, mostly
copy/pasta from the dense LLT dec for instance.

cheers,
Gael.

On Fri, Jul 23, 2010 at 11:40 AM, Daniel Lowengrub <lowdanie@xxxxxxxxx> wrote:
> Hi,
> I started working on a preconditioned conjugate gradient solver for
> the Sparse module.  Before working on the interface and eigen
> integration I decided to try and get the basic functionality setup.  I
> attached a file with the implementation of the algorithm and a small
> test program.  The test program uses the simple jacobi preconditioner.
>  Obviously better preconditioners will have to be implemented in order
> for the PCG algorithm to be feasible, but for now I'd like to focus on
> the algorithm itself.
> Any feedback on the implementation would be appreciated.
> The implementation used here is based on this paper:
> www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
>
> Daniel
>