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<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/<http://www.math.tamu.edu/%7Ebangerth/>
>
> ______________________________**_________________
> dealii mailing list
> http://poisson.dealii.org/**mailman/listinfo/dealii<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
