[petsc-users] Different precision from MatAssembly/MatView
Matthew Knepley
knepley at gmail.com
Sun Mar 15 17:04:44 CDT 2026
On Sun, Mar 15, 2026 at 5:24 PM Noam T. <dontbugthedevs at proton.me> wrote:
> 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.
>
>
> This is indeed the case.
>
>
> 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.
>
>
> I imagined this was simply a printing artifact. However, something is
> bugging me: the message from finite differences jacobian.
>
> 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.
>
> So I did a couple more checks, avoiding calls to "XView" type of
> functions, and getting a single value directly:
>
> (For comparison) Computed matrix in the jacobian function (used to
> assemble)
> ----
> 0.531360000000 0.066670000000 0.197333333327
> 0.066670000000 0.535360000000 -0.003333333327
> 0.197333333327 -0.003333333327 0.527413333333
> ----
>
> ** Calling "MatGetValues" and then printing them all:
>
> 5.313866666667E-01 6.667333333333E-02 1.973381999940E-01
> 6.667333333333E-02 5.353866666667E-01 -3.338499994017E-03
> 1.973381999940E-01 -3.338499994017E-03 5.274399026665E-01
>
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
1) Run in the debugger
2) Stop in MatSetValues()
3) Print out the input array v[] (you should get the first matrix in full
precision)
4) You could tace down to where each value is entered, but it might be
easier to
5) Call https://urldefense.us/v3/__https://petsc.org/main/manualpages/Mat/MatSeqAIJGetArray/__;!!G_uCfscf7eWS!Zfo59ELDNQEj7_52YmMFINuHw7PsOrHM9FLA6nEtPQoZb6T0J-cwsVCJhIILqNlPTSNx1bunO62Z1EToyltK$ and look
directly at the array right after the call
to MatSetValues()
Thanks,
Matt
> These are the same values shown as when calling MatView.
> They are all quite off. There is for example a difference of 5e-4 for
> entry (0,0).
> Just in case I messed up the precision, the closest single-precision value
> for entry (0,0) is 0x1.100e6p-1 = 0.53135997...
>
> ** 2-norm of the assembled matrix times a vector [1,1,1]
> 1.2294717692646715
>
> This is _exactly_ the norm using the matrix from MatView (i.e. the one
> from MatGetValues) times [1,1,1].
>
> The original/assembled matrix times [1,1,1] should be
> 1.229421704786441 -- difference of 5e-5
>
> Am I accessing stored/assembled values the wrong way so that they all are
> always off when printed? Or am I missing something else?
>
> Thanks
>
>
> On Sunday, March 15th, 2026 at 6:16 PM, Matthew Knepley <knepley at gmail.com>
> wrote:
>
> On Sun, Mar 15, 2026 at 12:58 PM Noam T. <dontbugthedevs at proton.me> wrote:
>
>> The results above of the SNES output were using the flags
>>
>> -ksp_type fgmres -pc_type lu
>>
>
> 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.
>
> 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.
>
> Thanks,
>
> Matt
>
>> Thanks
>> On Sunday, March 15th, 2026 at 5:12 PM, Matthew Knepley <
>> knepley at gmail.com> wrote:
>>
>> On Sun, Mar 15, 2026 at 11:19 AM Noam T. via petsc-users <
>> petsc-users at mcs.anl.gov> wrote:
>>
>>> A bit more information:
>>>
>>> Looking at the solution of the system of equations, knowing the exact
>>> RHS vector:
>>> ---
>>> [-0.00197, 0.00203, -0.00197]
>>> ---
>>>
>>> SNES gives the solution:
>>> ---
>>> 3.316262462871E-03, -4.189244774965E-03, 2.468317119413E-03
>>> ---
>>>
>>> which is indeed (closer to) the solution for the jacobian shown in
>>> MatView (the "lower precision" one). The "exact" solution would be:
>>> ---
>>> 3.282607248309093b-3, -4.241572005990172b-3,, 2.425030835360701b-3
>>> ---
>>>
>>> which is already different form the second decimal digit.
>>>
>>
>> Are you using an exact solution method, like LU?
>>
>> Thanks,
>>
>> Matt
>>
>>>
>>> ***
>>> 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.:
>>>
>>> Computed (printed from in-code):
>>> ---
>>> 0.536488494840 0.066474939660 0.198813816284
>>> 0.066474939660 0.529482312614 -0.002835445071
>>> 0.198813816284 -0.002835445071 0.530589055263
>>> ---
>>>
>>> MatView:
>>> ---
>>> row 0: (0, 0.538491) (1, 0.0662019) (2, 0.198808)
>>> row 1: (0, 0.0662019) (1, 0.526981) (2, -0.00257554)
>>> row 2: (0, 0.198808) (1, -0.00257554) (2, 0.529883)
>>> ---
>>>
>>> Only one digit is correct/the same in entries [1,2], [2,1] and [2,2].
>>>
>>> Thank you.
>>> Noam
>>>
>>>
>>> On Sunday, March 15th, 2026 at 3:29 PM, Noam T. <
>>> dontbugthedevs at proton.me> wrote:
>>>
>>> Hello,
>>>
>>> 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.
>>>
>>> In the jacobian function, after looping through all quadrature points,
>>> printing the matrix to stdout shows:
>>>
>>> ----
>>> 0.531360000000 0.066670000000 0.197333333327
>>> 0.066670000000 0.535360000000 -0.003333333327
>>> 0.197333333327 -0.003333333327 0.527413333333
>>> ----
>>>
>>> (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).
>>>
>>> This same matrix is the one added to the global jacobian during the
>>> quadrature loop, with "DMPlexMatSetClosure". After looping, calling
>>> "MatAssemblyBegin/End", then MatView:
>>>
>>> ----
>>> Mat Object: 1 MPI process
>>> type: seqaij
>>> row 0: (0, 0.531387) (1, 0.0666733) (2, 0.197338)
>>> row 1: (0, 0.0666733) (1, 0.535387) (2, -0.0033385)
>>> row 2: (0, 0.197338) (1, -0.0033385) (2, 0.52744)
>>> ----
>>>
>>> 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.
>>>
>>> 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):
>>>
>>> ----
>>> ||J - Jfd||_F/||J||_F = 4.90736e-05, ||J - Jfd||_F = 4.74265e-05
>>> Hand-coded Jacobian ----------
>>> Mat Object: 1 MPI process
>>> type: seqaij
>>> row 0: (0, 0.531387) (1, 0.0666733) (2, 0.197338)
>>> row 1: (0, 0.0666733) (1, 0.535387) (2, -0.0033385)
>>> row 2: (0, 0.197338) (1, -0.0033385) (2, 0.52744)
>>> Finite difference Jacobian ----------
>>> Mat Object: 1 MPI process
>>> type: seqaij
>>> row 0: (0, 0.53136) (1, 0.06667) (2, 0.197333)
>>> row 1: (0, 0.06667) (1, 0.53536) (2, -0.00333333)
>>> row 2: (0, 0.197333) (1, -0.00333333) (2, 0.527413)
>>> Hand-coded minus finite-difference Jacobian with tolerance 1e-05
>>> ----------
>>> Mat Object: 1 MPI process
>>> type: seqaij
>>> row 0: (0, 2.66554e-05)
>>> row 1: (1, 2.66554e-05)
>>> row 2: (2, 2.6558e-05)
>>> ----
>>>
>>> The main code is in Fortran, using double precision variables, in case
>>> it is relevant.
>>>
>>> Thank you.
>>>
>>> Noam
>>>
>>>
>>>
>>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Zfo59ELDNQEj7_52YmMFINuHw7PsOrHM9FLA6nEtPQoZb6T0J-cwsVCJhIILqNlPTSNx1bunO62Z1JIkGzX5$
>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Zfo59ELDNQEj7_52YmMFINuHw7PsOrHM9FLA6nEtPQoZb6T0J-cwsVCJhIILqNlPTSNx1bunO62Z1MKa9hqJ$ >
>>
>>
>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Zfo59ELDNQEj7_52YmMFINuHw7PsOrHM9FLA6nEtPQoZb6T0J-cwsVCJhIILqNlPTSNx1bunO62Z1JIkGzX5$
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Zfo59ELDNQEj7_52YmMFINuHw7PsOrHM9FLA6nEtPQoZb6T0J-cwsVCJhIILqNlPTSNx1bunO62Z1MKa9hqJ$ >
>
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Zfo59ELDNQEj7_52YmMFINuHw7PsOrHM9FLA6nEtPQoZb6T0J-cwsVCJhIILqNlPTSNx1bunO62Z1JIkGzX5$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Zfo59ELDNQEj7_52YmMFINuHw7PsOrHM9FLA6nEtPQoZb6T0J-cwsVCJhIILqNlPTSNx1bunO62Z1MKa9hqJ$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260315/8d8b6e4f/attachment-0001.html>
More information about the petsc-users
mailing list