<div dir="ltr"><div dir="ltr">On Mon, Feb 10, 2020 at 9:09 AM Olek Niewiarowski <<a href="mailto:aan2@princeton.edu">aan2@princeton.edu</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr">
<div style="font-family:Calibri,Arial,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Barry, </div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Thank you for your help and detailed suggestions. I will try to implement what you proposed and will follow-up with any questions. In the meantime, I just want to make sure I understand the use of SNESSetPicard:</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<table style="font-family:Times;font-size:medium">
<tbody>
<tr>
<td width="40"></td>
<td align="LEFT"><b>r</b></td>
<td>- vector to store function value</td>
</tr>
<tr>
<td width="40"></td>
<td align="LEFT"><b>b</b></td>
<td>- function evaluation routine - <b>my F(u) function </b></td>
</tr>
<tr>
<td width="40"></td>
<td align="LEFT"><b>Amat</b></td>
<td>- matrix with which A(x) x - b(x) is to be computed -<b> a MatCreateLRC() -- what's the best way of passing in scalar k?</b></td>
</tr>
<tr>
<td width="40"></td>
<td align="LEFT"><b>Pmat</b></td>
<td>- matrix from which preconditioner is computed (usually the same as Amat) - <b>
a regular Mat()</b></td>
</tr>
<tr>
<td width="40"></td>
<td align="LEFT"><b>J</b></td>
<td>- function to compute matrix value, see<span> </span><a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESJacobianFunction.html#SNESJacobianFunction" target="_blank">SNESJacobianFunction</a><span> </span>for details on its calling sequence -- <b> computes <span style="font-family:Calibri,Arial,Helvetica,sans-serif;background-color:rgb(255,255,255);display:inline">K
+ kaa'</span> </b></td>
</tr>
</tbody>
</table>
<br>
</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
By the way, t<span style="font-family:Calibri,Arial,Helvetica,sans-serif;background-color:rgb(255,255,255);display:inline">he manual page states "we do not recommend using this routine. It is far better to provide the nonlinear function
F() and some approximation to the Jacobian and use an approximate Newton solver :-)"</span></div></div></blockquote><div><br></div><div>This is still correct. Barry is suggesting Picard since the implementation for you is a little bit simpler. However, the Picard solver</div><div>is missing the K' piece of the Jacobian, so it will only ever convergence linearly, as opposed to quadratically for Newton in the convergence</div><div>basin. However, once you have Picard going, it will be simple to switch to Newton by just providing that extra piece.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">
<div style="font-family:Calibri,Arial,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Thanks again!</div>
<div style="font-family:Calibri,Arial,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<div id="gmail-m_6035646582007503743Signature">
<div style="font-family:Tahoma;font-size:13px">
<div>
<div><font style="font-size:13px"><font><b style="color:rgb(34,34,34);font-size:small;background-color:rgb(255,255,255)">Alexander (Olek) Niewiarowski</b>
</font></font>
<div style="color:rgb(34,34,34);background-color:rgb(255,255,255)"><font color="#666666" face="Helvetica" size="1">PhD Candidate, Civil & Environmental Engineering</font></div>
<div style="color:rgb(34,34,34);background-color:rgb(255,255,255)"><font color="#666666" face="Helvetica" size="1">Princeton University, 2020</font></div>
<div style="color:rgb(34,34,34);background-color:rgb(255,255,255)"><font color="#666666" face="Helvetica" size="1">Cell: +1 (610) 393-2978</font></div>
</div>
</div>
</div>
</div>
<div id="gmail-m_6035646582007503743appendonsend"></div>
<hr style="display:inline-block;width:98%">
<div id="gmail-m_6035646582007503743divRplyFwdMsg" dir="ltr"><font face="Calibri, sans-serif" style="font-size:11pt" color="#000000"><b>From:</b> Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>><br>
<b>Sent:</b> Thursday, February 6, 2020 13:25<br>
<b>To:</b> Olek Niewiarowski <<a href="mailto:aan2@princeton.edu" target="_blank">aan2@princeton.edu</a>><br>
<b>Cc:</b> Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>>; <a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b> Re: [petsc-users] Implementing the Sherman Morisson formula (low rank update) in petsc4py and FEniCS?</font>
<div> </div>
</div>
<div><font size="2"><span style="font-size:11pt">
<div><br>
If I remember your problem is K(u) + kaa' = F(u)<br>
<br>
You should start by creating a SNES and calling SNESSetPicard. Read its manual page. Your matrix should be a MatCreateLRC() for the Mat argument to SNESSetPicard and the Peat should be just your K matrix.<br>
<br>
If you run with -ksp_fd_operator -pc_type lu will be using K to precondition K + kaa' + d F(U)/dU . Newton's method should converge at quadratic order. You can use -ksp_fd_operator -pc_type anything else to use an iterative linear solver as the preconditioner
of K. <br>
<br>
If you really want to use Sherman Morisson formula then you would create a PC shell and do<br>
<br>
typedef struct{<br>
KSP innerksp;<br>
Vec u_1,u_2;<br>
} YourStruct;<br>
<br>
SNESGetKSP(&ksp);<br>
KSPGetPC(&pc);<br>
PCSetType(pc,PCSHELL);<br>
PCShellSetApply(pc,YourPCApply)<br>
PetscMemclear(yourstruct,si<br>
PCShellSetContext(pc,&yourstruct);<br>
<br>
Then <br>
<br>
YourPCApply(PC pc,Vec in, Vec out)<br>
{<br>
YourStruct *yourstruct;<br>
<br>
PCShellGetContext(pc,(void**)&yourstruct)<br>
<br>
if (!yourstruct->ksp) {<br>
PCCreate(comm,&yourstruct->ksp);<br>
KSPSetPrefix(yourstruct->ksp,"yourpc_");<br>
Mat A,B;<br>
KSPGetOperators(ksp,&A,&B);<br>
KSPSetOperators(yourstruct->ksp,A,B);<br>
create work vectors <br>
}<br>
Apply the solve as you do for the linear case with Sherman Morisson formula
<br>
}<br>
<br>
This you can run with for example -yourpc_pc_type cholesky<br>
<br>
Barry<br>
<br>
Looks complicated, conceptually simple.<br>
<br>
<br>
> 2) Call KSPSetOperators(ksp, K, K,)<br>
> <br>
> 3) Solve the first system KSPSolve(ksp, -F, u_1)<br>
> <br>
> 4) Solve the second system KSPSolve(ksp, a, u_2)<br>
> <br>
> 5) Calculate beta VecDot(a, u_2, &gamma); beta = 1./(1. + k*gamma);<br>
> <br>
> 6) Update the guess VecDot(u_2, F, &delta); VecAXPBYPCZ(u, 1.0, beta*delta, 1.0, u_1, u_2)<br>
<br>
No<br>
<br>
> On Feb 6, 2020, at 9:02 AM, Olek Niewiarowski <<a href="mailto:aan2@princeton.edu" target="_blank">aan2@princeton.edu</a>> wrote:<br>
> <br>
> Hi Matt, <br>
> <br>
> What you suggested in your last email was exactly what I did on my very first attempt at the problem, and while it "worked," convergence was not satisfactory due to the Newton step being fixed in step 6. This is the reason I would like to use the linesearch
in SNES instead. Indeed in your manual you "recommend most PETSc users work directly with SNES, rather than using PETSc for the linear problem within a nonlinear solver." Ideally I'd like to create a SNES solver, pass in the functions to evaluate K, F, a,
and k, and set up the underlying KSP object as in my first message. Is this possible?<br>
> <br>
> Thanks,<br>
> <br>
> Alexander (Olek) Niewiarowski<br>
> PhD Candidate, Civil & Environmental Engineering<br>
> Princeton University, 2020<br>
> Cell: +1 (610) 393-2978<br>
> From: Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br>
> Sent: Thursday, February 6, 2020 5:33<br>
> To: Olek Niewiarowski <<a href="mailto:aan2@princeton.edu" target="_blank">aan2@princeton.edu</a>><br>
> Cc: Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>>; <a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
> Subject: Re: [petsc-users] Implementing the Sherman Morisson formula (low rank update) in petsc4py and FEniCS?<br>
> <br>
> On Wed, Feb 5, 2020 at 8:53 PM Olek Niewiarowski <<a href="mailto:aan2@princeton.edu" target="_blank">aan2@princeton.edu</a>> wrote:<br>
> Hi Barry and Matt, <br>
> <br>
> Thank you for your input and for creating a new issue in the repo. <br>
> My initial question was more basic (how to configure the SNES's KSP solver as in my first message with a and k), but now I see there's more to the implementation. To reiterate, for my problem's structure, a good solution algorithm (on the algebraic level)
is the following "double back-substitution": <br>
> For each nonlinear iteration:<br>
> • define intermediate vectors u_1 and u_2 <br>
> • solve Ku_1 = -F --> u_1 = -K^{-1}F (this solve is cheap, don't actually need K^-1)<br>
> • solve Ku_2 = -a --> u_2 = -K^{-1}a (ditto)<br>
> • define \beta = 1/(1 + k a^Tu_2)<br>
> • \Delta u = u_1 + \beta k u_2^T F u_2<br>
> • u = u + \Delta u<br>
> This is very easy to setup: <br>
> <br>
> 1) Create a KSP object KSPCreate(comm, &ksp)<br>
> <br>
> 2) Call KSPSetOperators(ksp, K, K,)<br>
> <br>
> 3) Solve the first system KSPSolve(ksp, -F, u_1)<br>
> <br>
> 4) Solve the second system KSPSolve(ksp, a, u_2)<br>
> <br>
> 5) Calculate beta VecDot(a, u_2, &gamma); beta = 1./(1. + k*gamma);<br>
> <br>
> 6) Update the guess VecDot(u_2, F, &delta); VecAXPBYPCZ(u, 1.0, beta*delta, 1.0, u_1, u_2)<br>
> <br>
> Thanks,<br>
> <br>
> Matt<br>
> <br>
> I don't need the Jacobian inverse, [K−kaaT]-1 = K-1 - (kK-1 aaTK-1)/(1+kaTK-1a) just the solution Δu = [K−kaaT]-1F = K-1F - (kK-1 aFK-1a)/(1 + kaTK-1a)
<br>
> = u_1 + beta k u_2^T F u_2 (so I never need to invert K either). (To Matt's point on gitlab, K is a symmetric sparse matrix arising from a bilinear form. ) Also, eventually, I want to have more than one low-rank updates to K, but again, Sherman Morrisson
Woodbury should still work. <br>
> <br>
> Being new to PETSc, I don't know if this algorithm directly translates into an efficient numerical solver. (I'm also not sure if Picard iteration would be useful here.) What would it take to set up the KSP solver in SNES like I did below? Is it possible "out
of the box"? I looked at MatCreateLRC() - what would I pass this to? (A pointer to demo/tutorial would be very appreciated.) If there's a better way to go about all of this, I'm open to any and all ideas. My only limitation is that I use petsc4py exclusively
since I/future users of my code will not be comfortable with C.<br>
> <br>
> Thanks again for your help!<br>
> <br>
> <br>
> Alexander (Olek) Niewiarowski<br>
> PhD Candidate, Civil & Environmental Engineering<br>
> Princeton University, 2020<br>
> Cell: +1 (610) 393-2978<br>
> From: Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>><br>
> Sent: Wednesday, February 5, 2020 15:46<br>
> To: Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br>
> Cc: Olek Niewiarowski <<a href="mailto:aan2@princeton.edu" target="_blank">aan2@princeton.edu</a>>; <a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
> Subject: Re: [petsc-users] Implementing the Sherman Morisson formula (low rank update) in petsc4py and FEniCS?<br>
> <br>
> <br>
> <a href="https://gitlab.com/petsc/petsc/issues/557" target="_blank">https://gitlab.com/petsc/petsc/issues/557</a><br>
> <br>
> <br>
> > On Feb 5, 2020, at 7:35 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br>
> > <br>
> > Perhaps Barry is right that you want Picard, but suppose you really want Newton.<br>
> > <br>
> > "This problem can be solved efficiently using the Sherman-Morrison formula" Well, maybe. The main assumption here is that inverting K is cheap. I see two things you can do in a straightforward way:<br>
> > <br>
> > 1) Use MatCreateLRC() <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateLRC.html" target="_blank">
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateLRC.html</a> to create the Jacobian<br>
> > and solve using an iterative method. If you pass just K was the preconditioning matrix, you can use common PCs.<br>
> > <br>
> > 2) We only implemented MatMult() for LRC, but you could stick your SMW code in for MatSolve_LRC if you really want to factor K. We would<br>
> > of course help you do this.<br>
> > <br>
> > Thanks,<br>
> > <br>
> > Matt<br>
> > <br>
> > On Wed, Feb 5, 2020 at 1:36 AM Smith, Barry F. via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br>
> > <br>
> > I am not sure of everything in your email but it sounds like you want to use a "Picard" iteration to solve [K(u)−kaaT]Δu=−F(u). That is solve<br>
> > <br>
> > A(u^{n}) (u^{n+1} - u^{n}) = F(u^{n}) - A(u^{n})u^{n} where A(u) = K(u) - kaaT<br>
> > <br>
> > PETSc provides code to this with SNESSetPicard() (see the manual pages) I don't know if Petsc4py has bindings for this.<br>
> > <br>
> > Adding missing python bindings is not terribly difficult and you may be able to do it yourself if this is the approach you want.<br>
> > <br>
> > Barry<br>
> > <br>
> > <br>
> > <br>
> > > On Feb 4, 2020, at 5:07 PM, Olek Niewiarowski <<a href="mailto:aan2@princeton.edu" target="_blank">aan2@princeton.edu</a>> wrote:<br>
> > > <br>
> > > Hello, <br>
> > > 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.
<a href="http://fourier.eng.hmc.edu/e176/lectures/algebra/node6.html" target="_blank">http://fourier.eng.hmc.edu/e176/lectures/algebra/node6.html</a> ). 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 :
<br>
> > > <br>
> > > [K−kaaT]-1 = K-1 - (kK-1 aaTK-1)/(1+kaTK-1a)<br>
> > > 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
<a href="https://www.firedrakeproject.org/petsc-interface.html#defining-a-preconditioner" target="_blank">
https://www.firedrakeproject.org/petsc-interface.html#defining-a-preconditioner</a> :<br>
> > > • while (norm(delU) > alpha): # while not converged<br>
> > > • <br>
> > > • self.update_F() # call to method to update r.h.s form<br>
> > > • self.update_K() # call to update the jacobian form<br>
> > > • K = assemble(self.K) # assemble the jacobian matrix<br>
> > > • F = assemble(self.F) # assemble the r.h.s vector<br>
> > > • a = assemble(self.a_form) # assemble the a_form (see Sherman Morrison formula)<br>
> > > • <br>
> > > • for bc in self.mem.bc: # apply boundary conditions<br>
> > > • bc.apply(K, F) <br>
> > > • bc.apply(K, a) <br>
> > > • <br>
> > > • B = PETSc.Mat().create() <br>
> > > • <br>
> > > • # Assemble the bilinear form that defines A and get the concrete
<br>
> > > • # PETSc matrix <br>
> > > • A = as_backend_type(K).mat() # get the PETSc objects for K and a<br>
> > > • u = as_backend_type(a).vec() <br>
> > > • <br>
> > > • # Build the matrix "context" # see firedrake docs<br>
> > > • Bctx = MatrixFreeB(A, u, u, self.k) <br>
> > > • <br>
> > > • # Set up B <br>
> > > • # B is the same size as A <br>
> > > • B.setSizes(*A.getSizes()) <br>
> > > • B.setType(B.Type.PYTHON) <br>
> > > • B.setPythonContext(Bctx) <br>
> > > • B.setUp() <br>
> > > • <br>
> > > • <br>
> > > • ksp = PETSc.KSP().create() # create the KSP linear solver object<br>
> > > • ksp.setOperators(B) <br>
> > > • ksp.setUp() <br>
> > > • pc = ksp.pc <br>
> > > • pc.setType(pc.Type.PYTHON) <br>
> > > • pc.setPythonContext(MatrixFreePC()) <br>
> > > • ksp.setFromOptions() <br>
> > > • <br>
> > > • solution = delU # the incremental displacement at this iteration<br>
> > > • <br>
> > > • b = as_backend_type(-F).vec() <br>
> > > • delu = solution.vector().vec() <br>
> > > • <br>
> > > • ksp.solve(b, delu) <br>
> > > <br>
> > > • self.mem.u.vector().axpy(0.25, self.delU.vector()) # poor man's linesearch<br>
> > > • counter += 1 <br>
> > > Here is the corresponding petsc4py code adapted from the firedrake docs:<br>
> > > <br>
> > > • class MatrixFreePC(object): <br>
> > > • <br>
> > > • def setUp(self, pc): <br>
> > > • B, P = pc.getOperators() <br>
> > > • # extract the MatrixFreeB object from B <br>
> > > • ctx = B.getPythonContext() <br>
> > > • self.A = ctx.A <br>
> > > • self.u = ctx.u <br>
> > > • self.v = ctx.v <br>
> > > • self.k = ctx.k <br>
> > > • # Here we build the PC object that uses the concrete, <br>
> > > • # assembled matrix A. We will use this to apply the action
<br>
> > > • # of A^{-1} <br>
> > > • self.pc = PETSc.PC().create() <br>
> > > • self.pc.setOptionsPrefix("mf_") <br>
> > > • self.pc.setOperators(self.A) <br>
> > > • self.pc.setFromOptions() <br>
> > > • # Since u and v do not change, we can build the denominator
<br>
> > > • # and the action of A^{-1} on u only once, in the setup <br>
> > > • # phase. <br>
> > > • tmp = self.A.createVecLeft() <br>
> > > • self.pc.apply(self.u, tmp) <br>
> > > • self._Ainvu = tmp <br>
> > > • self._denom = 1 + self.k*self.v.dot(self._Ainvu) <br>
> > > • <br>
> > > • def apply(self, pc, x, y): <br>
> > > • # y <- A^{-1}x <br>
> > > • self.pc.apply(x, y) <br>
> > > • # alpha <- (v^T A^{-1} x) / (1 + v^T A^{-1} u) <br>
> > > • alpha = (self.k*self.v.dot(y)) / self._denom <br>
> > > • # y <- y - alpha * A^{-1}u <br>
> > > • y.axpy(-alpha, self._Ainvu) <br>
> > > • <br>
> > > • <br>
> > > • class MatrixFreeB(object): <br>
> > > • <br>
> > > • def __init__(self, A, u, v, k): <br>
> > > • self.A = A <br>
> > > • self.u = u <br>
> > > • self.v = v <br>
> > > • self.k = k <br>
> > > • <br>
> > > • def mult(self, mat, x, y): <br>
> > > • # y <- A x <br>
> > > • self.A.mult(x, y) <br>
> > > • <br>
> > > • # alpha <- v^T x <br>
> > > • alpha = self.v.dot(x) <br>
> > > • <br>
> > > • # y <- y + alpha*u <br>
> > > • y.axpy(alpha, self.u) <br>
> > > 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:<br>
> > > solver = PETScSNESSolver() # the FEniCS SNES wrapper<br>
> > > snes = solver.snes() # the petsc4py SNES object<br>
> > > ## ??<br>
> > > ksp = snes.getKSP()<br>
> > > # set ksp option similar to above<br>
> > > solver.solve()<br>
> > > <br>
> > > 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!).
<br>
> > > Many thanks in advance!<br>
> > > Alex<br>
> > <br>
> > <br>
> > <br>
> > -- <br>
> > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> > -- Norbert Wiener<br>
> > <br>
> > <a href="https://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
> <br>
> <br>
> <br>
> -- <br>
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> -- Norbert Wiener<br>
> <br>
> <a href="https://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
<br>
</div>
</span></font></div>
</div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>