[petsc-users] norm of KSPBuildResidual does not match norm computed from KSPBuildSolution

[Moral Sanchez, Elena]

Dear Barry,

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.

The way I was setting up the KSP is
    ksp.setOperators(ksp, lhs, None)

    precond = PETSc.PC().create(comm=comm)
    precond.setType(PETSc.PC.Type.JACOBI)
    precond.setOperators(A=nest_mass_matrix, P=None)
    precond.setUp()

    ksp.setPC(precond)
    ksp.setUp()

This turns out to behave as
    ksp.setOperators(A=nest_mass_matrix, P=nest_mass_matrix)
    precond = ksp.getPC()
    precond.setType(PETSc.PC.Type.JACOBI)
    ksp.setUp()

I managed to set up what I want with
    ksp.setOperators(A=lhs, P=nest_mass_matrix)
    precond = ksp.getPC()
    precond.setType(PETSc.PC.Type.JACOBI)
    ksp.setUp()

Here I am solving
    lhs x = rhs,
lhs is a matrix-free operator (python type) and my preconditioner is the diagonal of nest_mass_matrix, which is of type nest.

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

    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

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.

Cheers,
Elena


________________________________
From: Moral Sanchez, Elena
Sent: 04 November 2025 10:06:07
To: Barry Smith
Cc: PETSc
Subject: Re: [petsc-users] norm of KSPBuildResidual does not match norm computed from KSPBuildSolution


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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20251104/544104cc/attachment-0001.html>