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