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 <[email protected]> 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 <[email protected]> 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!
>
