On Mon 2008-10-20 11:23, Roy Stogner wrote: > > > On Mon, 20 Oct 2008, John Peterson wrote: > > > On Mon, Oct 20, 2008 at 11:14 AM, Roy Stogner <[EMAIL PROTECTED]> wrote: > >>>> Also, I would suggest to derive SparseMatrix from ShellMatrix rather than > >>>> vice versa. This could actually make SparseShellMatrix unnecessary. > >>>> This > >>>> is of course an advantage of the unification. > >>> > >>> This one I disagree with. Inheritance should always follow the "is a" > >>> organizational semantic, and a SparseMatrix (as we now define it) is > >>> most definitely *not* a ShellMatrix. > >> > >> Yes, but both of them are a MatrixBase. > >> > >> Yes, it is easier to add new classes when someone else has already > >> volunteered to do the work of writing them. ;-) > > > > In all seriousness, I still don't see the compelling argument for this > > extra base class. > > Because SparseMatrix is-not-a ShellMatrix (doesn't need to attach or > even acknowledge the existance of a user's matvec shell function), and > ShellMatrix is-not-a SparseMatrix (doesn't need sparsity pattern > details in the constructor, can't return or set individual matrix > entries), but both of them are-a MatrixBase (can perform matrix-vector > products on NumericVectors, at least, and with some tricks we can pass > them both to what are conceptually the same solvers, scale them, add > them together, etc).
Why can't a shell matrix implement these operations? Certainly a PETSc MATSHELL is a first-class matrix. All these operations are implementable by a shell matrix. http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/include/petscmat.h.html#line1309 In fact, it would be elegant to represent Tim's rank-1 perturbation as a matrix type that requires preallocation and allows the user to set values for the sparse stored part (i.e. it's a first-class matrix by any definition) but implements one extra function to set the update vectors. How about just adding SparseMatrix::petsc_mat() which returns a PETSc Mat for any matrix type. For a Shell matrix, this would return the internal object set up with MatCreateShell. For a PetscMatrix, this would be equivalent to PetscMatrix::mat() (which could be removed). This way, EpetraMatrix could implement petsc_mat() using MatShell which would enable one to use PETSc solvers with Epetra matrices. The converse would also be possible by adding SparseMatrix::epetra_mat() which could be implemented for PetscMatrix using Epetra's matrix-free abstraction (called RowMatrix IIRC). With this rearrangement, it might be worth dropping `Sparse' from SparseMatrix. If your shell matrix doesn't use preallocation, the implementation can just ignore it, I don't see it as justification for an awkward class hierarchy. To be clear, I see ShellMatrix derived from Matrix and implementing petsc_mat() (and similarly epetra_mat()) by calling MatShellSetOperation(...) for all implemented operations. The user would implement a ShellMatrix by deriving from it and implementing whatever functions are appropriate. Unfortunately with C++, we can't elegantly determine if a function is implemented (with C we check if a function pointer is non-NULL) so the derived class needs to keep a `table' consistent with which methods are implemented. (The type system is in the way here. We can't make the methods pure virtual because the compiler will complain that we haven't implemented all of them. So we put not_implemented in the base class, but we can't find out if the function is implemented without calling it, unless we keep this information elsewhere. See vtkAlgorithm::ProcessRequest() for one way of dealing with this.) It's not okay to add not_implemented methods to the MATSHELL object since the solver may use a more restrictive API if the matrix doesn't support certain operations. Implementation could look something like this. class Matrix { public: // not_implemented for all matrix operations (mult, multtranspose, etc) protected: Mat A_petsc; Epetra::RowMatrix A_epetra; }; class PetscMatrix : public Matrix { // implements matrix operations in terms of operations on A_petsc. }; class ShellMatrix : public Matrix { public: // ... set_implemented(const char *methods, int count) { // call MatShellSetOperation(A_petsc,...) to register methods // call Epetra functions for A_epetra } // Note: inherits not_implemented for mult, multtranspose, etc }; // In user code ... class MyShell : public ShellMatrix { public: MyShell(..) { static const char *methods[] = {"mult", "multtranspose", "multadd", "getdiagonal"}; // or an enum set_implemented(implemented, 4); // ... } // ... implement mult, multtranspose, multadd, getdiagonal private: // For Tim's example Matrix A; NumericVector u,v,Au; }; Jed
pgpIV4hT1MPjg.pgp
Description: PGP signature
------------------------------------------------------------------------- This SF.Net email is sponsored by the Moblin Your Move Developer's challenge Build the coolest Linux based applications with Moblin SDK & win great prizes Grand prize is a trip for two to an Open Source event anywhere in the world http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________ Libmesh-devel mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-devel
