<div style="font-family: Arial, sans-serif; font-size: 14px; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">That was the claim indeed. And it held for a while. Eventually figured out the reason for the differences and, not unexpectedly, the error was on my side. There was an additional operation on the small matrix that I consistently missed every single time I tested.</div><div style="font-family: Arial, sans-serif; font-size: 14px; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);"><br></div><div style="font-family: Arial, sans-serif; font-size: 14px; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">Values set with MatSetClosure / MatSetValues / etc are the expected ones, and from MatView they are also correct, just rounded for visualization/printing.</div><div style="font-family: Arial, sans-serif; font-size: 14px; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);"><br></div><div style="font-family: Arial, sans-serif; font-size: 14px; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">Thanks.</div><div style="font-family: Arial, sans-serif; font-size: 14px; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);"><br></div><div class="protonmail_quote">
On Sunday, March 15th, 2026 at 11:04 PM, Matthew Knepley <knepley@gmail.com> wrote:<br>
<blockquote class="protonmail_quote" type="cite">
<div dir="ltr"><div dir="ltr">On Sun, Mar 15, 2026 at 5:24 PM Noam T. <<a href="mailto:dontbugthedevs@proton.me" rel="noreferrer nofollow noopener">dontbugthedevs@proton.me</a>> wrote:</div><div class="gmail_quote gmail_quote_container"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px"></div><blockquote type="cite"><div dir="ltr"><div>This
is very easy to check. Multiply the solution by your matrix and
subtract the RHS. It must satisfy the tolerance you have for the solver.
If not, there is some other problem, but I imagine it does, which is
what we mean by the solution.</div></div></blockquote><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">This is indeed the case.</div><blockquote type="cite"><div><br></div><div dir="ltr"><div><div>To answer the
first question, there is no reduction in precision in MatSetValues(). I
am sure if you carry out the addition by hand you will see that the
result is within 1ulp as it should be from the standard.</div></div></div></blockquote><div style="font-family:Arial,sans-serif;font-size:14px"><div></div>
</div>
<div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">I imagined this was simply a printing artifact. However, something is bugging me: the message from finite differences jacobian.</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">The printed jacobian from finite differences is exactly the computed one (see below). If the differences from MatView were simply from printing, the error between finite difference and "in-code" should be (much closer to) zero, not 1e-5.</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">So I did a couple more checks, avoiding calls to "XView" type of functions, and getting a single value directly:</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">(For comparison) Computed matrix in the jacobian function (used to assemble)</div><div style="font-family:Arial,sans-serif;font-size:14px">----</div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:Menlo,Consolas,"Courier New",monospace">0.531360000000 0.066670000000 0.197333333327</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:Menlo,Consolas,"Courier New",monospace">0.066670000000 0.535360000000 -0.003333333327</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:Menlo,Consolas,"Courier New",monospace">0.197333333327 -0.003333333327 0.527413333333</span></div><div style="font-family:Arial,sans-serif;font-size:14px">----</div><div><br></div><div style="font-family:Arial,sans-serif;font-size:14px">** Calling "MatGetValues" and then printing them all:</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><span> 5.313866666667E-01 6.667333333333E-02 1.973381999940E-01</span><div><span> 6.667333333333E-02 5.353866666667E-01 -3.338499994017E-03</span></div><div><span> 1.973381999940E-01 -3.338499994017E-03 5.274399026665E-01</span></div></div></blockquote><div><br>So the claim is that you have some small matrix (the first) and you call MatSetValues() to put it in a big matrix, and you get different values in the big matrix? I am confident that this is impossible. So, in order to verify this, it should be easy to</div><div><br></div><div>1) Run in the debugger</div><div><br></div><div>2) Stop in MatSetValues()</div><div><br></div><div>3) Print out the input array v[] (you should get the first matrix in full precision)</div><div><br></div><div>4) You could tace down to where each value is entered, but it might be easier to</div><div><br></div><div>5) Call <a href="https://urldefense.us/v3/__https://petsc.org/main/manualpages/Mat/MatSeqAIJGetArray/__;!!G_uCfscf7eWS!aXSiC0ZBQh6RhRlIcJvFp6VmHB_7WBNFmlvzKWJ3T17gc57mGjdy7_EzYhgQEBKGnDaA5xeOXb-jQJXkF4ojUkqPVbrrFw--$" target="_blank" rel="noreferrer nofollow noopener">https://petsc.org/main/manualpages/Mat/MatSeqAIJGetArray/</a> and look directly at the array right after the call</div><div> to MatSetValues()</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 style="font-family:Arial,sans-serif;font-size:14px">These are the same values shown as when calling MatView.</div><div style="font-family:Arial,sans-serif;font-size:14px">They are all quite off. There is for example a difference of 5e-4 for entry (0,0).</div><div style="font-family:Arial,sans-serif;font-size:14px">Just in case I messed up the precision, the closest single-precision value for entry (0,0) is <span><span>0x1.100e6p-1 = </span>0.53135997...</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">** 2-norm of the assembled matrix times a vector [1,1,1]</div><div style="font-family:Arial,sans-serif;font-size:14px"><span>1.2294717692646715</span><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">This is _exactly_ the norm using the matrix from MatView (i.e. the one from MatGetValues) times [1,1,1].</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">The original/assembled matrix times [1,1,1] should be</div><div style="font-family:Arial,sans-serif;font-size:14px"><span>1.229421704786441 -- difference of 5e-5</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">Am I accessing stored/assembled values the wrong way so that they all are always off when printed? Or am I missing something else?</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">Thanks</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><br><div>
On Sunday, March 15th, 2026 at 6:16 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank" rel="noreferrer nofollow noopener">knepley@gmail.com</a>> wrote:<br>
<blockquote type="cite">
<div dir="ltr"><div dir="ltr">On Sun, Mar 15, 2026 at 12:58 PM Noam T. <<a href="mailto:dontbugthedevs@proton.me" rel="noreferrer nofollow noopener" target="_blank">dontbugthedevs@proton.me</a>> wrote:</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 style="font-family:Arial,sans-serif;font-size:14px">The results above of the SNES output were using the flags</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:Menlo,Consolas,"Courier New",monospace">-ksp_type fgmres -pc_type lu</span></div></blockquote><div><br></div><div>This is very easy to check. Multiply the solution by your matrix and subtract the RHS. It must satisfy the tolerance you have for the solver. If not, there is some other problem, but I imagine it does, which is what we mean by the solution.</div><div><br></div><div>To answer the first question, there is no reduction in precision in MatSetValues(). I am sure if you carry out the addition by hand you will see that the result is within 1ulp as it should be from the standard.</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 style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:Menlo,Consolas,"Courier New",monospace">Thanks</span></div><div>
On Sunday, March 15th, 2026 at 5:12 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com" rel="noreferrer nofollow noopener" target="_blank">knepley@gmail.com</a>> wrote:<br>
<blockquote type="cite">
<div dir="ltr"><div dir="ltr">On Sun, Mar 15, 2026 at 11:19 AM Noam T. via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" rel="noreferrer nofollow noopener" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:</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>A bit more information:</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Looking at the solution of the system of equations, knowing the exact RHS vector:</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">---</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace">[-0.00197, 0.00203, -0.00197]</span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">---</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">SNES gives the solution:</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">---</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace"> 3.316262462871E-03, -4.189244774965E-03, 2.468317119413E-03</span><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">---</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">which is indeed (closer to) the solution for the jacobian shown in MatView (the "lower precision" one). The "exact" solution would be:</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">---<br><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">3.282607248309093b-3, -4.241572005990172b-3,, 2.425030835360701b-3</span></div>---</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">which is already different form the second decimal digit.</div></blockquote><div><br></div><div>Are you using an exact solution method, like LU?</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 style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">***</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Requesting a low precision to have SNES do a few iterations, shows that the built jacobian seems to deviate more and more form the computed one after some iterations, e.g.:</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Computed (printed from in-code):</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">---</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace"> 0.536488494840 0.066474939660 0.198813816284</span><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> 0.066474939660 0.529482312614 -0.002835445071</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> 0.198813816284 -0.002835445071 0.530589055263</span></div><span></span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">---</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">MatView:</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">---</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace"><span>row 0: </span>(0, 0.538491) (1, 0.0662019) (2, 0.198808) </span><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"><span>row 1: </span>(0, 0.0662019) (1, 0.526981) (2, -0.00257554) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"><span>row 2: </span>(0, 0.198808) (1, -0.00257554) (2, 0.529883) </span></div><span></span>---</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Only one digit is correct/the same in entries [1,2], [2,1] and [2,2].</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Thank you.</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Noam</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div>
On Sunday, March 15th, 2026 at 3:29 PM, Noam T. <<a href="mailto:dontbugthedevs@proton.me" rel="noreferrer nofollow noopener" target="_blank">dontbugthedevs@proton.me</a>> wrote:<br>
<blockquote type="cite">
<div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">Hello,</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">Looking at the assembly of the Jacobian using MatView, I noticed values were somewhat different from the computed ones in the jacobian function (DMSNESSetJacobian). For the test, I used a single Q1 element, with 4 quadrature points, so the exact Jacobian matrix can be computed analytically. </span></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">In the jacobian function, after looping through all quadrature points, printing the matrix to stdout shows:</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">----</span><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:Menlo,Consolas,"Courier New",monospace"> 0.531360000000 0.066670000000 0.197333333327</span><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> 0.066670000000 0.535360000000 -0.003333333327</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> 0.197333333327 -0.003333333327 0.527413333333</span></div><span></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">----</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">(showing only the 3 free DoFs). The 2x2 submatrix in [0:1] has exactly four/five nonzero digits, as shown above. The rest of the elements have a periodic 3 ending (a couple of digits at the end are off, but that's fine).</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">This same matrix is the one added to the global jacobian during the quadrature loop, with "</span><span style="font-family:"Times New Roman",serif">DMPlexMatSetClosure". After looping, calling "MatAssemblyBegin/End", then MatView: </span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><span style="font-family:"Times New Roman",serif">----</span><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><span style="font-family:Menlo,Consolas,"Courier New",monospace">Mat Object: 1 MPI process</span><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> type: seqaij</span></div></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 0: (0, 0.531387) (1, 0.0666733) (2, 0.197338) </span><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 1: (0, 0.0666733) (1, 0.535387) (2, -0.0033385) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 2: (0, 0.197338) (1, -0.0033385) (2, 0.52744) </span></div><span></span><span style="font-family:"Times New Roman",serif">----</span><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">Values are close, but definitely not the same as computed. Are these values from "MatView" the ones actually stored in the matrix? It seems there is some precision loss in the process.</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">Interestingly, computing the jacobian through finite differences, shows the same result for the "hand-coded" jacobian, whereas the finite differences ones is the "exact" one (same as first one above):</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">----</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><span style="font-family:Menlo,Consolas,"Courier New",monospace"> ||J - Jfd||_F/||J||_F = 4.90736e-05, ||J - Jfd||_F = 4.74265e-05</span><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> Hand-coded Jacobian ----------</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">Mat Object: 1 MPI process</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> type: seqaij</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 0: (0, 0.531387) (1, 0.0666733) (2, 0.197338) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 1: (0, 0.0666733) (1, 0.535387) (2, -0.0033385) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 2: (0, 0.197338) (1, -0.0033385) (2, 0.52744) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> Finite difference Jacobian ----------</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">Mat Object: 1 MPI process</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> type: seqaij</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 0: (0, 0.53136) (1, 0.06667) (2, 0.197333) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 1: (0, 0.06667) (1, 0.53536) (2, -0.00333333) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 2: (0, 0.197333) (1, -0.00333333) (2, 0.527413) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> Hand-coded minus finite-difference Jacobian with tolerance 1e-05 ----------</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">Mat Object: 1 MPI process</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"> type: seqaij</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 0: (0, 2.66554e-05) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 1: (1, 2.66554e-05) </span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">row 2: (2, 2.6558e-05) </span></div><span></span></span><span style="font-family:"Times New Roman",serif"> ----</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">The main code is in Fortran, using double precision variables, in case it is relevant.</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">Thank you.</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px"><span style="font-family:"Times New Roman",serif">Noam</span></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div>
</blockquote><br>
</div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><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="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!aXSiC0ZBQh6RhRlIcJvFp6VmHB_7WBNFmlvzKWJ3T17gc57mGjdy7_EzYhgQEBKGnDaA5xeOXb-jQJXkF4ojUkqPVRtdch1Z$" rel="noreferrer nofollow noopener" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote><br>
</div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><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="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!aXSiC0ZBQh6RhRlIcJvFp6VmHB_7WBNFmlvzKWJ3T17gc57mGjdy7_EzYhgQEBKGnDaA5xeOXb-jQJXkF4ojUkqPVRtdch1Z$" rel="noreferrer nofollow noopener" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote><br>
</div></div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><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="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!aXSiC0ZBQh6RhRlIcJvFp6VmHB_7WBNFmlvzKWJ3T17gc57mGjdy7_EzYhgQEBKGnDaA5xeOXb-jQJXkF4ojUkqPVRtdch1Z$" target="_blank" rel="noreferrer nofollow noopener">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote><br>
</div>