From: Jed Brown <jed@59...>  20081020 18:42:42

On Mon 20081020 13:06, Roy Stogner wrote: > > > On Mon, 20 Oct 2008, Jed Brown wrote: >> >> 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. I think we have different ideas of what a shell matrix is. I say it just means that the user manages the required data structures. The essential operation (required to be a matrix) is mult, but there may be many others. > > It can't set individual matrix entries without adding some ugly data > structure on top of the existing userprovided 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? The user (implementing ShellMatrix) has some data structure, the library just doesn't get to see it. Presumably the user won't implement arbitrary insertion if their data structure doesn't support it. >> Certainly a PETSc MATSHELL is a firstclass matrix. All these >> operations are implementable by a shell matrix. >> >> http://wwwunix.mcs.anl.gov/petsc/petscas/snapshots/petsccurrent/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? It doesn't make sense for Tim's matrix, but a related factorization (factoring A in B = A + u v'). But LU factorization might make sense for someone else's shell matrix. >> 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. Why is it hard? They implement mult and whatever else, the returned PETSc matrix could have type MATSHELL. (Of course, if you don't implement MATOP_LUFACTOR_* you won't be able to use pc_type lu but the framework is general.) Of course converting the matrix to MPIAIJ or whatever might be hard and likely gets deeper than you'd like. You can always make a PCSHELL which is a proxy for the Laspack/Trilinos preconditioner. >> 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. I guess I don't see the necessity of distinguishing at the typelevel between matrices that have explicit entries and matrices that don't. As long as the appropriate methods are implemented in whatever you end up with at runtime, it's all the same. I don't see why such switch statements would ever be needed. >> 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 nonNULL) 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". Right, however at this time, the matrix hasn't been set up yet so the functions would need to throw a different exception early if they are implemented but not set up (to avoid possible segfault). Jed 