Doubts about SNESMF and matrix-free preconditioners

[Rafael Santos Coelho]

Hi, everyone (:

For a considerable amount of time, I've used PETSc a lot in my numerical
experiments, but then I stopped for a while and got rusty at it. So now, I'm
trying to get back on track by implementing my own serial matrix-free
preconditioner for the JFNK (Jacobian-free Newton-Krylov) method. Basicly,
it's a uniparametric LU-SGS (LU Symmetric Gauss-Seidel) preconditioner as
bellow:


   1. Let J(x), F(x) and di be the jacobian matrix, the nonlinear function
   and ith diagonal element of J(x), respectively
   2. Suppose J(x) = L + D + U and consider J(x)s = -F(x), where 's' is the
   Newton correction
   3. Let M = (D + wL)D^(-1)(D + wU), where w >0
   4. Now, we have M^(-1)J(x)s = -M^(-1)F(x), knowing that M^(-1)J(x)s =
   (F(x + hM^(-1)s) - F(x)) / h
   5. Let y = M^(-1)s = (D + wU)^(-1)D(D + wL)^(-1)s and z = (D + wL)^(-1)s
   <=> (D + wL)z = s
   6. To find z, we follow steps 7 and 8
   7. x' = x, z1 = s1 / d1,  x'1 = x1 + hz1
   8. for(i = 2; i <= n; i++) { zi = (si - w(Fi(x') - Fi(x)) / h) / di ; x'i
   = xi + hzi ; }
   9. Now, we have y = (D + wU)^(-1)Dz
   10. Let v = Dz and we got y = (D + wU)^(-1)v
   11. Let u = (D + wU)^(-1)v <=> (D + wU)u = v, and to find u we go through
   steps 12 and 13
   12. x' = x, un = vn / dn,  x'n = xn + hvn
   13. for(i = n - 1; i >= 1; i--) { ui = (vi - w(Fi(x') - Fi(x)) / h) / di
   ; x'i = xi + hui ; }
   14. Finally, we compute M^(-1)J(x)s = (F(x + hy) - F(x)) / h

As for the right-hand side, a similar approach can be done. I only need to
store diag[J(x)]. Well, I'm not sure how to code the steps above. However,
looking at "src/snes/examples/tutorials/ex6.c", things got cleared out a
little bit for me. I think it suffices to type "-snes_mf" on the
command-line, to write my own "MatrixFreePreconditioner" with all the steps
above and then

   1. ierr = PCSetType(pc, PCSHELL); CHKERRQ(ierr);
   2. ierr = PCShellSetApply(pc, MatrixFreePreconditioner); CHKERRQ(ierr);

Is my reasoning correct?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20081118/9801bba7/attachment.htm>