<div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Fri, Sep 28, 2018 at 11:13 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><div style="font-size:small">I'm working on a homotopy solver which follows a zero curve through solution space. A central aspect of this method is to calculate the vector tangent to the curve, predict the next point, then correct using iteration of, e.g. Newton's method, orthogonal to the tangent vector.</div><div style="font-size:small"><br></div><div style="font-size:small">Previously, we discussed the possibilities of implementing this within PETSc's SNES. Within my FormJacobian() function, I construct the linear system which defines the tangent vector, solve it, then add the vector to the nullspace of the Jacobian. I think that in principle this can work, but I suspect I'm doing something wrong.</div><div style="font-size:small"><br></div><div style="font-size:small">Here's a summary of the code within FormJacobian():</div><div><ul style="font-size:small"><li>Set values and assemble Jacobian matrix A - this is working fine</li><li>Set values and assemble RHS vector b for linear system defining tangent vector n - this is working fine</li><li>SNESGetKSP(snes,&my_ksp) - I thought it made sense to use the KSP associated with the SNES, hoping that PCs which use a factorization could be reused when the SNES calls KSPSolve() to calculate the update</li><li>KSPSetFromOptions(my_ksp) - not sure this is necessary but one of my problems is setting options for this KSP from the command line and even with this call it doesn't seem to be working properly</li><li>MatSetNullSpace(A,NULL) - remove any existing null space from Jacobian</li><li>KSPSetOperators(my_ksp,A,P) - P is the other matrix in FormJacobian()</li><li>VecSet(n,0) - set initial guess to zero</li><li>KSPSolve(my_ksp,b,n) - these solves appear to work, i.e. use the options passed from the command line with -ksp_XXX or -pc_XXX</li><li>VecNormalize(n,NULL)</li><li>MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,1,&n,&nullsp)</li><li>MatSetNullSpace(A,nullsp)</li><li>MatNullSpaceDestroy(&nullsp)</li><li>return</li></ul></div><div>The immediate problem is that the subsequent KSPSolve(), i.e. the one called internally by SNESSolve(), behaves strangely. For example, if I use -pc_type none -ksp_monitor -ksp_monitor_true_residual, then the KSPSolve() that I call within FormJacobian() looks correct - "preconditioned" norm and true norm are identical, and both converge as I expect (i.e. slowly but geometrically). However, the subsequent KSPSolve(), internal to the SNESSolve(), has large differences between the preconditioned norm and the true norm. In addition, the KSP does not converge in the true residual, but I'll have a hard time debugging that without knowing how to properly set the options.</div></div></div></blockquote><div><br></div><div>We need to clear up the usage first. If you want EXACTLY the same solver for both solvers, then reuse</div><div>the KSP, otherwise do not do it. Does it work then?</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</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><div>I hope someone can help me see what I'm doing wrong.</div></div></div><br><div class="gmail_quote"><div dir="ltr">On Sun, Jul 22, 2018 at 9:09 PM zakaryah <<a href="mailto:zakaryah@gmail.com" target="_blank">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">Thanks Matt and Barry,</div><div style="font-size:small"><br></div><div 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 style="font-size:small"><br></div><div 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 style="font-size:small"><br></div><div 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 style="font-size:small"><br></div><div 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 style="font-size:small"><br></div><div style="font-size:small">FormJacobian(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx) {</div><div style="font-size:small"> KSP ksp;</div><div style="font-size:small"> MatNullSpace unull;</div><div style="font-size:small"> user_struct *user = (user_struct*)ctx;</div><div style="font-size:small"><br></div><div style="font-size:small"> calculate Amat</div><div style="font-size:small"> assemble Amat</div><div style="font-size:small"><br></div><div style="font-size:small"> if (Pmat != Amat) {</div><div style="font-size:small"> assemble Amat</div><div style="font-size:small"> } </div><div style="font-size:small"><br></div><div style="font-size:small"> SNESGetKSP(snes,&ksp);</div><div style="font-size:small"> KSPSetOperators(ksp,Amat,Pmat);</div><div style="font-size:small"> KSPSolve(ksp,user->b,user->u);</div><div style="font-size:small"> </div><div style="font-size:small"> MatNullSpaceCreate(PetscObjectComm((PetscObject)Amat), PETSC_FALSE, 1, &(user->u),&unull);</div><div style="font-size:small"> MatSetNullSpace(Amat,unull);</div><div style="font-size:small"> MatNullSpaceDestroy(&unull);</div><div style="font-size:small">}</div><div style="font-size:small"><br></div><div style="font-size:small">Does this look right?</div><div style="font-size:small"><br></div></div>
</blockquote></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><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.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>