<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="display:none;"><!-- P {margin-top:0;margin-bottom:0;} --></style>
</head>
<body dir="ltr">
<style type="text/css" style="display:none;"><!-- P {margin-top:0;margin-bottom:0;} --></style>
<div id="divtagdefaultwrapper" style="font-size:12pt;color:#000000;font-family:Calibri,Helvetica,sans-serif;" dir="ltr">
<p></p>
<div>Dear Barry,<br>
<br>
I just realized that the operator from KSPGetOperators()[0] does not match my operator. In fact, KSPGetOperators()[0] returns the operator used for the preconditioner. This explains all the issues I reported.<br>
<br>
The way I was setting up the KSP is<br>
ksp.setOperators(ksp, lhs, None)<br>
<br>
precond = PETSc.PC().create(comm=comm)<br>
precond.setType(PETSc.PC.Type.JACOBI)<br>
precond.setOperators(A=nest_mass_matrix, P=None)<br>
precond.setUp()<br>
<br>
ksp.setPC(precond)<br>
ksp.setUp()<br>
<br>
This turns out to behave as<br>
ksp.setOperators(A=nest_mass_matrix, P=nest_mass_matrix)<br>
precond = ksp.getPC()<br>
precond.setType(PETSc.PC.Type.JACOBI)<br>
ksp.setUp()<br>
<br>
I managed to set up what I want with<br>
ksp.setOperators(A=lhs, P=nest_mass_matrix)<br>
precond = ksp.getPC()<br>
precond.setType(PETSc.PC.Type.JACOBI)<br>
ksp.setUp()<br>
<br>
Here I am solving </div>
<div><span> </span>lhs x = rhs, </div>
<div>lhs is a matrix-free operator (python type) and my preconditioner is the diagonal of nest_mass_matrix, which is of type nest.<br>
<br>
I find this extremely confusing, especially because from the output of KSPView I could not detect the problem. For illustration, after the fix, PCView prints<br>
<br>
Mat Object: 1 MPI process<br>
type: python<br>
Python: __main__.LHSOperator<br>
Mat Object: 1 MPI process<br>
type: nest<br>
Matrix object:<br>
type=nest, rows=3, cols=3<br>
MatNest structure:<br>
(0,0) : type=mpiaij, rows=176, cols=176<br>
(0,1) : NULL<br>
(0,2) : NULL<br>
(1,0) : NULL<br>
(1,1) : type=mpiaij, rows=172, cols=172<br>
(1,2) : NULL<br>
(2,0) : NULL<br>
(2,1) : NULL<br>
(2,2) : type=mpiaij, rows=176, cols=176<br>
<br>
This brings me to the question: is the Jacobi preconditioner using the operator lhs or the nested operator? I think that the output of PCView and KSPView should be more clear on this.<br>
<br>
Cheers,<br>
Elena</div>
<br>
<p></p>
</div>
<hr style="display:inline-block;width:98%" tabindex="-1">
<div id="divRplyFwdMsg" dir="ltr"><font face="Calibri, sans-serif" style="font-size:11pt" color="#000000"><b>From:</b> Moral Sanchez, Elena<br>
<b>Sent:</b> 04 November 2025 10:06:07<br>
<b>To:</b> Barry Smith<br>
<b>Cc:</b> PETSc<br>
<b>Subject:</b> Re: [petsc-users] norm of KSPBuildResidual does not match norm computed from KSPBuildSolution</font>
<div> </div>
</div>
<div>
<div id="divtagdefaultwrapper" style="font-size:12pt;color:#000000;font-family:Calibri,Helvetica,sans-serif;" dir="ltr">
<p></p>
<div>
<p>Dear Barry,</p>
<p>Thanks for the fast answer. Unfortunately in my case the discrepancy is huge. With the flags</p>
<blockquote>
<p> -ksp_monitor_true_residual -ksp_norm_type unpreconditioned</p>
</blockquote>
<p>this is the output:</p>
<blockquote>
<p> 0 KSP unpreconditioned resid norm 5.568889644229e-01 true resid norm 5.568889644229e-01 ||r(i)||/||b|| 1.000000000000e+00<br>
1 KSP unpreconditioned resid norm 2.831772665189e-01 true resid norm 2.831772665189e-01 ||r(i)||/||b|| 5.084986139245e-01<br>
2 KSP unpreconditioned resid norm 1.875950094147e-01 true resid norm 1.875950094147e-01 ||r(i)||/||b|| 3.368625011435e-01</p>
</blockquote>
<p>and this is the output of my own monitor function:</p>
<blockquote>
<p>Iter 0/10 | res = 5.57e-01/1.00e-08 | 0.0 s<br>
difference KSPBuildSolution and u: 0.0<br>
UNPRECONDITIONED norm: 0.5568889644229376<br>
PRECONDITIONED norm: 2.049041078011257<br>
KSPBuildResidual 2-norm: 0.5568889644229299<br>
difference KSPBuildResidual and b-A(KSPBuildSolution): 6.573603152700697e-13<br>
<br>
Iter 1/10 | res = 2.83e-01/1.00e-08 | 0.0 s<br>
difference KSPBuildSolution and u: 0.0<br>
UNPRECONDITIONED norm: 0.7661983589104541<br>
PRECONDITIONED norm: 2.7387602134717137<br>
KSPBuildResidual 2-norm: 0.2831772665189212<br>
difference KSPBuildResidual and b-A(KSPBuildSolution): 0.1700718741085172<br>
<br>
Iter 2/10 | res = 1.88e-01/1.00e-08 | 0.0 s<br>
difference KSPBuildSolution and u: 0.0<br>
UNPRECONDITIONED norm: 0.7050518160900253<br>
PRECONDITIONED norm: 2.421773833445645<br>
KSPBuildResidual 2-norm: 0.18759500941469456<br>
difference KSPBuildResidual and b-A(KSPBuildSolution): 0.19327058976599623</p>
</blockquote>
<p>Here u is the vector in the KSPSolve. </p>
<p>After the first iteration, the residual computed from KSPBuildSolution and the residual from KSPBuildResidual diverge. They are the same when I run the same code without preconditioner.</p>
<p>Another observation is that after convergence (wrt. unpreconditioned norm == 2-norm of KSPBuildResidual) the solution with and without preconditioner looks quite different. How is this possible if my preconditioner is SPD?</p>
<p><br>
</p>
<p>By the way, where can I find your implementation of "My monitor" in src/snes/tutorials/ex5.c? I tried to look at the Gitlab repository but could not find it.</p>
<p>Thanks for the help.</p>
<p>Cheers,</p>
<p>Elena</p>
</div>
<p></p>
<p><br>
</p>
<p><br>
</p>
<p><br>
</p>
<p><br>
</p>
<p><br>
</p>
<div class="moz-cite-prefix">On 11/4/25 03:01, Barry Smith wrote:<br>
</div>
<blockquote type="cite">
<div> 0 KSP unpreconditioned resid norm 1.265943996096e+00 true resid norm 1.265943996096e+00 ||r(i)||/||b|| 1.000000000000e+00</div>
<div>My monitor 0 1.265943996096e+00</div>
</blockquote>
</div>
</div>
</body>
</html>