[petsc-users] Diagnosing a difference between "unpreconditioned" and "true" residual norms

[Barry Smith]

   Patrick,

     I have only seen this when the "linear" operator turned out to not actually be linear or at least not linear in double precision. 
Are you using differencing or anything in your MatShell that might make it not be a linear operator in full precision?

Since your problem is so small you can compute the Jacobian explicitly via finite differencing and then use that matrix plus your shell preconditioner. I beat if you do this you will see the true and non-true residual norms remain the same, this would likely mean something is wonky with your shell matrix.

  Barry

> On Sep 9, 2016, at 9:32 AM, Patrick Sanan <patrick.sanan at gmail.com> wrote:
> 
> I am debugging a linear solver which uses a custom operator and
> preconditioner, via MATSHELL and PCSHELL. Convergence seems to be
> fine, except that I unexpectedly see a difference between the
> "unpreconditioned" and "true" residual norms when I use
> -ksp_monitor_true_residual with a right-preconditioned Krylov method
> (FGMRES or right-preconditioned GMRES).
> 
>  0 KSP unpreconditioned resid norm 9.266794204683e+08 true resid norm
> 9.266794204683e+08 ||r(i)||/||b|| 1.000000000000e+00
>  1 KSP unpreconditioned resid norm 2.317801431974e+04 true resid norm
> 2.317826550333e+04 ||r(i)||/||b|| 2.501217248530e-05
>  2 KSP unpreconditioned resid norm 4.453270507534e+00 true resid norm
> 2.699824780158e+01 ||r(i)||/||b|| 2.913439880638e-08
>  3 KSP unpreconditioned resid norm 1.015490793887e-03 true resid norm
> 2.658635801018e+01 ||r(i)||/||b|| 2.868991953738e-08
>  4 KSP unpreconditioned resid norm 4.710220776105e-07 true resid norm
> 2.658631616810e+01 ||r(i)||/||b|| 2.868987438467e-08
> KSP Object:(mgk_) 1 MPI processes
>  type: fgmres
>    GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
> Orthogonalization with no iterative refinement
>    GMRES: happy breakdown tolerance 1e-30
>  maximum iterations=10000, initial guess is zero
>  tolerances:  relative=1e-13, absolute=1e-50, divergence=10000.
>  right preconditioning
>  using UNPRECONDITIONED norm type for convergence test
> PC Object:(mgk_) 1 MPI processes
>  type: shell
>    Shell: Custom PC
>  linear system matrix = precond matrix:
>  Mat Object:  Custom Operator   1 MPI processes
>    type: shell
>    rows=256, cols=256
>      has attached null space
> 
> I have dumped the explicit operator and preconditioned operator, and I
> can see that the operator and the preconditioned operator each have a
> 1-dimensional nullspace (a constant-pressure nullspace) which I have
> accounted for by constructing a normalized, constant-pressure vector
> and supplying it to the operator via a MatNullSpace.
> 
> If I disregard the (numerically) zero singular value, the operator has
> a condition number of 1.5669e+05 and the preconditioned operator has a
> condition number of 1.01 (strong preconditioner).
> 
> Has anyone seen this sort of behavior before and if so, is there a
> common culprit that I am overlooking? Any ideas of what to test next
> to try to isolate the issue?
> 
> As I understand it, the unpreconditioned and true residual norms
> should be identical in exact arithmetic, so I would suspect that
> somehow I've ended up with a "bad Hessenberg matrix" in some way as I
> perform this solve (or maybe I have a more subtle bug).