[petsc-users] Implementing the Sherman Morisson formula (low rank update) in petsc4py and FEniCS?

[Olek Niewiarowski]

Hello,
I am a FEniCS user but new to petsc4py. I am trying to modify the KSP solver through the SNES object to implement the Sherman-Morrison formula (e.g.  http://fourier.eng.hmc.edu/e176/lectures/algebra/node6.html ). I am solving a nonlinear system of the form [K(u)−kaaT]Δu=−F(u). Here the jacobian matrix K is modified by the term kaaT, where k is a scalar.  Notably, K is sparse, while the term kaaT results in a full matrix. This problem can be solved efficiently using the Sherman-Morrison formula :

[K−kaaT]-1 = K-1  - (kK-1 aaTK-1)/(1+kaTK-1a)
I have managed to successfully implement this at the linear solve level (by modifying the KSP solver) inside a custom Newton solver in python by following an incomplete tutorial at https://www.firedrakeproject.org/petsc-interface.html#defining-a-preconditioner :
  *               while (norm(delU) > alpha):  # while not converged
  *
  *                   self.update_F()  # call to method to update r.h.s form
  *                   self.update_K()  # call to update the jacobian form
  *                   K = assemble(self.K)  # assemble the jacobian matrix
  *                   F = assemble(self.F)  # assemble the r.h.s vector
  *                   a = assemble(self.a_form)  # assemble the a_form (see Sherman Morrison formula)
  *
  *                   for bc in self.mem.bc:  # apply boundary conditions
  *                       bc.apply(K, F)
  *                       bc.apply(K, a)
  *
  *                   B = PETSc.Mat().create()
  *
  *                   # Assemble the bilinear form that defines A and get the concrete
  *                   # PETSc matrix
  *                   A = as_backend_type(K).mat()  # get the PETSc objects for K and a
  *                   u = as_backend_type(a).vec()
  *
  *                   # Build the matrix "context"  # see firedrake docs
  *                   Bctx = MatrixFreeB(A, u, u, self.k)
  *
  *                   # Set up B
  *                   # B is the same size as A
  *                   B.setSizes(*A.getSizes())
  *                   B.setType(B.Type.PYTHON)
  *                   B.setPythonContext(Bctx)
  *                   B.setUp()
  *
  *
  *                   ksp = PETSc.KSP().create()   # create the KSP linear solver object
  *                   ksp.setOperators(B)
  *                   ksp.setUp()
  *                   pc = ksp.pc
  *                   pc.setType(pc.Type.PYTHON)
  *                   pc.setPythonContext(MatrixFreePC())
  *                   ksp.setFromOptions()
  *
  *                   solution = delU    # the incremental displacement at this iteration
  *
  *                   b = as_backend_type(-F).vec()
  *                   delu = solution.vector().vec()
  *
  *                   ksp.solve(b, delu)
  *
  *                   self.mem.u.vector().axpy(0.25, self.delU.vector())  # poor man's linesearch
  *                   counter += 1
Here is the corresponding petsc4py code adapted from the firedrake docs:


  1.  class MatrixFreePC(object):
  2.
  3.      def setUp(self, pc):
  4.          B, P = pc.getOperators()
  5.          # extract the MatrixFreeB object from B
  6.          ctx = B.getPythonContext()
  7.          self.A = ctx.A
  8.          self.u = ctx.u
  9.          self.v = ctx.v
  10.         self.k = ctx.k
  11.         # Here we build the PC object that uses the concrete,
  12.         # assembled matrix A.  We will use this to apply the action
  13.         # of A^{-1}
  14.         self.pc = PETSc.PC().create()
  15.         self.pc.setOptionsPrefix("mf_")
  16.         self.pc.setOperators(self.A)
  17.         self.pc.setFromOptions()
  18.         # Since u and v do not change, we can build the denominator
  19.         # and the action of A^{-1} on u only once, in the setup
  20.         # phase.
  21.         tmp = self.A.createVecLeft()
  22.         self.pc.apply(self.u, tmp)
  23.         self._Ainvu = tmp
  24.         self._denom = 1 + self.k*self.v.dot(self._Ainvu)
  25.
  26.     def apply(self, pc, x, y):
  27.         # y <- A^{-1}x
  28.         self.pc.apply(x, y)
  29.         # alpha <- (v^T A^{-1} x) / (1 + v^T A^{-1} u)
  30.         alpha = (self.k*self.v.dot(y)) / self._denom
  31.         # y <- y - alpha * A^{-1}u
  32.         y.axpy(-alpha, self._Ainvu)
  33.
  34.
  35. class MatrixFreeB(object):
  36.
  37.     def __init__(self, A, u, v, k):
  38.         self.A = A
  39.         self.u = u
  40.         self.v = v
  41.         self.k = k
  42.
  43.     def mult(self, mat, x, y):
  44.         # y <- A x
  45.         self.A.mult(x, y)
  46.
  47.         # alpha <- v^T x
  48.         alpha = self.v.dot(x)
  49.
  50.         # y <- y + alpha*u
  51.         y.axpy(alpha, self.u)

However, this approach is not efficient as it requires many iterations due to the Newton step being fixed, so I would like to implement it using SNES and use line search. Unfortunately, I have not been able to find any documentation/tutorial on how to do so. Provided I have the FEniCS forms for F, K, and a, I'd like to do something along the lines of:
solver  = PETScSNESSolver() # the FEniCS SNES wrapper
snes = solver.snes()  # the petsc4py SNES object
## ??
ksp = snes.getKSP()
 # set ksp option similar to above
solver.solve()

I would be very grateful if anyone could could help or point me to a reference or demo that does something similar (or maybe a completely different way of solving the problem!).
Many thanks in advance!
Alex
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200204/30cb6f80/attachment-0001.html>