[petsc-users] Example of SNES using matrix free jacobian

[Daniel Pino Munoz]

Dear all,

I have a linear problem that I am currently solving with a KSP matrix-free.

I would like to move on to a non linear problem, so figure I could start 
by solving the same linear problem using SNES. So I am setting the 
problem as follows:

SNESCreate(PETSC_COMM_WORLD, &snes);
MatCreateShell(PETSC_COMM_WORLD, n_dofs, n_dofs, PETSC_DETERMINE, 
PETSC_DETERMINE, &ctx, &J);
MatCreateShell(PETSC_COMM_WORLD, n_dofs, n_dofs, PETSC_DETERMINE, 
PETSC_DETERMINE, &ctx, &B);
MatShellSetOperation(J, MATOP_MULT, (void (*)(void))(Multiplication));
MatCreateVecs(J, &x_sol, &b);
VecDuplicate(x_sol, &r);
SNESSetFromOptions(snes);
SNESSetFunction(snes, r, &(computeResidual), &ctx);
SNESSetUseMatrixFree(snes, PETSC_FALSE, PETSC_TRUE);
SNESGetLineSearch(snes, &linesearch);
SNESGetKSP(snes, &ksp);
KSPSetOperators(ksp, J, J);
KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);

I tested it with a small problem (compiled in debug) and it works.

When I compiled it in Release, it crashes with a segfault. I tried 
running the Debug version through valgrind, but even for a small 
problem, it is too slow. So I was wondering if you guys could see any 
rocky mistake on the lines I used above?

Otherwise, is there any example that uses a SNES combined with a matrix 
free KSP operator?

Thank you,

   Daniel