how is snes_test_jacobian_view computing the matrix difference?
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
At first glance the output is surprising. I believe your understanding is correct and something is going very wrong to produce the output you report. Is there any way you can share the code that produces this output (if you send it to [email protected] the code will never be shared outside the PETSc developers who monitor that email). Barry
On Mar 2, 2026, at 10:09 AM, Matteo Semplice via petsc-users <[email protected]> wrote:
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
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 < [email protected]> 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
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 <[email protected]> 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
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 <[email protected]> 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 <[email protected] <mailto:[email protected]>> 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
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 <[email protected]> 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 <[email protected]> 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
Thanks for reporting the problem and testing the fix. I have an MR to fix the problem https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/9... Barry
On Mar 3, 2026, at 3:44 AM, Matteo Semplice via petsc-users <[email protected]> wrote:
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 <[email protected]> <mailto:[email protected]> 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 <[email protected] <mailto:[email protected]>> 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
participants (3)
-
Barry Smith -
Matteo Semplice -
Stefano Zampini