<div dir="ltr"><div class="gmail_default" style="font-size:small">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.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">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:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">y = (I-(x_1 e^T)/(x_1^T e + 1))x_2</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">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.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">I suppose I also need a shell preconditioner, which preconditions only M.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Thanks in advance for any help you can provide.</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Jul 3, 2018 at 11:28 AM, zakaryah <span dir="ltr"><<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_default" style="font-size:small">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.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">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:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">A = M + L</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">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.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">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?</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Thanks for your help!</div></div>
</blockquote></div><br></div>