On Mon, 20 Oct 2008, Jed Brown wrote:
>> 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? It can return individual matrix entries, but I'd want to put big red warning flags on the method that does so, since users of it would naturally be expecting O(1) and not O(N) cost. It can't set individual matrix entries without adding some ugly data structure on top of the existing user-provided shell function and making some assumptions about the user not changing their shell function "behind your back" (even though such changes would be implicit in most matrices for adaptive FEM). If the user just hands you a pointer to a function that scales a NumericVector, how do you get into that function to change individual entries properly afterward? > 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 That's interesting. How do you do an LU factorization of a shell matrix without just turning it into a sparse (or in Tim's case, dense) matrix first? > 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. This is actually part of why I'm pushing for "MatrixBase" - there are other common matrix types like his that could be put into the library without requiring a user function but which don't fall well into a SparseMatrix paradigm. > How about just adding SparseMatrix::petsc_mat() which returns a PETSc > Mat for any matrix type. That's going to be hard for the Laspack and Trilinos matrix types. ;-) But I get your point - if the underlying package already has polymorphism in it, we want their bottom layer encapsulated by our bottom layer to avoid reinventing their wheel. > 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. Awkward how? Is it really so bad to have a ShellMatrix that has constructors which take a shell function and a SparseMatrix that has constructors which take sparsity information? I think it would be more awkward if we ended up with more "am I really a sparse matrix or am I really a shell matrix" if/switch statements inside a hybrid class. > To be clear, I see ShellMatrix derived from Matrix and implementing > petsc_mat() (and similarly epetra_mat()) by calling > MatShellSetOperation(...) for all implemented operations. As long as the user never has to touch the internals, that makes sense. > 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. In many cases we could just try calling the method, catch a not_implemented exception, and fall back on other code if we get one. But that definitely doesn't qualify as "elegantly". Making users specify what methods they're overriding as in your example looks better. Not that it's soon enough to help us, but does anyone know if they're adding any better reflection capabilities to C++0x? --- Roy ------------------------------------------------------------------------- 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
