Dear Barry, Thank you for the clear answer. Indeed, it would be useful to know that they are both stored in PC, that PCSetOperators may overwrite KSPSetOperators and that the order (1 for A, 2 for PC) is the order that appears in KSPView.
In any case, it is clear to me now. Thanks for the help! Cheers, Elena ________________________________ From: Barry Smith <[email protected]> Sent: 04 November 2025 16:08:21 To: Moral Sanchez, Elena Cc: PETSc Subject: Re: [petsc-users] norm of KSPBuildResidual does not match norm computed from KSPBuildSolution The preconditioner is always built from the second of the two matrices passed in KSPSetOperators() or PCSetOperators(). In the first case below both matrices are the same (the nest matrix) and so diagonal is extracted from the nest matrix (which is the only matrix). In the second case below the first matrix is custom and the second is the next matrix, again the preconditioner is constructed from the second one. linear system matrix = precond matrix: Mat Object: 1 MPI process type: nest rows=524, cols=524 Matrix object: type=nest, rows=3, cols=3 MatNest structure: (0,0) : type=mpiaij, rows=176, cols=176 (0,1) : NULL (0,2) : NULL (1,0) : NULL (1,1) : type=mpiaij, rows=172, cols=172 (1,2) : NULL (2,0) : NULL (2,1) : NULL (2,2) : type=mpiaij, rows=176, cols=176 Mat Object: 1 MPI process type: python Python: __main__.LHSOperator Mat Object: 1 MPI process type: nest Matrix object: type=nest, rows=3, cols=3 MatNest structure: (0,0) : type=mpiaij, rows=176, cols=176 (0,1) : NULL (0,2) : NULL (1,0) : NULL (1,1) : type=mpiaij, rows=172, cols=172 (1,2) : NULL (2,0) : NULL (2,1) : NULL (2,2) : type=mpiaij, rows=176, cols=176 The arguments to KSPSetOperators() and PCSetOperators() are stored in the same place (inside the PC) In fact, PetscErrorCode KSPSetOperators(KSP ksp, Mat Amat, Mat Pmat) { .... PetscCall(PCSetOperators(ksp->pc, Amat, Pmat)); so when you call precond.setOperators(A=nest_mass_matrix, P=None) ksp.setPC(precond) you are overwriting the A = lhs that you passed into ksp.setOperators(ksp, lhs, None). Given the names ksp.setOperators() and precond.setOperators() I can see how this can be confusing. It is reasonable to conclude that ksp.setOperators is providing the linear system and that pc.setOperators() is providing the matrix with which to build the preconditioner, but that is not the case. But the code has been this way for 31 years. I am not sure what we can change with the documentation. Perhaps in KSPSetOperators it can say that it sets them into the PC. Barry On Nov 4, 2025, at 4:06 AM, Moral Sanchez, Elena <[email protected]> wrote: Dear Barry, Thanks for the fast answer. Unfortunately in my case the discrepancy is huge. With the flags -ksp_monitor_true_residual -ksp_norm_type unpreconditioned this is the output: 0 KSP unpreconditioned resid norm 5.568889644229e-01 true resid norm 5.568889644229e-01 ||r(i)||/||b|| 1.000000000000e+00 1 KSP unpreconditioned resid norm 2.831772665189e-01 true resid norm 2.831772665189e-01 ||r(i)||/||b|| 5.084986139245e-01 2 KSP unpreconditioned resid norm 1.875950094147e-01 true resid norm 1.875950094147e-01 ||r(i)||/||b|| 3.368625011435e-01 and this is the output of my own monitor function: Iter 0/10 | res = 5.57e-01/1.00e-08 | 0.0 s difference KSPBuildSolution and u: 0.0 UNPRECONDITIONED norm: 0.5568889644229376 PRECONDITIONED norm: 2.049041078011257 KSPBuildResidual 2-norm: 0.5568889644229299 difference KSPBuildResidual and b-A(KSPBuildSolution): 6.573603152700697e-13 Iter 1/10 | res = 2.83e-01/1.00e-08 | 0.0 s difference KSPBuildSolution and u: 0.0 UNPRECONDITIONED norm: 0.7661983589104541 PRECONDITIONED norm: 2.7387602134717137 KSPBuildResidual 2-norm: 0.2831772665189212 difference KSPBuildResidual and b-A(KSPBuildSolution): 0.1700718741085172 Iter 2/10 | res = 1.88e-01/1.00e-08 | 0.0 s difference KSPBuildSolution and u: 0.0 UNPRECONDITIONED norm: 0.7050518160900253 PRECONDITIONED norm: 2.421773833445645 KSPBuildResidual 2-norm: 0.18759500941469456 difference KSPBuildResidual and b-A(KSPBuildSolution): 0.19327058976599623 Here u is the vector in the KSPSolve. 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. 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? 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. Thanks for the help. Cheers, Elena On 11/4/25 03:01, Barry Smith wrote: 0 KSP unpreconditioned resid norm 1.265943996096e+00 true resid norm 1.265943996096e+00 ||r(i)||/||b|| 1.000000000000e+00 My monitor 0 1.265943996096e+00
