<!DOCTYPE html><html><head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
</head>
<body>
<p>Dear Barry,</p>
<p> indeed your suggestion fixes things! (I assume that you
meant PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));)</p>
<p>Now I get a consistent output.</p>
<p>Hand-coded Jac:</p>
<p><span style="font-family:monospace"><span style="color:#000000;background-color:#ffffff;"> row 84:
(84, 60.) (85, 0.) (3446, -80.) (3447, 0.)
(6808, 20.) (6809, 0.) </span><span style="color:#000000;background-color:#ffffff;"> </span><br>
<span style="color:#000000;background-color:#ffffff;"> row
85: (84, 0.00914434) (85, 1.) (3446, 0.)
(3447, 0.) (6808, 0.) (6809, 0.) </span><br>
<span style="color:#000000;background-color:#ffffff;"></span><br>
</span></p>
<p>Finite-difference Jac:</p>
<p><span style="font-family:monospace"><span style="color:#000000;background-color:#ffffff;"> row 84:
(84, 60.) (85, 0.) (3446, -80.) (3447, 0.)
(6808, 20.) (6809, 0.) </span><span style="color:#000000;background-color:#ffffff;"> </span><br>
<span style="color:#000000;background-color:#ffffff;"> row
85: (84, 0.) (85, 0.999971) (3446, 0.)
(3447, 0.) (6808, 0.) (6809, 0.) </span><br>
</span></p>
<p>Difference of the two matrices:<br>
<span style="font-family:monospace"><span style="color:#000000;background-color:#ffffff;"> row 84:
(84, 0.) (85, 0.) (3446, 0.) (3447, 0.)
(6808, 0.) (6809, 0.) </span><span style="color:#000000;background-color:#ffffff;"> </span><br>
<span style="color:#000000;background-color:#ffffff;"> row
85: (84, 0.00914434) (85, 2.9293e-05) (3446, 0.)
(3447, 0.) (6808, 0.) (6809, 0.) </span><br>
</span>
<br>
</p>
<p>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:</p>
<p><span style="font-family:monospace"><span style="font-weight:bold;color:#000000;background-color:#ffffff;">diff
--git a/src/snes/interface/snes.c b/src/snes/interface/snes.c</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="font-weight:bold;color:#000000;background-color:#ffffff;">index
f66c23a1417..eee3d070d4c 100644</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="font-weight:bold;color:#000000;background-color:#ffffff;">---
a/src/snes/interface/snes.c</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="font-weight:bold;color:#000000;background-color:#ffffff;">+++
b/src/snes/interface/snes.c</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#18b2b2;background-color:#ffffff;">@@ -2768,8
+2768,6 @@</span><span style="color:#000000;background-color:#ffffff;">
PetscErrorCode SNESTestJacobian(SNES snes, PetscReal *Jnorm,
PetscReal *diffNorm</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;"> Vec
x = snes->vec_sol, f;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;">
PetscReal nrm, gnorm;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;">
PetscReal threshold = 1.e-5;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">- MatType
mattype;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscInt m, n, M, N;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;"> void
*functx;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;">
PetscBool complete_print = PETSC_FALSE,
threshold_print = PETSC_FALSE, flg, istranspose;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;">
PetscBool silent = diffNorm != PETSC_NULLPTR ?
PETSC_TRUE : PETSC_FALSE;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#18b2b2;background-color:#ffffff;">@@
-2832,14 +2830,7 @@</span><span style="color:#000000;background-color:#ffffff;">
PetscErrorCode SNESTestJacobian(SNES snes, PetscReal *Jnorm,
PetscReal *diffNorm</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;">
PetscCall(MatComputeOperator(jacobian, MATAIJ, &A));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;"> }</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;"> </span><span style="color:#000000;background-color:#ffffff;"> </span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatGetType(A, &mattype));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatGetSize(A, &M, &N));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatGetLocalSize(A, &m, &n));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatCreate(PetscObjectComm((PetscObject)A),
&B));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatSetType(B, mattype));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatSetSizes(B, m, n, M, N));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatSetBlockSizesFromMats(B, A, A));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatSetUp(B));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#18b218;background-color:#ffffff;">+
PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;">
PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR,
PETSC_FALSE));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;"> </span><span style="color:#000000;background-color:#ffffff;"> </span><br>
<span style="color:#000000;background-color:#ffffff;">
PetscCall(SNESGetFunction(snes, NULL, NULL, &functx));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#18b2b2;background-color:#ffffff;">@@
-2865,11 +2856,7 @@</span><span style="color:#000000;background-color:#ffffff;">
PetscErrorCode SNESTestJacobian(SNES snes, PetscReal *Jnorm,
PetscReal *diffNorm</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;">
const PetscInt *bcols;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;">
const PetscScalar *bvals;</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;"> </span><span style="color:#000000;background-color:#ffffff;"> </span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatCreate(PetscObjectComm((PetscObject)A),
&C));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatSetType(C, mattype));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatSetSizes(C, m, n, M, N));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatSetBlockSizesFromMats(C, A, A));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#b21818;background-color:#ffffff;">-
PetscCall(MatSetUp(C));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#18b218;background-color:#ffffff;">+
PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES,
&C));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;">
PetscCall(MatSetOption(C,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));</span><span style="color:#000000;background-color:#ffffff;">
</span><br>
<span style="color:#000000;background-color:#ffffff;"> </span><span style="color:#000000;background-color:#ffffff;"> </span><br>
<span style="color:#000000;background-color:#ffffff;">
PetscCall(MatAYPX(B, -1.0, A,
DIFFERENT_NONZERO_PATTERN));</span><br>
<span style="color:#000000;background-color:#ffffff;">
</span><br>
</span>
<br>
</p>
<p>Thanks a lot!</p>
<p>Matteo</p>
<div class="moz-cite-prefix">On 03/03/2026 03:28, Barry Smith wrote:<br>
</div>
<blockquote type="cite" cite="mid:C5164AE4-27F4-4822-BEFB-05B5F35064BA@petsc.dev">
<div> Stefano hit the nail on the head!</div>
<div><br>
</div>
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.
<div><br>
</div>
<div> SNESTestJacobian() has the code:</div>
<div><br>
</div>
<div>```</div>
<div> PetscCall(MatGetType(A, &mattype));
<div> PetscCall(MatGetSize(A, &M, &N));</div>
<div> PetscCall(MatGetLocalSize(A, &m, &n));</div>
<div> PetscCall(MatCreate(PetscObjectComm((PetscObject)A),
&B));</div>
<div> PetscCall(MatSetType(B, mattype));</div>
<div> PetscCall(MatSetSizes(B, m, n, M, N));</div>
<div> PetscCall(MatSetBlockSizesFromMats(B, A, A));</div>
<div> PetscCall(MatSetUp(B));</div>
<div> PetscCall(MatSetOption(B,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));</div>
<div><br>
</div>
<div> PetscCall(SNESGetFunction(snes, NULL, NULL,
&functx));</div>
<div> PetscCall(SNESComputeJacobianDefault(snes, x, B, B,
functx));</div>
<div>```</div>
<div><br>
</div>
<div> It does not use MatDuplicate() hence the B (finite
difference matrix) does not inherit the MatView_MPI_DA()
function.</div>
<div><br>
</div>
<div> 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 </div>
<div><br>
</div>
<div>```</div>
<div> PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES,
&D));</div>
<div>```</div>
<div><br>
</div>
<div> 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.</div>
<div><br>
</div>
<div> Barry</div>
<div><br>
</div>
<div><br>
</div>
<div><br id="lineBreakAtBeginningOfMessage">
<div><br>
<blockquote type="cite">
<div>On Mar 2, 2026, at 6:00 PM, Matteo Semplice via
petsc-users <a class="moz-txt-link-rfc2396E" href="mailto:petsc-users@mcs.anl.gov"><petsc-users@mcs.anl.gov></a> wrote:</div>
<br class="Apple-interchange-newline">
<div>
<div>
<p>It's indeed a DMDA (with 2 dofs per point), but the
code is run with</p>
<p>-snes_test_jacobian -snes_test_jacobian_view</p>
<p>Would the hand-coded and the f.d. jacobian be
different matrix type in this case?</p>
<p>Matteo</p>
<div class="moz-cite-prefix">Il 02/03/26 21:32,
Stefano Zampini ha scritto:<br>
</div>
<blockquote type="cite" cite="mid:CAGPUisiMN2k2A-2VVhokn6uaDxmDYFcpgZ_cMka=z3gU8Wj-XQ@mail.gmail.com">
<div dir="ltr">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. </div>
<br>
<div class="gmail_quote gmail_quote_container">
<div dir="ltr" class="gmail_attr">Il giorno lun 2
mar 2026 alle ore 18:10 Matteo Semplice via
petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" moz-do-not-send="true" class="moz-txt-link-freetext">petsc-users@mcs.anl.gov</a>>
ha scritto:<br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>
<p>Dear all,</p>
<p> 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.</p>
<font face="monospace">0 SNES Function norm
2.945516399303e-01<br>
---------- Testing Jacobian -------------<br>
Testing hand-coded Jacobian, if (for
double precision runs) ||J - Jfd||_F/||J||_F
is<br>
O(1.e-8), the hand-coded Jacobian is
probably correct.<br>
||J - Jfd||_F/||J||_F = 0.000121529, ||J -
Jfd||_F = 0.395752<br>
Hand-coded Jacobian ----------<br>
Mat Object: 4 MPI processes<br>
row 84: (84, 60.) (85, 0.) <font color="#ed5c57">(3446, -80.) </font>
(3447, 0.) <font color="#ed5c57">
(6808, 20.) </font> (6809, 0.) <br>
row 85: (84, 0.00914434) (85,
1.) (3446, 0.) (3447, 0.)
(6808, 0.) (6809, 0.) <br>
<br>
Finite difference Jacobian ----------<br>
Mat Object: 4 MPI processes<br>
type: mpiaij<br>
row 84: (84, 60.) <font color="#ed5c57">(1806, -80.) </font> <font color="#ed5c57"> (3528, 20.) </font> <br>
row 85: (85, 0.999971)<br>
<br>
Hand-coded minus finite-difference
Jacobian with tolerance 1e-05 ----------<br>
Mat Object: 4 MPI processes<br>
type: mpiaij<br>
<font color="#ed5c57">row 84:</font> <br>
row 85: (84, 0.00914434) (85,
2.9293e-05)</font>
<p>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?</p>
<p>I would have expected an output like</p>
<p>row 84: (1806,80.) (3446,-80.), etc</p>
<p>(This was on <span style="font-family:monospace"><span style="background-color: rgb(255, 255, 255);">PETSc Release Version
3.24.2</span>)</span> </p>
<p>Thanks</p>
<p> Matteo</p>
<pre cols="72">--
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</pre>
</div>
</blockquote>
</div>
<div><br clear="all">
</div>
<div><br>
</div>
<span class="gmail_signature_prefix">-- </span><br>
<div dir="ltr" class="gmail_signature">Stefano</div>
</blockquote>
<pre class="moz-signature" cols="72">--
---
Professore Associato in Analisi Numerica
Dipartimento di Scienza e Alta Tecnologia
Università degli Studi dell'Insubria
Via Valleggio, 11 - Como</pre>
</div>
</div>
</blockquote>
</div>
<br>
</div>
</div>
</blockquote>
<pre class="moz-signature" cols="72">--
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</pre>
</body>
</html>