trying to keep this on -devel... On Wed 2008-10-22 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
finite-dimensional 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 (matrix-size) 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
PETSc-native 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 PETSc-native 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
pgpA4nt1tZCg5.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
