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: [EMAIL PROTECTED] ** _______________________________________________ Numpy-discussion mailing list [email protected] http://projects.scipy.org/mailman/listinfo/numpy-discussion
