<div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Mon, Jul 16, 2018 at 11:33 PM zakaryah <<a href="mailto:zakaryah@gmail.com">zakaryah@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div 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 style="font-size:small"><br></div><div 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 style="font-size:small"><br></div><div style="font-size:small">J(z) dz = -F(z)</div><div style="font-size:small"><br></div><div 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 style="font-size:small"><br></div><div style="font-size:small">A^T = [J^T e_k]</div><div style="font-size:small"><br></div><div style="font-size:small">and a trivial RHS:</div><div style="font-size:small"><br></div><div style="font-size:small">b^T = [0 ... 0 b_k]</div><div style="font-size:small"><br></div><div 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 style="font-size:small"><br></div><div style="font-size:small">A u = b (1)</div><div style="font-size:small"><br></div><div style="font-size:small">and a specific solution for the Newton iterate solves</div><div style="font-size:small"><br></div><div style="font-size:small">A v = c (2)</div><div style="font-size:small"><br></div><div style="font-size:small">where c^T = [-F^T c_{n+1}] </div><div style="font-size:small"><br></div><div 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 style="font-size:small"><br></div><div 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 style="font-size:small"><br></div><div 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></blockquote><div><br></div><div>Just calculate the nullspace inside FormJacobian().</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div 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></blockquote><div><br></div><div>Not really, unless you use some type of factorization for the preconditioner.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div 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>
</blockquote></div><br clear="all"><div>Yes, I think it should be fine.</div><div><br></div><div> Matt</div><div><br></div>-- <br><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.caam.rice.edu/~mk51/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div>