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