"""Approximate idealkrestriction AMG."""
from copy import deepcopy
import numpy as np
from scipy.sparse import issparse
from ..multilevel import MultilevelSolver
from ..relaxation.smoothing import change_smoothers
from ..strength import (classical_strength_of_connection,
symmetric_strength_of_connection, evolution_strength_of_connection,
distance_strength_of_connection, algebraic_distance,
affinity_distance, energy_based_strength_of_connection)
from ..util.utils import filter_matrix_rows, asfptype
from ..classical.interpolate import (direct_interpolation, classical_interpolation,
injection_interpolation, one_point_interpolation,
local_air)
from .split import RS, PMIS, PMISc, CLJP, CLJPc
from .cr import CR
[docs]
def air_solver(A,
strength=('classical', {'theta': 0.3, 'norm': 'min'}),
CF=('RS', {'second_pass': True}),
interpolation='one_point',
restrict=('air', {'theta': 0.05, 'degree': 2}),
presmoother=None,
postsmoother=('fc_jacobi', {'omega': 1.0, 'iterations': 1,
'withrho': False, 'f_iterations': 2,
'c_iterations': 1}),
filter_operator=None,
max_levels=20, max_coarse=20,
keep=False, **kwargs):
"""Create a multilevel solver using approximate ideal restriction (AIR) AMG.
Parameters
----------
A : csr_array
Square (non)symmetric matrix in CSR format.
strength : {'symmetric', 'classical', 'evolution', 'distance', \
'algebraic_distance','affinity', 'energy_based', None}
Method used to determine the strength of connection between unknowns
of the linear system. Method-specific parameters may be passed in
using a tuple, e.g. strength=('symmetric',{'theta' : 0.25 }). If
strength=None, all nonzero entries of the matrix are considered strong.
CF : {string} : default 'RS' with second pass
Method used for coarse grid selection (C/F splitting)
Supported methods are RS, PMIS, PMISc, CLJP, CLJPc, and CR.
interpolation : str
Options include 'direct', 'classical', 'inject' and 'one-point'.
restrict : str
Option is 'air' for local approximate ideal restriction (lAIR),
with inner options specifying degree, strength tolerance, etc..
presmoother : str
Method used for presmoothing at each level. Method-specific parameters
may be passed in using a tuple.
postsmoother : str
Postsmoothing method with the same usage as presmoother.
postsmoother=('fc_jacobi', ... ) with 2 F-sweeps, 1 C-sweep is default.
filter_operator : (bool, tol) : default None
Remove small entries in operators on each level if True. Entries are
considered *small* if `|a_ij| < tol |a_ii|`.
max_levels : {integer} : default 20
Maximum number of levels to be used in the multilevel solver.
max_coarse : {integer} : default 20
Maximum number of variables permitted on the coarse grid.
keep : bool
Flag to indicate keeping strength of connection matrix (C) in hierarchy.
**kwargs : dict
Extra keywords passed to the Multilevel class.
Returns
-------
MultilevelSolver
Multigrid hierarchy of matrices and prolongation operators.
See Also
--------
pyamg.multilevel.MultilevelSolver,
pyamg.aggregation.smoothed_aggregation_solver,
pyamg.aggregation.rootnode_solver,
ruge_stuben_solver
Notes
-----
`coarse_solver` is an optional argument and is the solver used at the
coarsest grid. The default is a pseudo-inverse. Most simply,
coarse_solver can be one of {'splu', 'lu', 'cholesky, 'pinv',
'gauss_seidel'}. Additionally, `coarse_solver` may be a tuple
`(fn, args)`, where `fn` is a string such as 'splu' or a callable
function, and args is a dictionary of arguments to be passed to fn.
See [3]_ for additional details.
References
----------
.. [1] Manteuffel, T. A., Münzenmaier, S., Ruge, J., & Southworth, B. S.
(2019). Nonsymmetric reduction-based algebraic multigrid. SIAM
Journal on Scientific Computing, 41(5), S242-S268.
.. [2] Manteuffel, T. A., Ruge, J., & Southworth, B. S. (2018).
Nonsymmetric algebraic multigrid based on local approximate ideal
restriction (lAIR). SIAM Journal on Scientific Computing, 40(6),
A4105-A4130.
.. [3] Trottenberg, U.; Oosterlee, C. W. & Schüller, A. (2001),
Multigrid, Vol. 33, Academic Press.
Examples
--------
>>> from pyamg.gallery import poisson
>>> from pyamg import air_solver
>>> A = poisson((10,),format='csr')
>>> ml = air_solver(A,max_coarse=3)
"""
# preprocess A
A = asfptype(A)
if A.shape[0] != A.shape[1]:
raise ValueError('expected square matrix')
if np.iscomplexobj(A.data):
raise ValueError('AIR solver not verified for complex matrices')
levels = [MultilevelSolver.Level()]
levels[-1].A = A
while len(levels) < max_levels and levels[-1].A.shape[0] > max_coarse:
bottom = extend_hierarchy(levels, strength, CF, interpolation, restrict,
filter_operator, keep)
if bottom:
break
ml = MultilevelSolver(levels, **kwargs)
change_smoothers(ml, presmoother, postsmoother)
return ml
def extend_hierarchy(levels, strength, CF, interpolation, restrict, filter_operator, keep):
"""Extend the multigrid hierarchy."""
def unpack_arg(v):
if isinstance(v, tuple):
return v[0], v[1]
return v, {}
# Filter operator. Need to keep original matrix on finest level for
# computing residuals
if (filter_operator is not None) and (filter_operator[1] != 0):
if len(levels) == 1:
A = deepcopy(levels[-1].A)
else:
A = levels[-1].A
filter_matrix_rows(A, filter_operator[1], diagonal=True, lump=filter_operator[0])
else:
A = levels[-1].A
# Check if matrix was filtered to be diagonal --> coarsest grid
if A.nnz == A.shape[0]:
return 1
# Compute the strength-of-connection matrix C, where larger
# C[i,j] denote stronger couplings between i and j.
fn, kwargs = unpack_arg(strength)
if fn == 'symmetric':
C = symmetric_strength_of_connection(A, **kwargs)
elif fn == 'classical':
C = classical_strength_of_connection(A, **kwargs)
elif fn == 'distance':
C = distance_strength_of_connection(A, **kwargs)
elif fn in ('ode', 'evolution'):
C = evolution_strength_of_connection(A, **kwargs)
elif fn == 'energy_based':
C = energy_based_strength_of_connection(A, **kwargs)
elif fn == 'algebraic_distance':
C = algebraic_distance(A, **kwargs)
elif fn == 'affinity':
C = affinity_distance(A, **kwargs)
elif fn is None:
C = A
else:
raise ValueError(f'Unrecognized strength of connection method: {fn}')
# Generate the C/F splitting
fn, kwargs = unpack_arg(CF)
if fn == 'RS':
splitting = RS(C, **kwargs)
elif fn == 'PMIS':
splitting = PMIS(C, **kwargs)
elif fn == 'PMISc':
splitting = PMISc(C, **kwargs)
elif fn == 'CLJP':
splitting = CLJP(C, **kwargs)
elif fn == 'CLJPc':
splitting = CLJPc(C, **kwargs)
elif fn == 'CR':
splitting = CR(C, **kwargs)
else:
raise ValueError(f'Unknown C/F splitting method {CF}')
# Make sure all points were not declared as C- or F-points
num_fpts = np.sum(splitting)
if (num_fpts == len(splitting)) or (num_fpts == 0):
return 1
# Generate the interpolation matrix that maps from the coarse-grid to the
# fine-grid
fn, kwargs = unpack_arg(interpolation)
if fn == 'classical':
P = classical_interpolation(A, C, splitting, **kwargs)
elif fn == 'direct':
P = direct_interpolation(A, C, splitting, **kwargs)
elif fn == 'one_point':
P = one_point_interpolation(A, C, splitting, **kwargs)
elif fn == 'inject':
P = injection_interpolation(A, splitting, **kwargs)
else:
raise ValueError(f'Unknown interpolation method {fn}')
# Build restriction operator
fn, kwargs = unpack_arg(restrict)
if fn == 'air':
R = local_air(A, splitting, **kwargs)
else:
raise ValueError(f'Unknown restriction method {fn}')
# Store relevant information for this level
if keep:
levels[-1].C = C # strength of connection matrix
levels[-1].splitting = splitting.astype(bool) # C/F splitting
levels[-1].P = P # prolongation operator
levels[-1].R = R # restriction operator
# RAP = R*(A*P)
A = R @ A @ P
# Make sure coarse-grid operator is in correct sparse format
if issparse(P) and P.format == 'csr' and issparse(A) and A.format != 'csr':
A = A.tocsr()
elif issparse(P) and P.format == 'bsr' and issparse(A) and A.format != 'bsr':
A = A.tobsr()
levels.append(MultilevelSolver.Level())
levels[-1].A = A
return 0
