[petsc-users] Implementing a homotopy solver

[Smith, Barry F.]

   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!
>