[petsc-users] matrix-free methods

[coco at dmi.unict.it]

Dear all,

I am trying to use a matrix-free method to solve a linear system with  
Petsc. In particular, I noticed that if I solve the same linear system  
with the same KSP solver, the same preconditioner, and the same  
options, I obtain different convergence results if I solve the system  
with the full- or free-matrix method. Is that plausible?

Here is a piece of code:

//Method #1:
ierr=MatCreate(PETSC_COMM_WORLD,&M); CHKERRQ(ierr);
ierr=MatSetType(M, MATMPIAIJ); CHKERRQ(ierr);
ierr=MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,M); CHKERRQ(ierr);
[...filling the matrix...]
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,M,M,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr =  
KSPSetTolerances(ksp,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
ierr = KSPSolve(ksp,RHS,U);CHKERRQ(ierr);

//Method #2:
extern PetscErrorCode UserMult(Mat,Vec,Vec);
ierr=MatCreateShell(PETSC_COMM_WORLD,N,M,n,M,ctx,&M1); CHKERRQ(ierr);
ierr=MatShellSetOperation(M,MATOP_MULT,(void(*)(void)) UserMult);  
CHKERRQ(ierr);
ierr=MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr=MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp1);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp1,M1,M1,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr =  
KSPSetTolerances(ksp1,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
ierr = KSPSolve(ksp1,RHS,U1);CHKERRQ(ierr);

In particular, the method #1 converges in 15 iterations, while the  
method #2 in more than 1000.

Thank you in advance.
Armando