Hi, everyone (:<br><br>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:<br>
<br><ol><li>Let J(x), F(x) and di be the jacobian matrix, the nonlinear function and ith diagonal element of J(x), respectively<br></li><li>Suppose J(x) = L + D + U and consider J(x)s = -F(x), where 's' is the Newton correction<br>
</li><li>Let M = (D + wL)D^(-1)(D + wU), where w >0<br></li><li>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<br></li><li>Let y = M^(-1)s = (D + wU)^(-1)D(D + wL)^(-1)s and z = (D + wL)^(-1)s <=> (D + wL)z = s</li>
<li>To find z, we follow steps 7 and 8</li><li>x' = x, z1 = s1 / d1, x'1 = x1 + hz1</li><li>for(i = 2; i <= n; i++) { zi = (si - w(Fi(x') - Fi(x)) / h) / di ; x'i = xi + hzi ; }</li><li>Now, we have y = (D + wU)^(-1)Dz</li>
<li>Let v = Dz and we got y = (D + wU)^(-1)v</li><li>Let u = (D + wU)^(-1)v <=> (D + wU)u = v, and to find u we go through steps 12 and 13</li><li>x' = x, un = vn / dn, x'n = xn + hvn</li><li>for(i = n - 1; i >= 1; i--) { ui = (vi - w(Fi(x') - Fi(x)) / h) / di ; x'i = xi + hui ; }</li>
<li>Finally, we compute M^(-1)J(x)s = (F(x + hy) - F(x)) / h</li></ol>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<br>
<ol><li>ierr = PCSetType(pc, PCSHELL); CHKERRQ(ierr);</li><li>ierr = PCShellSetApply(pc, MatrixFreePreconditioner); CHKERRQ(ierr);<br></li></ol>Is my reasoning correct?<br><br>