[Numpy-discussion] Matrix Class [was numpy release]

[Bill Spotz]

On Apr 29, 2008, at 2:50 AM, Robert Cimrman wrote:

> FYI:
> In scipy, the iterative solvers (cg and like) should work with
> LinearOperator (a class with matvec, rmatvec, and eventually matmat
> methods) passed instead of the A matrix argument. There is a function
> aslinearoperator() is scipy.sparse.linalg, that constructs a
> LinearOperator instance from a numpy array, matrix and scipy sparse
> matrix. It could be generalized to accept a function for a matrix-free
> computation, too. Also LinearOperator.__mul__, __rmul__ could be  
> easily
> defined via matvec, rmatvec.
>
> IMHO the all iterative solvers could use such a concept, to shield  
> them
> from what the system matrix A or the preconditioner actually are and  
> to
> be able writing them in a visually appealing way (i.e. close to what  
> one
> writes on a paper).


I agree.  This is very close to what we use in Sandia's C++ Trilinos  
project: http://trilinos.sandia.gov.  The linear algebra services  
package (Epetra: http://trilinos.sandia.gov/packages/epetra) supports  
the following class hierarchies:

Epetra_Operator <- Epetra_RowMatrix <- Epetra_CrsMatrix
Epetra_MultiVector <- Epetra_Vector

(Where "Crs" stands for "compressed row storage").  A typical solver  
package (e.g. AztecOO: http://trilinos.sandia.gov/packages/aztecoo)  
defines solver routines that take an Epetra_Operator for A (and  
preconditioner P) and an Epetra_MultiVector for b and x.  The most  
common way to call them is with an Epetra_CrsMatrix and an  
Epetra_Vector, but the flexibility is there.

** Bill Spotz                                              **
** Sandia National Laboratories  Voice: (505)845-0170      **
** P.O. Box 5800                 Fax:   (505)284-0154      **
** Albuquerque, NM 87185-0370    Email: wfspotz at sandia.gov **