Thanks Matt - I was trying to write it in linear algebra notation, but the nature of the problem might produce some confusion. I was also originally a bit confused by a fairly major typo in the paper I was following, but I'm almost sure I have it clear now. I'll try to be concise - the meat is in the last two paragraphs.
My system of nonlinear equations, F(x)=0, has dimension m. Traditional methods for solving, like Newton's method, tend to reach points at which the Jacobian becomes singular. I'm trying the homotopy trick, which involves introducing an additional scalar variable, s, and solving the system H(x,s) = sF(x) + (1-s)(x-x_0), starting at s=0 and following the zero curve to, hopefully, a solution at s=1. The homotopy map H only has m components, so its Jacobian J is m x (m+1), and has rank m. I need to find a vector u in its null space: J u = 0. Note that u is size m+1. Rather than solve J u = 0, I make an augmented matrix A, which is (m+1)x(m+1), and consists of J in its first m rows, then a somewhat arbitrary vector in the last row. Then the vector u is the solution to A u = b, where b is zero everywhere except its last element. This solve seems to be working fine. The SNES should then solve for a particular solution y of the Newton update, J y = - H(y). A (unique) particular solution can also be found by solving A y = [-H(y) 0]^T. Finally - and this is the step I can't figure out how to implement properly, I want to find the norm-minimizing solution z for the Newton update: z = y - u(y^T u)/(u^T u). I doubt there is any way to do this elegantly within the existing routines, because u is not in the nullspace of A. Rather, it's in the nullspace of J, which is the submatrix I care about. I think what I need is a way to just apply the orthogonal projection in the preceding paragraph, after the SNES solves the linear system but before the Newton update is applied. I don't know the SNES internals well at all. Is it possible to hook this in somehow, or put it in an existing function? Thanks again for all your help. On Oct 1, 2018 6:49 AM, "Matthew Knepley" <[email protected]> wrote: On Sun, Sep 30, 2018 at 6:06 PM zakaryah <[email protected]> wrote: > OK, thanks. > > I'm using a composite DM, DM_packer. To make a separate KSP, I did the > following in the FormJacobian() routine, after assembling the Jacobian > matrix A and the RHS for the tangent vector, b: > > - KSPCreate(PETSC_COMM_WORLD,&ksp) > - KSPSetDM(ksp,DM_packer) > - KSPSetDMActive(ksp,PETSC_FALSE) - because I want to set the operators > - KSPSetOperators(ksp,A,P) > - VecSet(n,0) - set initial guess to zero > - KSPSolve(ksp,b,n) - this solve works correctly > - VecNormalize(n,NULL) > - > > MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,1,&n,&nullsp) > - MatSetNullSpace(A,nullsp) > > Then, with -snes_type newtonls -pc_type none -ksp_monitor > -ksp_monitor_true_residual -ksp_view, the output for the KSPSolve described > above, i.e. for the tangent vector, looks correct: > > 0 KSP preconditioned resid norm 1.000000000000e+03 true resid norm > 1.000000000000e+03 ||r(i)||/||b|| 1.000000000000e+00 > ... > > 185 KSP preconditioned resid norm 9.900713131874e-03 true resid norm > 9.900713131904e-03 ||r(i)||/||b|| 9.900713131904e-06 > > Linear solve converged due to CONVERGED_RTOL iterations 185 > > KSP Object: 1 MPI processes > > type: gmres > > restart=30, using Classical (unmodified) Gram-Schmidt > Orthogonalization with no iterative refinement > > happy breakdown tolerance 1e-30 > > maximum iterations=10000, initial guess is zero > > tolerances: relative=1e-05, absolute=1e-50, divergence=10000. > > left preconditioning > > using PRECONDITIONED norm type for convergence test > > PC Object: 1 MPI processes > > type: none > > linear system matrix = precond matrix: > > Mat Object: 1 MPI processes > > type: seqaij > > rows=78247, cols=78247 > > total: nonzeros=6063481, allocated nonzeros=6063481 > > total number of mallocs used during MatSetValues calls =0 > > using I-node routines: found 26083 nodes, limit used is 5 > > > But the next KSP, i.e. the one in the SNES, doesn't converge in the true > residual: > > > 0 KSP preconditioned resid norm 2.045599092896e-04 true resid norm > 2.803828296212e-04 ||r(i)||/||b|| 1.000000000000e+00 > > 191 KSP preconditioned resid norm 2.009941278534e-09 true resid norm > 1.636010142734e-04 ||r(i)||/||b|| 5.834915586465e-01 > > > KSP Object: 1 MPI processes > > type: gmres > > restart=30, using Classical (unmodified) Gram-Schmidt > Orthogonalization with no iterative refinement > > happy breakdown tolerance 1e-30 > > maximum iterations=10000, initial guess is zero > > tolerances: relative=1e-05, absolute=1e-50, divergence=10000. > > left preconditioning > > using PRECONDITIONED norm type for convergence test > > PC Object: 1 MPI processes > > type: none > > linear system matrix = precond matrix: > > Mat Object: 1 MPI processes > > type: seqaij > > rows=78247, cols=78247 > > total: nonzeros=6063481, allocated nonzeros=6063481 > > total number of mallocs used during MatSetValues calls =0 > > has attached null space > > using I-node routines: found 26083 nodes, limit used is 5 > > > My guess about what's going on is that the tangent vector n isn't really > in the nullspace of A. > I do not understand the rest, but this is the problem. Maybe it would make more sense if you wrote things in linear algebraic notation. Thanks, Matt > Rather, it's in the nullspace of the m x (m+1) submatrix of A. So, An=c > e_{m+1}, where c is an arbitrary constant and e_{m+1} is the m+1 th basis > vector. The nonlinear function also has an added row, F_{m+1}, which is > set to zero in the FormFunction() routine. I don't care about the value of > F_{m+1}, but I suppose that if it's included in the true residual, and NOT > in the "preconditioned" residual, even with pc_type none, then it will be > tricky to diagnose the performance. Should I be using my own routine to > evaluate the residual, so that F_{m+1} is not included? > -- 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/>
