"""Classical AMG (Ruge-Stuben AMG)."""
from warnings import warn
from scipy.sparse import csr_matrix, isspmatrix_csr, SparseEfficiencyWarning
import numpy as np
from pyamg.multilevel import MultilevelSolver
from pyamg.relaxation.smoothing import change_smoothers
from pyamg.strength import classical_strength_of_connection, \
symmetric_strength_of_connection, evolution_strength_of_connection, \
distance_strength_of_connection, energy_based_strength_of_connection, \
algebraic_distance, affinity_distance
from pyamg.classical.interpolate import direct_interpolation, classical_interpolation
from . import split
from .cr import CR
[docs]
def ruge_stuben_solver(A,
strength=('classical', {'theta': 0.25}),
CF=('RS', {'second_pass': False}),
interpolation='classical',
presmoother=('gauss_seidel', {'sweep': 'symmetric'}),
postsmoother=('gauss_seidel', {'sweep': 'symmetric'}),
max_levels=30, max_coarse=10, keep=False, **kwargs):
"""Create a multilevel solver using Classical AMG (Ruge-Stuben AMG).
Parameters
----------
A : csr_matrix
Square matrix in CSR format
strength : str
Valid strings are ['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 : str or tuple, default 'RS'
Method used for coarse grid selection (C/F splitting)
Supported methods are RS, PMIS, PMISc, CLJP, CLJPc, and CR.
interpolation : str, default 'classical'
Method for interpolation. Options include 'direct', 'classical'.
presmoother : str or dict
Method used for presmoothing at each level. Method-specific parameters
may be passed in using a tuple, e.g.
presmoother=('gauss_seidel',{'sweep':'symmetric}), the default.
postsmoother : str or dict
Postsmoothing method with the same usage as presmoother
max_levels : int, default 30
Maximum number of levels to be used in the multilevel solver.
max_coarse : int, default 20
Maximum number of variables permitted on the coarse grid.
keep : bool, default False
Flag to indicate keeping strength of connection (C) in the
hierarchy for diagnostics.
Returns
-------
ml : MultilevelSolver
Multigrid hierarchy of matrices and prolongation operators
Examples
--------
>>> from pyamg.gallery import poisson
>>> from pyamg import ruge_stuben_solver
>>> A = poisson((10,),format='csr')
>>> ml = ruge_stuben_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
----------
.. [2001TrOoSc] Trottenberg, U., Oosterlee, C. W., and Schuller, A.,
"Multigrid" San Diego: Academic Press, 2001. Appendix A
See Also
--------
aggregation.smoothed_aggregation_solver, MultilevelSolver,
aggregation.rootnode_solver
"""
levels = [MultilevelSolver.Level()]
# convert A to csr
if not isspmatrix_csr(A):
try:
A = csr_matrix(A)
warn('Implicit conversion of A to CSR',
SparseEfficiencyWarning)
except Exception as e:
raise TypeError('Argument A must have type csr_matrix, '
'or be convertible to csr_matrix') from e
# preprocess A
A = A.asfptype()
if A.shape[0] != A.shape[1]:
raise ValueError('expected square matrix')
levels[-1].A = A
while len(levels) < max_levels and levels[-1].A.shape[0] > max_coarse:
bottom = _extend_hierarchy(levels, strength, CF, interpolation, keep)
if bottom:
break
ml = MultilevelSolver(levels, **kwargs)
change_smoothers(ml, presmoother, postsmoother)
return ml
# internal function
def _extend_hierarchy(levels, strength, CF, interpolation, keep):
"""Extend the multigrid hierarchy."""
def unpack_arg(v):
if isinstance(v, tuple):
return v[0], v[1]
return v, {}
A = levels[-1].A
# 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 = split.RS(C, **kwargs)
elif fn == 'PMIS':
splitting = split.PMIS(C, **kwargs)
elif fn == 'PMISc':
splitting = split.PMISc(C, **kwargs)
elif fn == 'CLJP':
splitting = split.CLJP(C, **kwargs)
elif fn == 'CLJPc':
splitting = split.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
# Return early, do not add another coarse level
num_fpts = np.sum(splitting)
if (num_fpts == len(splitting)) or (num_fpts == 0):
return True
# 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)
else:
raise ValueError(f'Unknown interpolation method {interpolation}')
# Generate the restriction matrix that maps from the fine-grid to the
# coarse-grid
R = P.T.tocsr()
# 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
# Form next level through Galerkin product
levels.append(MultilevelSolver.Level())
A = R * A * P
levels[-1].A = A
return False
