From: Jed Brown <jed@59...>  20081022 16:52:50

trying to keep this on devel... On Wed 20081022 12:38, Tim Kroeger wrote: > Your inheritance hierarchy makes user code that uses the shell matrices > functionality solver dependent: The user has to explicitly instantiate > e.g. PetscDyadMatrix. You might later implement shell matrices for other > solver packages, and I would find it more appropriate to have a > DyadMatrix class as the user API and let everything that depends on the > solver package be hidden from the user. There are two ways to think of what it means to be a PetscMatrix, (a) the internal storage is PETSc native and (b) a matrix which can be used with PETSc solvers. Because of shell matrices, these are not the same. The existing PetscMatrix plays both roles, but the only thing needed by the solver is the ability to get a PETSc Mat object. Suppose the user implements a ShellMatrix (call it UserShellMatrix) which implements some libmesh matrix operations (multiply, and maybe others). If we (the library) can determine what operations have been implemented (for a ShellMatrix or any other libmesh Matrix), producing a PETSc MatShell is entirely routine and certainly shouldn't be done by the user. Creating the MatShell does not require knowledge of anything internal to UserShellMatrix, it only needs to know what has been implemented. Once the PETSc Mat object has been made, there is no harm in wrapping it in a PetscMatrix. What about just adding a constructor (my abstract class Matrix is an arbitrary linear transformation on a finitedimensional vector space, i.e. it must implement mult, nothing else is required). PetscMatrix::PetscMatrix(Matrix *B) { MatCreateShell(comm,B.m,B.n,B.M,b.N,(void*)B,&_mat); MatShellSetOperation(_mat,MATOP_MULT,(void(*)(void))matshell_mult); // set wrappers for any other operations supported by B } // wrappers look like this and are not exported PetscErrorCode matshell_mult(Mat A,Vec x,Vec y) { ShellMatrix B; MatShellGetContext(A,(void**)&B); B.mult(x,y); } No implementation is exposed, there is no (matrixsize) duplicate storage, and you can use any PETSc solver (KSP) with any matrix type (shell or otherwise, provided they implement mult and possibly mult_transpose). For preconditioning with a shell matrix, the user will either provide a PCShell (could be done similarly) or a different PETScnative preconditioning matrix. The scary thing here would be the implicit cast (creating a new PetscMatrix which is a MatShell) occurs when a Matrix (which might even be a PetscMatrix) is passed to a Petsc solver. This is why I think it's cleaner to put the getter for the PETSc Mat in Matrix since this would enable the solvers to take any Matrix. In this case, we'd have Mat Matrix::petsc_mat() { Mat mat; MatCreateShell(comm,B.m,B.n,B.M,b.N,(void*)this,&mat); MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))matshell_mult); // set wrappers for any other operations supported by B return mat; } Mat PetscMatrix::petsc_mat() { return _mat; } This is great except for the (small) memory leak. We can wrap the Mat object in a reference counted pointer which calls MatDestroy() when needed or just put the _mat object in Matrix (which I don't see as a huge break in abstraction) in which it is guaranteed to be unique for any Matrix. With this design, every solver takes a plain Matrix (a PetscMatrix just means that it uses a PETScnative format). The inheritance is very simple and it's easy to add new native formats (with no modification it will automatically work with all existing solvers supporting the concept of shell matrices) and new solver packages (this requires one function be added to Matrix but I think it's simple enough and infrequent enough to be okay). Jed 