thanks a lot for your answer. But defining the matrix_free_mult function
static, we can not call the virtual function vmult from inside which we
want to implement only in derived classes ...
I think this is a design flaw in PETSc: the function that is called is
only given a handle to the matrix and vectors but no other information.
You need this other information to identify which object this function
should belong to. I sort of saw this coming when you first proposed this
(because the MatShellSetOperation function, unlike similar functions in
POSIX, pthreads, etc that establish call-backs, does not allow you to
attach a void* that can be used to pass the information about the object
back to the callback function). But I wanted to see how you solve this
problem because I had no good idea initially either ;-)
Anyway, what you need is a way for the matrix_mult function (whether
it's global or a static member function) to establish what the proper
MatrixFree object is that belongs to the matrix object it is called
with. Here is a crude approach:
class MatrixFree {
virtual void vmult (...) = 0;
};
std::map<Mat, MatrixFree*> petsc_mat_to_deal_II_matrix_free_map;
void matrix_mult (Mat m, Vec src, Vec dst)
{
Assert (petsc_mat_to_deal_II_matrix_free_map.find(m)
!= petsc_mat_to_deal_II_matrix_free_map.end(),
ExcMessage ("Can't find corresponding deal.II object "
"for matrix!"));
petsc_mat_to_deal_II_matrix_free_map[m].vmult (dst, src);
}
void MatrixFree::do_reinit (const unsigned int m,
const unsigned int n,
const unsigned int local_rows,
const unsigned int local_columns)
{
...
MatCreateShell(communicator, local_rows, local_columns, m, n,
(void*)(NULL), &matrix);
MatShellSetOperation (matrix, MATOP_MULT,
&dealii::PETScWrappers::MPI::MatrixFree::matrix_free_mult_to_mult);
// now also register for the callback which MatrixFree object
// belongs to this PETSc Mat object
Assert (petsc_mat_to_deal_II_matrix_free_map.find(m)
== petsc_mat_to_deal_II_matrix_free_map.end(),
ExcMessage ("Matrix is already registered!"));
petsc_mat_to_deal_II_matrix_free_map[m] = matrix;
}
You'll have to remove the element from this map again when you destroy
the shell object, and you need to guard all operations on this map with
a mutex for it to be threadsafe (I can help with that, so don't worry
about it for now). Obviously, one would ultimately want the map and the
callback functions to be static members of the class, but that's just
syntax.
But maybe there is an easier approach for the library. I only want to
use PETSc here because there are so many linear solvers implemented
which work in parallel. But if we have the MatrixFree class and also the
parallel::distributed::Vector class in the library now, this is a bit
artificial, it should be more efficient and easier to modify the linear
solvers in the library to work with parallel::distributed::Vector's by
just adding an instantiation
template class Solver<parallel::distributed::Vector<double> >
That would be independently useful. Feel free to propose a patch.
Best
W.
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
