Sorry again about the formatting - I struggle with notation in email. I think I see what you are saying - my augmented matrix should differ in the last row between the one used to calculate the nullspace, u, and the one the SNES uses as the Jacobian. If I construct the SNES Jacobian A so that u is truly in its nullspace, then I can just add u to the nullspace of A and let PETSc take care of the projection? I will try this - thanks for the excellent suggestion.
On Fri, Oct 5, 2018 at 4:42 PM Matthew Knepley <[email protected]> wrote: > On Fri, Oct 5, 2018 at 4:15 PM zakaryah <[email protected]> wrote: > >> 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? >> > > I do not claim to fully understand the above (it would be easier for me to > write the equations explicitly on their own lines). However, the issue with > nullspaces is soluble I think. So you say that > > J \in R^{m \cross (m+1)} is full rank > > so that its nullspace has dimension 1 and is spanned by u, > > J u = 0. > > Now you augment J to get a square matrix A by adding a row, > > / J \ u = / J u \ = / 0 \ > \ v^T / \ v^T u / \ 0 / > > if v is orthogonal to u. Thus, just choose v from the m-1 dimensional > space orthogonal to u. > > Thanks, > > Matt > > >> 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/> >> >> >> > > -- > 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/> >
