Dear Martin, Just one more comment (you probably have figured it out in the meantime): If you work in parallel and supply PETSc vectors to the MatrixFree/FEEvaluation concept, you will need to first copy them into a parallel::distributed::Vector with the layout given by the one in the MatrixFree class (there is no such function yet, but you could have a look at how things are done for the usual dealii::Vector<double>). This is because FEEvaluation reads from vector elements in local MPI index space, but all vectors except parallel::distributed::Vector<Number> are accessed with operator() which expects global indices.
Best regards, Katharina On 8 jun 2012, at 23.08, Martin Steigemann wrote: > 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
_______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
