<div dir="ltr"><div class="gmail_default" style=""><div class="gmail_default" 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 class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" 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 class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Here's a summary of the code within FormJacobian():</div><div class="gmail_default" style=""><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 class="gmail_default" style="">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 class="gmail_default" style=""><br></div><div class="gmail_default" style="">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">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 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>
</blockquote></div>