<div dir="ltr"><div class="gmail_default" style="font-size:small">​Thanks Matt and Barry,</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Matt - if I do the calculation in FormJacobian(), which makes by far the most sense and is as per your suggestion, do I need to set the operators of the SNES's KSP back to whatever they were before I set them?  <span style="background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">The linear system I want to solve within FormJacobian() involves the Jacobian matrix itself, and I want to remove the "nullspace" from that same matrix within FormFunction().</span></div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Barry - I'm trying to implement a homotopy solver.  In short, I have a system of n nonlinear equations in n variables, F(x), which is hard to solve because the Jacobian tends to become singular using standard methods​.  I want to add an auxiliary variable, lambda, to create a homotopy: H(lambda,x) = lambda*F(x) + (1-lambda)G(x), where G is "easy to solve", and the idea is that the Jacobian of the n+1 variable system will not become singular along the curve H(lambda,x) = 0.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">The method involves adding an equation to H so that the Jacobian H' is square.  The "submatrix" refers to the n x (n+1) matrix which represents the Jacobian without the added equation, whereas my FormJacobian() forms the entire (n+1) x (n+1) matrix H'.  I only refer to the submatrix because it has a nullspace, and I want to find it by solving a linear system designed for this purpose, H' u = b, where b is not the zero vector.  H' has no nullspace, but I want to remove the projection of u from my SNES solution vector, as u IS in the nullspace of the submatrix.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">The RHS vector b can be calculated outside the SNES solve.  I guess my FormJacobian() should look like this:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">FormJacobian(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx) {</div><div class="gmail_default" style="font-size:small">  KSP ksp;</div><div class="gmail_default" style="font-size:small">  MatNullSpace unull;</div><div class="gmail_default" style="font-size:small">  user_struct *user = (user_struct*)ctx;</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">  calculate Amat</div><div class="gmail_default" style="font-size:small">  assemble Amat</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">  if (Pmat != Amat) {</div><div class="gmail_default" style="font-size:small">    assemble Amat</div><div class="gmail_default" style="font-size:small">  }  </div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">  SNESGetKSP(snes,&ksp);</div><div class="gmail_default" style="font-size:small">  KSPSetOperators(ksp,Amat,Pmat);</div><div class="gmail_default" style="font-size:small">  KSPSolve(ksp,user->b,user->u);</div><div class="gmail_default" style="font-size:small">  </div><div class="gmail_default" style="font-size:small">  MatNullSpaceCreate(PetscObjectComm((PetscObject)Amat), PETSC_FALSE, 1, &(user->u),&unull);</div><div class="gmail_default" style="font-size:small">  MatSetNullSpace(Amat,unull);</div><div class="gmail_default" style="font-size:small">  MatNullSpaceDestroy(&unull);</div><div class="gmail_default" style="font-size:small">}</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Does this look right?</div><div class="gmail_default" style="font-size:small"><br></div></div>