Probable bug in the Krylov solver

[436ahsan]

I have been testing AIR solver performance in some of my experiments from two different air branches :gen_air and master_air and it turns out after setting all the parameters same when I use AIR as pre-conditioner to GMRES, i.e. accelerated AIR solver then the solver results are different. For example:

Accelerated AIR solver results using gen_air branch
n | iter | rho | rhoasymp | OpCx | GdCx | Work

279000 | 10 | 0.1 | 0.11 | 1.2 | 1.9 | 1.2

Accelerated AIR solver results using master_air branch
n | iter | rho | rhoasymp | OpCx | GdCx | Work

279000 | 15 | 0.14 | 0.91 | 1.2 | 1.9 | 1.4

It turns out, this result differences only arises when I was using accel=‘gmres’ in the multilevel solver, so there could be a possible bug in the Krylov solvers. Standalone AIR solvers seems to be doing fine from both gen_air and master_air branches and the results are exactly same for both case.

[bensworth]

@lukeolson I havent done a diff on the krylov routines yet, but has anything changed in their implementation in recent years that you recall? This was a strange issue, hierarchies appeared to match exactly and Ahsan realized GMRES was the only obvious culprit. rhoasymp means final convergence factor, so note how the GMRES in main has slowed down a ton while in the old gen AIR branch is more or less the same as initial

[lukeolson]

There have been changes, including the stopping criteria.

One thing to try is the same experiment but with unpreconditioned gmres.

[jbschroder]

@436ahsan , nice work tracking this down.

Have you compared the residual norm histories between the two different GMRES? Are the residual norms different?

Then, Luke raises a good point. I would change your code to directly call GMRES

x = A, b, x0=None, tol=1e-5, restrt=None, maxiter=None, xtype=None,
M=None, callback=None, residuals=None, orthog='householder',
**kwargs):

instead of the solver_name.solve(...) interface that your using. The multigrid preconditioner M can be obtained with

M = solver_name.aspreconditioner()

You'll have to fill in x0, tol, maxiter, and so on to match your current experiments.

Then, you can run GMRES with no preconditioning. If there is a difference between the two GMRES, this will be a good way to track it down.

[436ahsan]

@jbschroder In my earlier tests, yes, my residuals norms are different. I am now looking into the suggested experiments and keep you all updated on this issue.

[bensworth]

@436ahsan , can you try again on the air_pr branch and gen_air branch? If i copy your test and run

from pyamg.classical import air_solver as air
# from pyamg.classical import AIR_solver as air
from pyamg.gallery import poisson
from scipy.sparse.linalg import gmres
import numpy as np
A = poisson((100, 100), format='csr')     # matrix
rng = np.random.default_rng(2021)
b = rng.random((A.shape[0],))
ml = air(A, presmoother=None, postsmoother=('fc_jacobi', {'omega': 1.0, 'withrho': False, 'iterations': 1, 'f_iterations': 2, 'c_iterations': 0} ))
# ml = air(A, presmoother=None, postsmoother=('FC_jacobi', {'omega': 1.0, 'withrho': False, 'iterations': 1, 'F_iterations': 2, 'C_iterations': 0} ))
M = ml.aspreconditioner(cycle='V')       # preconditioner
x, info = gmres(A, b, tol=1e-8, maxiter=30, M=M) # solve with CG
x

on these two branches I get identical results.

Also, I manually set the postsmoother due to different defaults in the two branches. Here if I use accel in pyamg or external, I still get more or les identical results.