[petsc-users] Implementing a homotopy solver

[zakaryah]

​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!
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180706/15c1f681/attachment.html>