Everyone knows that when you code Newtons method you should compute the function and its Jacobian at the same time. How can one do this in PETSc?
The update in Newton’s method is computed as u^{n+1} = u^n – lambda * approx-inverse[J(u^n)] * F(u^n)]. The reason PETSc doesn’t default to computing both the function and Jacobian at the same time is • In order to do the line search, F (u^n – lambda * step) may need to be computed for several lambda, the Jacobian is not needed for each of those and one does not know in advance which will be the final lambda until after the function value is computed, so many extra Jacobians may be computed. • In the final step if || F(u^p)|| satisfies the convergence criteria then a Jacobian need not be computed. You are free to have your “FormFunction” compute as much of the Jacobian at that point as you like, keep the information in the user context (the final argument to FormFunction and FormJacobian) and then retreive the information in your FormJacobian() function.
