"""Approximate idealkrestriction AMG."""
from copy import deepcopy
import numpy as np
from scipy.sparse import isspmatrix_csr, isspmatrix_bsr
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
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_matrix
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 : {string} : default 'one_point'
Options include 'direct', 'classical', 'inject' and 'one-point'.
restrict : {string} : default distance-2 AIR, with theta = 0.05.
Option is 'air' for local approximate ideal restriction (lAIR),
with inner options specifying degree, strength tolerance, etc.
presmoother : {string or dict} : default None
Method used for presmoothing at each level. Method-specific parameters
may be passed in using a tuple.
postsmoother : {string or dict}
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} : default False
Flag to indicate keeping strength of connection matrix (C) in
hierarchy.
Returns
-------
ml : multilevel_solver
Multigrid hierarchy of matrices and prolongation operators
Examples
--------
>>> from pyamg.gallery import poisson
>>> from pyamg import air_solver
>>> A = poisson((10,),format='csr')
>>> ml = air_solver(A,max_coarse=3)
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', 'lu', ...] or a callable
function, and args is a dictionary of arguments to be passed to fn.
See [2001TrOoSc]_ 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.
See Also
--------
aggregation.smoothed_aggregation_solver, multilevel_solver,
aggregation.rootnode_solver, ruge_stuben_solver
"""
# preprocess A
A = A.asfptype()
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 (isspmatrix_csr(P) and (not isspmatrix_csr(A))):
A = A.tocsr()
elif (isspmatrix_bsr(P) and (not isspmatrix_bsr(A))):
A = A.tobsr()
levels.append(MultilevelSolver.Level())
levels[-1].A = A
return 0
