[petsc-users] Different precision from MatAssembly/MatView
Matthew Knepley
knepley at gmail.com
Tue Mar 17 11:54:24 CDT 2026
On Tue, Mar 17, 2026 at 12:50 PM Noam T. <dontbugthedevs at proton.me> wrote:
> 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.
>
> Values set with MatSetClosure / MatSetValues / etc are the expected ones,
> and from MatView they are also correct, just rounded for
> visualization/printing.
>
Cool. Do you think there is diagnostic output I could give that would make
looking into the process easier? If you are using the Plex assembly, I have
-dm_plex_print_fem 5
to print out all the element vec/mats and local assembly. There is also
-petscds_print_integrate 5
for lower level integration output. Let me know and I can add something.
Thanks,
Matt
> Thanks.
>
> On Sunday, March 15th, 2026 at 11:04 PM, Matthew Knepley <
> knepley at gmail.com> wrote:
>
> 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!ecH8MVZGPPyWCLBwVWOE1kIET70KehuSuvpgwWWHXgtCx2q0ct8VJ5ErEr5YKL78tvU0Fs5KIoC-_H-0ue1e$ 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!ecH8MVZGPPyWCLBwVWOE1kIET70KehuSuvpgwWWHXgtCx2q0ct8VJ5ErEr5YKL78tvU0Fs5KIoC-_PFIFdfe$
>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!ecH8MVZGPPyWCLBwVWOE1kIET70KehuSuvpgwWWHXgtCx2q0ct8VJ5ErEr5YKL78tvU0Fs5KIoC-_Euofjy6$ >
>>>
>>>
>>>
>>
>> --
>> 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!ecH8MVZGPPyWCLBwVWOE1kIET70KehuSuvpgwWWHXgtCx2q0ct8VJ5ErEr5YKL78tvU0Fs5KIoC-_PFIFdfe$
>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!ecH8MVZGPPyWCLBwVWOE1kIET70KehuSuvpgwWWHXgtCx2q0ct8VJ5ErEr5YKL78tvU0Fs5KIoC-_Euofjy6$ >
>>
>>
>>
>
> --
> 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!ecH8MVZGPPyWCLBwVWOE1kIET70KehuSuvpgwWWHXgtCx2q0ct8VJ5ErEr5YKL78tvU0Fs5KIoC-_PFIFdfe$
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!ecH8MVZGPPyWCLBwVWOE1kIET70KehuSuvpgwWWHXgtCx2q0ct8VJ5ErEr5YKL78tvU0Fs5KIoC-_Euofjy6$ >
>
>
>
--
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!ecH8MVZGPPyWCLBwVWOE1kIET70KehuSuvpgwWWHXgtCx2q0ct8VJ5ErEr5YKL78tvU0Fs5KIoC-_PFIFdfe$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!ecH8MVZGPPyWCLBwVWOE1kIET70KehuSuvpgwWWHXgtCx2q0ct8VJ5ErEr5YKL78tvU0Fs5KIoC-_Euofjy6$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260317/8ff6463d/attachment-0001.html>
More information about the petsc-users
mailing list