<div dir="ltr"><div class="gmail_default" style="font-size:small">I'm still working on this, and I have a separate question about the implementation, which is related to calculating and removing the nullspace from the Jacobian.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">I have a system of n nonlinear equations in n+1 variables, due to the auxiliary variable (homotopy variable).  Given a guess, z, I'd like to calculate a corrector (Newton iteration) dz which solves</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">J(z) dz = -F(z)</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">where J is the n x (n+1) Jacobian of F.  Since J has a rank-one nullspace, I'd like to calculate that nullspace in each iteration of the SNES, then find the minimal norm solution dz.  The tricks for finding the nullspace without messing up the nice structure of the n x n submatrix of J are in Watson's paper: add a "dummy" row to J to form a square matrix A:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">A^T = [J^T e_k]</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">and a trivial RHS:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">b^T = [0 ... 0 b_k]</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">where the value of b_k is essentially arbitrary but can be chosen to orient the solution in the direction of increasing arc length.  A vector in the rank-one nullspace of J, u, solves</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">A u = b              (1)</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">and a specific solution for the Newton iterate solves</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">A v = c              (2)</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">where c^T = [-F^T c_{n+1}] </div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">and c_{n+1} is also an arbitrary constant.  The desired Newton update dz is just the Graham-Schmidt orthogonalization of v wrt u.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">From the previous discussion I think A can be implemented as a shell matrix as a straightforward way to efficiently solve the linear systems.  I have three questions about the implementation of this within the SNES:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">1) What is the best way to "hook" a calculation of the nullspace of J, via eq. (1), within the SNES?  c is formed in FormFunction(), A is formed in FormJacobian(), but I need an intermediate step to form b and solve for u, then, I assume, set it as the nullspace of A (see question 3).</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">2) Is there a way to reuse the KSP for (1) to efficiently solve (2), given the method which answers question 1) ?</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">3) Can I set u to be in the nullspace of A?  It's really in the nullspace of the submatrix J, while the full matrix A has no nullspace, but the point is that I don't care about the last row of A - that's the whole trick of the homotopy.  I guess my question is whether anything in the solver will complain that Au != 0, because for all other purposes I want to project u out of the solution to (2) and it seems like treating it as a nullspace will achieve that.</div></div>