# [petsc-users] Implementing a homotopy solver

zakaryah zakaryah at gmail.com
Fri Jul 6 10:48:52 CDT 2018

```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.

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