Wolfgang and Zhengyong, thanks a lot for your answers, that is exactly what I am looking for, but missed this matrix-free object of PETSc. I will put the functionality into a PETScWrappers::MatrixFree interface the way you both suggest.
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
