[petsc-users] how is snes_test_jacobian_view computing the matrix difference?
Matteo Semplice
matteo.semplice at uninsubria.it
Tue Mar 3 02:44:42 CST 2026
Dear Barry,
indeed your suggestion fixes things! (I assume that you
meant PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));)
Now I get a consistent output.
Hand-coded Jac:
row 84: (84, 60.) (85, 0.) (3446, -80.) (3447,
0.) (6808, 20.) (6809, 0.)
row 85: (84, 0.00914434) (85, 1.) (3446, 0.)
(3447, 0.) (6808, 0.) (6809, 0.)
Finite-difference Jac:
row 84: (84, 60.) (85, 0.) (3446, -80.) (3447,
0.) (6808, 20.) (6809, 0.)
row 85: (84, 0.) (85, 0.999971) (3446, 0.) (3447,
0.) (6808, 0.) (6809, 0.)
Difference of the two matrices:
row 84: (84, 0.) (85, 0.) (3446, 0.) (3447, 0.)
(6808, 0.) (6809, 0.)
row 85: (84, 0.00914434) (85, 2.9293e-05) (3446, 0.)
(3447, 0.) (6808, 0.) (6809, 0.)
I have applied your suggestion on the top of the release branch and I
had to fix also the creation of matrix C. Here below is the complete
diff, in case you want to apply it to the official repository:
diff --git a/src/snes/interface/snes.c b/src/snes/interface/snes.c
index f66c23a1417..eee3d070d4c 100644
--- a/src/snes/interface/snes.c
+++ b/src/snes/interface/snes.c
@@ -2768,8 +2768,6 @@PetscErrorCode SNESTestJacobian(SNES snes,
PetscReal *Jnorm, PetscReal *diffNorm
Vec x = snes->vec_sol, f;
PetscReal nrm, gnorm;
PetscReal threshold = 1.e-5;
- MatType mattype;
- PetscInt m, n, M, N;
void *functx;
PetscBool complete_print = PETSC_FALSE, threshold_print =
PETSC_FALSE, flg, istranspose;
PetscBool silent = diffNorm != PETSC_NULLPTR ? PETSC_TRUE :
PETSC_FALSE;
@@ -2832,14 +2830,7 @@PetscErrorCode SNESTestJacobian(SNES snes,
PetscReal *Jnorm, PetscReal *diffNorm
PetscCall(MatComputeOperator(jacobian, MATAIJ, &A));
}
- PetscCall(MatGetType(A, &mattype));
- PetscCall(MatGetSize(A, &M, &N));
- PetscCall(MatGetLocalSize(A, &m, &n));
- PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
- PetscCall(MatSetType(B, mattype));
- PetscCall(MatSetSizes(B, m, n, M, N));
- PetscCall(MatSetBlockSizesFromMats(B, A, A));
- PetscCall(MatSetUp(B));
+ PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));
PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR,
PETSC_FALSE));
PetscCall(SNESGetFunction(snes, NULL, NULL, &functx));
@@ -2865,11 +2856,7 @@PetscErrorCode SNESTestJacobian(SNES snes,
PetscReal *Jnorm, PetscReal *diffNorm
const PetscInt *bcols;
const PetscScalar *bvals;
- PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
- PetscCall(MatSetType(C, mattype));
- PetscCall(MatSetSizes(C, m, n, M, N));
- PetscCall(MatSetBlockSizesFromMats(C, A, A));
- PetscCall(MatSetUp(C));
+ PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &C));
PetscCall(MatSetOption(C, MAT_NEW_NONZERO_ALLOCATION_ERR,
PETSC_FALSE));
PetscCall(MatAYPX(B, -1.0, A, DIFFERENT_NONZERO_PATTERN));
Thanks a lot!
Matteo
On 03/03/2026 03:28, Barry Smith wrote:
> Stefano hit the nail on the head!
>
> The matrices are of the same type it is just that MatView() for
> DMDA AIJ matrices MatView_MPI_DA() causes them to printed out in
> “natural ordering” related to the grid instead of PETSc parallel
> ordering. It is perhaps a case of Barry being too clever for his own
> good.
>
> SNESTestJacobian() has the code:
>
> ```
> PetscCall(MatGetType(A, &mattype));
> PetscCall(MatGetSize(A, &M, &N));
> PetscCall(MatGetLocalSize(A, &m, &n));
> PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
> PetscCall(MatSetType(B, mattype));
> PetscCall(MatSetSizes(B, m, n, M, N));
> PetscCall(MatSetBlockSizesFromMats(B, A, A));
> PetscCall(MatSetUp(B));
> PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR,
> PETSC_FALSE));
>
> PetscCall(SNESGetFunction(snes, NULL, NULL, &functx));
> PetscCall(SNESComputeJacobianDefault(snes, x, B, B, functx));
> ```
>
> It does not use MatDuplicate() hence the B (finite difference
> matrix) does not inherit the MatView_MPI_DA() function.
>
> Could you please try the following. Edit src/snes/interface/snes.c
> and locate the routine SNESTestJacobian() and replace the block of
> code that creates the B matrix above with
>
> ```
> PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &D));
> ```
>
> Then run make all in the PETSc home directory. Then rerun your
> program. The plan is now both matrices will be displayed in the same
> ordering (natural ordering on the multidimensional grid) and so the
> three matrices displayed should make sense.
>
> Barry
>
>
>
>
>> On Mar 2, 2026, at 6:00 PM, Matteo Semplice via petsc-users
>> <petsc-users at mcs.anl.gov> wrote:
>>
>> It's indeed a DMDA (with 2 dofs per point), but the code is run with
>>
>> -snes_test_jacobian -snes_test_jacobian_view
>>
>> Would the hand-coded and the f.d. jacobian be different matrix type
>> in this case?
>>
>> Matteo
>>
>> Il 02/03/26 21:32, Stefano Zampini ha scritto:
>>> Maybe you are using a DMDA, with one matrix coming from it and the
>>> other not? DMDA changes the View to dump in natural ordering.
>>>
>>> Il giorno lun 2 mar 2026 alle ore 18:10 Matteo Semplice via
>>> petsc-users <petsc-users at mcs.anl.gov> ha scritto:
>>>
>>> Dear all,
>>>
>>> the following output from -snes_test_jacobian_view lokks
>>> clear enough to pinpoint my mistake (I must be putting jacobian
>>> elements in the wrong place), but it nevertheless left me suprised.
>>>
>>> 0 SNES Function norm 2.945516399303e-01
>>> ---------- Testing Jacobian -------------
>>> Testing hand-coded Jacobian, if (for double precision runs)
>>> ||J - Jfd||_F/||J||_F is
>>> O(1.e-8), the hand-coded Jacobian is probably correct.
>>> ||J - Jfd||_F/||J||_F = 0.000121529, ||J - Jfd||_F = 0.395752
>>> Hand-coded Jacobian ----------
>>> Mat Object: 4 MPI processes
>>> row 84: (84, 60.) (85, 0.) (3446, -80.) (3447, 0.)
>>> (6808, 20.) (6809, 0.)
>>> row 85: (84, 0.00914434) (85, 1.) (3446, 0.)
>>> (3447, 0.) (6808, 0.) (6809, 0.)
>>>
>>> Finite difference Jacobian ----------
>>> Mat Object: 4 MPI processes
>>> type: mpiaij
>>> row 84: (84, 60.) (1806, -80.) (3528, 20.)
>>> row 85: (85, 0.999971)
>>>
>>> Hand-coded minus finite-difference Jacobian with tolerance
>>> 1e-05 ----------
>>> Mat Object: 4 MPI processes
>>> type: mpiaij
>>> row 84:
>>> row 85: (84, 0.00914434) (85, 2.9293e-05)
>>>
>>> The difference of the two matrices shoudl be nonzero in row 84
>>> for all elements 1806, 3528, 3446, 6808. Is it intended that the
>>> difference is shown only for matrix elements that are present in
>>> both the hand-coded and the FD jacobian?
>>>
>>> I would have expected an output like
>>>
>>> row 84: (1806,80.) (3446,-80.), etc
>>>
>>> (This was on PETSc Release Version 3.24.2)
>>>
>>> Thanks
>>>
>>> Matteo
>>>
>>> --
>>> Prof. Matteo Semplice
>>> Università degli Studi dell’Insubria
>>> Dipartimento di Scienza e Alta Tecnologia – DiSAT
>>> Professore Associato
>>> Via Valleggio, 11 – 22100 Como (CO) – Italia
>>> tel.: +39 031 2386316
>>>
>>>
>>>
>>> --
>>> Stefano
>> --
>> ---
>> Professore Associato in Analisi Numerica
>> Dipartimento di Scienza e Alta Tecnologia
>> Università degli Studi dell'Insubria
>> Via Valleggio, 11 - Como
>
--
Prof. Matteo Semplice
Università degli Studi dell’Insubria
Dipartimento di Scienza e Alta Tecnologia – DiSAT
Professore Associato
Via Valleggio, 11 – 22100 Como (CO) – Italia
tel.: +39 031 2386316
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260303/fc027b66/attachment-0001.html>
More information about the petsc-users
mailing list