Dear Zhengyong,
I have some problems with this matrix-free PETSc objects. As Wolfgang
suggested, I have implemented a MatrixFree base class in the following way:
namespace PETScWrappers
{
namespace MPI
{
class MatrixFree : public MatrixBase
{
public:
MatrixFree ();
/* some other constructors and functions */
virtual
void vmult (VectorBase &dst,
const VectorBase &src) const = 0;
private:
MPI_Comm communicator;
int matrix_free_mult (Mat A, Vec src, Vec dst);
void do_reinit (const unsigned int m,
const unsigned int n,
const unsigned int local_rows,
const unsigned int local_columns);
/* some other things */
};
}
}
where the matrix_free_mult () is the function we give to the PETSc MatShell
object and that just calls vmult() which has to be implemented in derived
classes:
int MatrixFree::matrix_free_mult (Mat A, Vec src, Vec dst)
{
VectorBase x(&(*src)), y(&(*dst));
vmult (y, x);
return (0);
}
the do_reinit function is called everytime we want to reinit and follows your
code:
void MatrixFree::do_reinit (const unsigned int m,
const unsigned int n,
const unsigned int local_rows,
const unsigned int local_columns)
{
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
int ierr = MatDestroy (matrix);
#else
int ierr = MatDestroy (&matrix);
#endif
AssertThrow (ierr == 0, ExcPETScError(ierr));
ierr = MatCreateShell(communicator, local_rows, local_columns, m, n,
(void*)(NULL), &matrix);
AssertThrow (ierr == 0, ExcPETScError(ierr));
ierr = MatShellSetOperation (matrix, MATOP_MULT,
(void(*)(void))&dealii::PETScWrappers::MPI::MatrixFree::matrix_free_mult_to_mult);
AssertThrow (ierr == 0, ExcPETScError(ierr));
ierr = MatSetFromOptions (matrix);
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
(matrix is the Mat object in the MatrixBase class) but this does not work, gcc
4.6.3 compiles the code, but PETSc breaks down when calling the MatMult
function. If I define a multiplication
matrix_mult (Mat A, Vec src, Vec dst) { ... }
somewhere outside the class and set
ierr = MatShellSetOperation (matrix, MATOP_MULT,
(void(*)(void))matrix_mult);
everything works fine ...
The type for the matrix_operation in MatShellSetOperation function is void
(*f)(void) , but typecasting seems not to work for class functions here ...
maybe there is something completely wrong, any ideas? Sorry if this is a stupid
question.
Many thanks,
Martin
----- Original Message -----
From: RenZhengYong
To: [email protected]
Sent: Friday, June 08, 2012 4:51 PM
Subject: Re: [deal.II] MatrixFree with PETSc
Hi, Martin,
I used the matrix-free version of petsc.
Using iterative solvers, you should offer a subroutine to compute the product
of system matrix with a guess vector. This means you should offer your own
subroutine:
int your_own_mat_vec_product(Mat A, Vec X, Vec Y);
Then do a typedef:
typedef int (*USERMULT)(Mat A, Vec X, Vec Y);
where Mat and Vec are the matrix and vector in petsc.
Then, in you own subroutine which call petsc to solve your Ax=b, the order of
the calls should be like this:
void solve_Ax_b_by_matrix_free_petsc(int argc, char** argv, /*options for
pestc*/
Matrix P, /*a possible preconditioner*/
USERMULT usermult, /*a function pointer to your own mat_vec_product*/
Vector B, /*RHS vector*/,
Vector X. /*solution*/) /*Matrix, Vector are your own data structure*/
{
PetscMPIInt cpu_size;
PetscInitialize(&argc, &argv, (*char*)(0), (char*)(0));
MPI_Comm_size(PETSC_COMM_SELF, &cpu_size); /*working on a cpu*/
Mat _A, _P;
int m=P.rows();
/*Creat A in a matrix-free style*/
MatCreatShell( PETSC_COMM_SELF, m, m, m, m, (void*)(NULL), &_A);
MatShellSetOperation(_A, MATOP_MULT, (void(*)()_ usermult);
MatSetFromOptions(_A); /*must be called*/
/*translate P to _P which is in Mat format*/
Vec _X, _B;
VecCreate(PETSC_COMM_SLEF, &_X);
VecSetSizes(_X, PETSC_DECIDE, m);
VecSetFromOptions(_X); /*must be called*/
VecDuplicate(_X, &_B);
/*Creat Ksp and PC*/
KSP _ksp; PC _pc;
KSPCreate(PETSC_COMM_SELF, &_ksp);
KSPSetType(_ksp, KSPGMRES); /*default is GMRES solver*/
KSPSetOperators(_ksp, _A, _P, SAME_PRECONDITIONER); /*set preconditioner
_P; if user-defined preconditioner, set _P to _A */
KSPGetPC(_ksp, &_pc);
PCSetType(_pc, PCJACOBI); /*default is Jacobi*/
KSPSetTolerances(_ksp, 1.e-14, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
KSPSetFromOptions(_ksp); /*must be called*/
/*fullfill _X, _B with X and B*/
/*finish the assemble of X and B, no modification allowed*/
VecAssemblyBegin(_X);
VecAssemblyEnd (_X);
VecAssemblyBegin(_B);
VecAssemblyEnd (_B);
/*finish the assemble of _A, _P*/
MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd (_A,MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(_P,MAT_FINAL_ASSEMBLY);
MatAssemblyEnd (_P,MAT_FINAL_ASSEMBLY);
KSPSolve(_ksp, _B, _X);
/*destroy*/
VecDestroy(_X);
VecDestroy(_B);
MatDestroy(_A);
MatDestroy(_P);
KSPDestroy(_ksp); /*PC is also destroyed*/
PetscFinalize();
return;
}
Then, you call to solve_Ax_b_by_matrix_free_petsc in main() should be
int main(int argc, char** argv)
{
Matrix P;
Vector B, X;
solve_Ax_b_by_matrix_free_petsc(argc, argv, P, your_own_mat_vec_product, B,
X);
/*Then, you will obtain your solution in vector X;*/
return 0;
}
Good luck.
Zhengyong
On Fri, Jun 8, 2012 at 4:07 PM, Wolfgang Bangerth <[email protected]>
wrote:
Martin,
I am trying to use this very cool new MarixFree-technique together with
PETSc and distributed triangulations, or in other words, I want to
combine step-37 with step-40. But one of the problems is, that
PETScWrappers::SolverBase only accept PETScWrappers::MatrixBase as
matrix-typ and I see no other way, than to derive the LaplaceOperator
(the class that implements the matrix operations) from
PETScWrappers::MatrixBase and to implement all the methods needed, but
maybe I missed something and someone knows a much easier way ....??
I am not a PETSc expert but I imagine that PETSc has a way to create a
'Mat' handle in such a way that it represents a matrix-free object that can be
given to solvers. In fact, I just found how:
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateShell.html
The MatrixBase class stores the 'Mat' handle of such objects, and so I
would suggest that you derive a class that represents such objects in the
following way:
class PETScWrappers::MatrixFree : PETScWrappers::MatrixBase
{
// overload existing virtual functions
virtual
void vmult (Vector &dst, Vector &src) const = 0;
};
This class would set up things in such a way that it calls MatCreateShell
in the beginning and then registers some callback with MatShellSetOperation in
such a way that whenever the callback is called, it in turn calls the vmult()
function above. Then you can implement your LaplaceOperator class in terms of
this MatrixFree interface. Ideally, the vmult function is the only one that
you'd need to implement in a derived class.
Best
W.
PS: Of course we'd love to have this PETScWrappers::MatrixFree class if you
decide to go this route :-)
------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.tamu.edu/~bangerth/
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
--
Zhengyong Ren
AUG Group, Institute of Geophysics
Department of Geosciences, ETH Zurich
NO H 47 Sonneggstrasse 5
CH-8092, Zürich, Switzerland
Tel: +41 44 633 37561
e-mail: [email protected]
Gmail: [email protected]
------------------------------------------------------------------------------
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
