Hi Matt - thanks for all your help. Let's say I want exactly the same solver for the tangent vector and the SNES update, so I should reuse the KSP.
My attempt to do this looks like the summary of FormJacobian() in the previous message: - assemble Jacobian - assemble RHS - get KSP from the SNES passed to FormJacobian() - KSPSetFromOptions() - KSPSetOperators() - KSPSolve() I'm not sure that's the right approach, but it doesn't work - the KSPSolve() in the summary, i.e. the one for the tangent vector, seems to work fine. But the next KSPSolve() looks strange - it seems to use a preconditioner even with -pc_type none, etc. This makes me think I am doing something seriously wrong. On Sat, Sep 29, 2018, 8:16 AM Matthew Knepley <[email protected]> wrote: > On Fri, Sep 28, 2018 at 11:13 PM zakaryah <[email protected]> wrote: > >> 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. >> >> 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. >> >> Here's a summary of the code within FormJacobian(): >> >> - Set values and assemble Jacobian matrix A - this is working fine >> - Set values and assemble RHS vector b for linear system defining >> tangent vector n - this is working fine >> - 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 >> - 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 >> - MatSetNullSpace(A,NULL) - remove any existing null space from >> Jacobian >> - KSPSetOperators(my_ksp,A,P) - P is the other matrix in >> FormJacobian() >> - VecSet(n,0) - set initial guess to zero >> - 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 >> - VecNormalize(n,NULL) >> - >> >> MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,1,&n,&nullsp) >> - MatSetNullSpace(A,nullsp) >> - MatNullSpaceDestroy(&nullsp) >> - return >> >> 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. >> > > We need to clear up the usage first. If you want EXACTLY the same solver > for both solvers, then reuse > the KSP, otherwise do not do it. Does it work then? > > Thanks, > > Matt > > >> I hope someone can help me see what I'm doing wrong. >> >> On Sun, Jul 22, 2018 at 9:09 PM zakaryah <[email protected]> wrote: >> >>> Thanks Matt and Barry, >>> >>> 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? 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(). >>> >>> 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. >>> >>> 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. >>> >>> The RHS vector b can be calculated outside the SNES solve. I guess my >>> FormJacobian() should look like this: >>> >>> FormJacobian(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx) { >>> KSP ksp; >>> MatNullSpace unull; >>> user_struct *user = (user_struct*)ctx; >>> >>> calculate Amat >>> assemble Amat >>> >>> if (Pmat != Amat) { >>> assemble Amat >>> } >>> >>> SNESGetKSP(snes,&ksp); >>> KSPSetOperators(ksp,Amat,Pmat); >>> KSPSolve(ksp,user->b,user->u); >>> >>> MatNullSpaceCreate(PetscObjectComm((PetscObject)Amat), PETSC_FALSE, 1, >>> &(user->u),&unull); >>> MatSetNullSpace(Amat,unull); >>> MatNullSpaceDestroy(&unull); >>> } >>> >>> Does this look right? >>> >>> > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener > > https://www.cse.buffalo.edu/~knepley/ > <http://www.cse.buffalo.edu/~knepley/> >
