[SciPy-dev] the current state and future directions of the sparse solvers

[Nathan Bell]

On Tue, Apr 8, 2008 at 3:46 PM, Ondrej Certik <ondrej@certik.cz> wrote:
>  >  Supporting a common interface to the various solvers would be quite
>  >  useful.  However, like (2) we should not support interfaces to
>  >  optional libraries/methods.
>
>  You mean support in scipy by default, right? I would like to have the
>  unified interface to all solvers (including petsc and others).

I think that a common interface to UMFPACK/PETSc/etc. should live a scikit.

Ideally, the scikit would include those libraries itself.

>  OK, that sounds very good. I am going to study how to start such a
>  scikit and we are going to move our additional solvers we have in
>  sfepy to it (pysparse, soon petsc4py) -- are you ok with that Robert?
>  :) Plus umfpack.  Plus the eigensolvers, like primme, blzpack and
>  others.

I just recently saw a talk by some people that work on providing
abstract interfaces to the solvers in Trilinos:
http://trilinos.sandia.gov/packages/docs/r8.0/packages/stratimikos/doc/html/index.html
http://trilinos.sandia.gov/packages/belos/

I have no experience with these, but they may provide some guidance.

>  And then we'll see how it goes and if some part is useful (like
>  pysparse or blzpack) and with BSD license, it can later be moved to
>  scipy directly.

Absolutely.

>  Do you have some idea about the interface? I noticed the
>  linalg/interface.py, that's the way to go. This is only for the
>  iterative solvers though.

LinearOperator was designed to formalize the interface expected by
isolve/iterative.py

I haven't thought about it much since then, so it's possible that
additions/modifications will be necessary.


>  We need to devise some unified interface for the direct solvers and
>  also for eigensolvers. Ok, I'll try to follow the arpack example and
>  then we'll see if it's ok, or it needs some changes for example in the
>  blzpack case.
>
>  For example the lobpcg now in scipy has a different interface, doesn't it?

Only slightly. Since lobpcg operates on block vectors I added a
.matmat() to LinearOperator.  If a operator provides a matvec()
method, but no matmat(), then the matvec() is applied to each column.

In the future some sparse matrices may have optimized matmat() for blockvectors.


-- 
Nathan Bell wnbell@gmail.com
http://graphics.cs.uiuc.edu/~wnbell/