<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<style type="text/css" style="display:none;"> P {margin-top:0;margin-bottom:0;} </style>
</head>
<body dir="ltr">
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
</div>
<div style="margin: 0px; font-size: 15px; font-family: Helvetica, Arial, sans-serif, serif, EmojiFont; color: rgb(34, 34, 34); background-color: rgb(255, 255, 255)">
<span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont">Hello, </span></div>
<div style="margin: 0px; font-size: 15px; font-family: Helvetica, Arial, sans-serif, serif, EmojiFont; color: rgb(34, 34, 34); background-color: rgb(255, 255, 255)">
<span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont">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
<span style="font-family: Helvetica, Arial, sans-serif, serif, EmojiFont; background-color: rgb(255, 255, 255); display: inline !important">
(e.g. </span><a href="http://fourier.eng.hmc.edu/e176/lectures/algebra/node6.html" target="_blank" rel="noopener noreferrer" style="margin: 0px; font-family: Helvetica, Arial, sans-serif, serif, EmojiFont; background-color: rgb(255, 255, 255)">http://fourier.eng.hmc.edu/e176/lectures/algebra/node6.html</a><span style="font-family: Helvetica, Arial, sans-serif, serif, EmojiFont; background-color: rgb(255, 255, 255); display: inline !important"> )</span>.</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"> I
am solving a nonlinear system of the form</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"> </span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">[</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal"><b>K(u)</b></span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">−</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">k<b>aa</b></span><sup style="font-family: MJXc-TeX-main-B, MJXc-TeX-main-Bx, MJXc-TeX-main-Bw; word-spacing: normal"><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"><i>T</i></span></sup><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">]</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">Δ</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal"><b>u</b></span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">=</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">−</span><span style="margin: 0px; word-spacing: normal"><span style="margin: 0px; font-size: 11pt"><font face="MJXc-TeX-main-B, MJXc-TeX-main-Bx, MJXc-TeX-main-Bw"><span style="margin: 0px; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"><b>F(u)</b></span></font></span></span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">. </span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">Here
the jacobian matrix</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal"> </span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">K</span><span style="margin: 0px"><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"> </span></span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont">is
modified by the term </span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">k<b>aa</b></span><sup style="font-family: MJXc-TeX-main-B, MJXc-TeX-main-Bx, MJXc-TeX-main-Bw; word-spacing: normal"><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"><i>T</i></span></sup><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">,
w</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">here k is a scalar. </span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal"> Notably,</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal"> </span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">K</span><span style="margin: 0px"><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"> </span></span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont">is
sparse, while the term</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"> k<b>aa</b></span><sup style="font-family: MJXc-TeX-main-B, MJXc-TeX-main-Bx, MJXc-TeX-main-Bw"><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"><i>T</i></span></sup><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal"> r</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont">esults
in a full matrix. </span><span style="margin: 0px; font-size: 11pt">This problem can be solved efficiently using the Sherman-Morrison formula : </span></div>
<div style="margin: 0px; font-size: 15px; font-family: Helvetica, Arial, sans-serif, serif, EmojiFont; color: rgb(34, 34, 34); background-color: rgb(255, 255, 255)">
<span style="margin: 0px; font-size: 11pt"><br>
</span></div>
<div style="margin: 0px; font-size: 15px; font-family: Helvetica, Arial, sans-serif, serif, EmojiFont; color: rgb(34, 34, 34); background-color: rgb(255, 255, 255)">
<span style="margin: 0px; font-size: 14.6667px"><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">[</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal"><b>K</b></span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">−</span><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">k<b>aa</b></span><sup style="word-spacing: normal; font-family: MJXc-TeX-main-B, MJXc-TeX-main-Bx, MJXc-TeX-main-Bw"><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"><i>T</i></span></sup><span style="margin: 0px; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; word-spacing: normal">]<sup>-1<span> </span></sup>=<span> </span><b>K</b><sup>-1 </sup> -
(k<span style="margin: 0px; background-color: rgb(255, 255, 255); display: inline !important"><b>K</b></span><sup style="font-family: Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255)">-1</sup><span style="margin: 0px; background-color: rgb(255, 255, 255)"> <b>a</b><span style="margin: 0px; background-color: rgb(255, 255, 255)"><b>a</b></span><sup style="background-color: rgb(255, 255, 255); font-family: MJXc-TeX-main-B, MJXc-TeX-main-Bx, MJXc-TeX-main-Bw"><span style="margin: 0px; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"><i>T</i></span></sup></span><span style="margin: 0px; background-color: rgb(255, 255, 255); display: inline !important"><b>K</b></span><sup style="font-family: Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255)">-1</sup><span style="margin: 0px; background-color: rgb(255, 255, 255)">)</span>/(1+k<span style="margin: 0px; background-color: rgb(255, 255, 255)"><span style="margin: 0px; background-color: rgb(255, 255, 255)"><b>a</b></span><sup style="background-color: rgb(255, 255, 255); font-family: MJXc-TeX-main-B, MJXc-TeX-main-Bx, MJXc-TeX-main-Bw"><span style="margin: 0px; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont"><i>T</i></span></sup></span><span style="margin: 0px; background-color: rgb(255, 255, 255); display: inline !important"><b>K</b></span><sup style="font-family: Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255)">-1</sup><span style="margin: 0px; background-color: rgb(255, 255, 255)"><b>a</b>)</span></span></span></div>
<div style="margin: 0px; font-size: 15px; font-family: "Segoe UI", "Segoe UI Web (West European)", "Segoe UI", -apple-system, system-ui, Roboto, "Helvetica Neue", sans-serif; color: rgb(32, 31, 30); background-color: rgb(255, 255, 255)">
<span tabindex="0" style="margin: 0px; padding: 1px 0px; line-height: 0; outline: 0px; display: inline-block; text-align: left; max-width: none; max-height: none; min-width: 0px; min-height: 0px"><span style="margin: 0px; display: inline-block; border-collapse: separate; border-spacing: 0px"><span style="margin: 0px; display: inline-block; box-sizing: content-box !important"><span style="margin: 0px; display: inline-block; box-sizing: content-box !important"><span style="margin: 0px; padding: 0.455em 0px 0.568em; font-weight: normal; font-size: 11pt; line-height: normal; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; color: rgb(34, 34, 34); background-color: rgb(255, 255, 255); word-spacing: normal; display: block; box-sizing: content-box !important">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" rel="noopener noreferrer" style="margin: 0px">https://www.firedrakeproject.org/petsc-interface.html#defining-a-preconditioner</a> :</span><span style="margin: 0px; padding: 0.455em 0px 0.568em; line-height: normal; display: block; box-sizing: content-box !important">
<li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"><span style="margin: 0px; border: 0px none"> </span><span style="margin: 0px; border: 0px none; font-weight: bold; color: rgb(0, 102, 153)">while</span><span style="margin: 0px; border: 0px none"> (norm(</span><span style="margin: 0px; border: 0px none">delU) > alpha):
# while not converged</span></span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.update_F() # call to method to update r.h.s form</span></span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.update_K() # call to update the jacobian form</span></span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> K = assemble(<span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.K) # assemble the jacobian matrix</span></span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> F = assemble(<span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.F) # assemble the r.h.s vector</span></span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> a = assemble(<span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.a_form) # assemble the a_form (see Sherman Morrison formula)</span></span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; font-weight: bold; color: rgb(0, 102, 153)">for</span><span style="margin: 0px; border: 0px none"> bc </span><span style="margin: 0px; border: 0px none; font-weight: bold; color: rgb(0, 102, 153)">in</span><span style="margin: 0px; border: 0px none"> </span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.mem.bc:
# apply boundary conditions</span></span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> bc.apply(K, F) </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> bc.apply(K, a) </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> B = PETSc.Mat().create() </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># Assemble the bilinear form that defines A and get the concrete</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># PETSc matrix</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> A = as_backend_type(K).mat() # get the PETSc objects for K and a<br>
</span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> u = as_backend_type(a).vec() </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># Build the matrix "context"</span><span style="margin: 0px; border: 0px none"> # see firedrake docs</span></span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> Bctx = MatrixFreeB(A, u, u, <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.k) </span></span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># Set up B</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># B is the same size as A</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> B.setSizes(*A.getSizes()) </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> B.setType(B.Type.PYTHON) </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> B.setPythonContext(Bctx) </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> B.setUp() </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> ksp = PETSc.KSP().create() # create the KSP linear solver object</span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> ksp.setOperators(B) </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> ksp.setUp() </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> pc = ksp.pc </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> pc.setType(pc.Type.PYTHON) </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> pc.setPythonContext(MatrixFreePC()) </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> ksp.setFromOptions() </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> solution = <span style="margin: 0px; border: 0px none">delU # the incremental displacement at this iteration</span></span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> b = as_backend_type(-F).vec() </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> delu = solution.vector().vec() </span></li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> ksp.solve(b, delu) <br>
</span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<br>
</li><li style="background-color: rgb(248, 248, 248); color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.mem.u.vector().axpy(</span><span style="margin: 0px; border: 0px none">0.25</span><span style="margin: 0px; border: 0px none">, </span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.delU.vector())
# poor man's linesearch</span></span></li><li style="background-color: white; color: rgb(92, 92, 92); font-family: Consolas, "Courier New", Courier, mono, serif; font-size: 12px; font-weight: normal; word-spacing: normal; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> counter += <span style="margin: 0px; border: 0px none">1</span><span style="margin: 0px; border: 0px none"> <br>
</span></span></li><span style="margin: 0px; text-align: start; display: inline !important"><span style="margin: 0px; padding: 0.455em 0px 0.568em; font-weight: normal; font-size: 11pt; font-family: Arial, Helvetica, sans-serif, serif, EmojiFont; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255); word-spacing: normal; display: block; box-sizing: content-box !important"><span style="margin: 0px; font-size: medium; font-family: Calibri, Arial, Helvetica, sans-serif, serif, EmojiFont; background-color: rgb(255, 255, 255); display: inline !important">
<div style="margin: 0px; font-size: 15px; color: black; background-color: rgb(255, 255, 255)">
<font size="3">Here is the corresponding petsc4py code adapted from the firedrake docs:</font></div>
<div style="margin: 0px; font-size: 15px; color: black; background-color: rgb(255, 255, 255)">
<font size="3"><br>
</font></div>
<div style="margin: 0px; font-size: 15px; color: black; background-color: rgb(255, 255, 255)">
<font size="3">
<ol start="1" style="color: rgb(92, 92, 92); font-size: 12px; font-family: Consolas, "Courier New", Courier, mono, serif; background-color: white; list-style: decimal; margin: 0px 0px 1px 45px; border-style: none">
<li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"><span style="margin: 0px; border: 0px none; font-weight: bold; color: rgb(0, 102, 153)">class</span><span style="margin: 0px; border: 0px none"> MatrixFreePC(object): </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; font-weight: bold; color: rgb(0, 102, 153)">def</span><span style="margin: 0px; border: 0px none"> setUp(</span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">, pc): </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> B, P = pc.getOperators() </span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># extract the MatrixFreeB object from B</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> ctx = B.getPythonContext() </span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.A = ctx.A </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.u = ctx.u </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.v = ctx.v </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.k = ctx.k </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># Here we build the PC object that uses the concrete,</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># assembled matrix A. We will use this to apply the action</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># of A^{-1}</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.pc = PETSc.PC().create() </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.pc.setOptionsPrefix(</span><span style="margin: 0px; border: 0px none; color: blue">"mf_"</span><span style="margin: 0px; border: 0px none">) </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.pc.setOperators(</span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.A) </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.pc.setFromOptions() </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># Since u and v do not change, we can build the denominator</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># and the action of A^{-1} on u only once, in the setup</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># phase.</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> tmp = <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.A.createVecLeft() </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.pc.apply(</span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.u, tmp) </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">._Ainvu = tmp </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">._denom = </span><span style="margin: 0px; border: 0px none">1</span><span style="margin: 0px; border: 0px none"> + </span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.k*</span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.v.dot(</span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">._Ainvu) </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; font-weight: bold; color: rgb(0, 102, 153)">def</span><span style="margin: 0px; border: 0px none"> apply(</span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">, pc, x, y): </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># y <- A^{-1}x</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.pc.apply(x, y) </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># alpha <- (v^T A^{-1} x) / (1 + v^T A^{-1} u)</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> alpha = (<span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.k*</span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.v.dot(y)) / </span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">._denom </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># y <- y - alpha * A^{-1}u</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> y.axpy(-alpha, <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">._Ainvu) </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"><span style="margin: 0px; border: 0px none; font-weight: bold; color: rgb(0, 102, 153)">class</span><span style="margin: 0px; border: 0px none"> MatrixFreeB(object): </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; font-weight: bold; color: rgb(0, 102, 153)">def</span><span style="margin: 0px; border: 0px none"> __init__(</span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">, A, u, v, k): </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.A = A </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.u = u </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.v = v </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.k = k </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; font-weight: bold; color: rgb(0, 102, 153)">def</span><span style="margin: 0px; border: 0px none"> mult(</span><span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">, mat, x, y): </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># y <- A x</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.A.mult(x, y) </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># alpha <- v^T x</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> alpha = <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.v.dot(x) </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> </span></li><li style="background-color: rgb(248, 248, 248); list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> <span style="margin: 0px; border: 0px none; color: rgb(0, 130, 0)"># y <- y + alpha*u</span><span style="margin: 0px; border: 0px none"> </span></span></li><li style="background-color: white; list-style: decimal-leading-zero; margin: 0px; padding: 0px 3px 0px 10px; border-style: none none none solid; border-left-width: 3px; border-left-color: rgb(108, 226, 108); line-height: 14px">
<span style="margin: 0px; border: 0px none; color: black"> y.axpy(alpha, <span style="margin: 0px; border: 0px none">self</span><span style="margin: 0px; border: 0px none">.u) </span></span></li></ol>
</font></div>
</span></span><span style="margin: 0px; padding: 0.455em 0px 0.568em; display: block; box-sizing: content-box !important">Ho<span style="margin: 0px">wever, 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<span> </span><b>F</b>,<span> </span><b>K</b>, and<span> </span><b>a</b>,
I'd like to do something along the lines of:<br>
</span>
<div style="margin: 0px">solver = PETScSNESSolver() # the FEniCS SNES wrapper<br>
</div>
<span style="margin: 0px">snes = solver.snes() # the petsc4py SNES object</span>
<div style="margin: 0px">## ??<br>
<div style="margin: 0px">ksp = snes.getKSP()</div>
<div style="margin: 0px"> # set ksp option similar to above</div>
<div style="margin: 0px">solver.solve()</div>
<div style="margin: 0px">
<div style="margin: 0px; font-size: 12pt; font-family: Calibri, Arial, Helvetica, sans-serif, serif, EmojiFont; color: rgb(0, 0, 0)">
<br>
</div>
<div style="margin: 0px; font-size: 12pt; font-family: Calibri, Arial, Helvetica, sans-serif, serif, EmojiFont; color: rgb(0, 0, 0)">
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!). </div>
</div>
</div>
</span><span style="margin: 0px; padding: 0.455em 0px 0.568em; font-size: 12pt; font-family: Calibri, Arial, Helvetica, sans-serif, serif, EmojiFont; color: rgb(0, 0, 0); display: block; box-sizing: content-box !important">Many thanks in advance!</span></span></span></span></span></span></span></div>
<div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Alex</div>
</div>
</body>
</html>