[petsc-users] Implementing a homotopy solver

[Smith, Barry F.]

   Since it is simpler to try first without the Sherman-Morrison formula, just try it. If the convergence is slow with A then you can consider adding it.

    Barry


> On Jul 6, 2018, at 11:30 AM, zakaryah <zakaryah at gmail.com> wrote:
> 
> ​Thanks for your help, Barry.
> 
> I agree about the preconditioning.  I still don't understand why I don't need a particular solver for my shell matrix.  My reasoning is that KSP is easy with M but difficult with A, since A has a dense row and column, whereas M is entirely sparse.  Sherman-Morrison seems to be an efficient way of dealing with this but I could be wrong.​
> 
> On Fri, Jul 6, 2018 at 12:19 PM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
> 
>    I am assuming the matrix M is provided as a regular PETSc AIJ (for example) sparse matrix.
> 
>     I am not sure you need bother with the Sherman-Morrison formula. Just write a shell matrix that performs the MatMult with A = M + L and call KSPSetOperators(ksp,A,M); PETSc will solve the linear system with A and construct the preconditioner using the matrix M; since A is only a rank one perturbation of M; the preconditioner constructed from M will be a good preconditioner for A.
> 
>    Barry
> 
> 
> > On Jul 6, 2018, at 10:48 AM, zakaryah <zakaryah at gmail.com> wrote:
> > 
> > I'm still hoping to get some advice on this.  I've done some more reading and it seems clear that I'll need to implement a shell matrix for A, so that the user context consists of the sparse matrix M and the rank one matrix L, which can be represented with a vector, u: L = u e^T, where e is the appropriate basis vector.  The part I'm most unsure of is how to implement the solver for linear systems involving the matrix A.
> > 
> > The solver for Ay=b solves two linear systems with the sparse matrix: Mx_1=u and Mx_2=b, using standard PETSc KSP routines for the sparse matrix M.  Then the Sherman-Morrison formula is used to calculate y:
> > 
> > y = (I-(x_1 e^T)/(x_1^T e + 1))x_2
> > 
> > My question is how to implement this specialized solver for the shell matrix A.  Do I use MATOP_SOLVE?  MATOP_MAT_SOLVE?  How do I know which MATOPs need to be provided?  Is there any documentation or tutorial for this type of problem, or for shell matrices in general?  I've only seen examples which provide MATOP_MULT.
> > 
> > I suppose I also need a shell preconditioner, which preconditions only M.
> > 
> > Thanks in advance for any help you can provide.
> > 
> > On Tue, Jul 3, 2018 at 11:28 AM, zakaryah <zakaryah at gmail.com> wrote:
> > I'm hoping to implement a homotopy solver for sparse systems in PETSc.  I have a question about the details of implementing the linear algebra steps so that I can take advantage of all the SNES tools.
> > 
> > My question doesn't have much to do with homotopy maps, per se.  The idea, which as far as I know comes from Layne Watson's 1986 paper, is to decompose the Jacobian of the homotopy map, A, into a sum of two matrices with special properties:
> > 
> > A = M + L
> > 
> > where M is sparse, symmetric, and invertible, and L is rank one.  Therefore, linear systems with M should be relatively easy to solve, using preconditioning and Krylov subspace methods.  The Newton update, which solves Az = b, can be found using the Sherman-Morrison formula.
> > 
> > I have two questions.  First, is it possible to implement this using tools that already exist in PETSc?  If not, is the best approach to write a shell preconditioner?  Second, would a homotopy solver like this be useful to the community?
> > 
> > Thanks for your help!
> > 
> 
>