Source code for pyamg.classical.classical

"""Classical AMG (Ruge-Stuben AMG)"""

__docformat__ = "restructuredtext en"

from warnings import warn
import scipy
from scipy.sparse import csr_matrix, isspmatrix_csr, isspmatrix_bsr, eye

from pyamg.multilevel import multilevel_solver
from pyamg.relaxation.smoothing import change_smoothers
from pyamg.strength import classical_strength_of_connection, \
                           symmetric_strength_of_connection, \
                           evolution_strength_of_connection

from interpolate import direct_interpolation

__all__ = ['ruge_stuben_solver']

[docs]def ruge_stuben_solver(A, strength=('classical',{'theta':0.25}), CF='RS', presmoother=('gauss_seidel',{'sweep':'symmetric'}), postsmoother=('gauss_seidel',{'sweep':'symmetric'}), max_levels=10, max_coarse=500, keep=False, **kwargs): """Create a multilevel solver using Classical AMG (Ruge-Stuben AMG) Parameters ---------- A : csr_matrix Square matrix in CSR format strength : ['symmetric', 'classical', 'evolution', 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' Method used for coarse grid selection (C/F splitting) Supported methods are RS, PMIS, PMISc, CLJP, and CLJPc presmoother : {string 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 : {string or dict} Postsmoothing method with the same usage as presmoother max_levels: {integer} : default 10 Maximum number of levels to be used in the multilevel solver. max_coarse: {integer} : default 500 Maximum number of variables permitted on the coarse grid. keep: {bool} : default False Flag to indicate keeping extra operators in the hierarchy for diagnostics. For example, if True, then strength of connection (C) and tentative prolongation (T) are kept. Returns ------- ml : multilevel_solver 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. References ---------- .. [1] Trottenberg, U., Oosterlee, C. W., and Schuller, A., "Multigrid" San Diego: Academic Press, 2001. Appendix A See Also -------- aggregation.smoothed_aggregation_solver, multilevel_solver, aggregation.rootnode_solver """ levels = [ multilevel_solver.level() ] # convert A to csr if not isspmatrix_csr(A): try: A = csr_matrix(A) warn("Implicit conversion of A to CSR", scipy.sparse.SparseEfficiencyWarning) except: raise TypeError('Argument A must have type csr_matrix, or be convertible to csr_matrix') # 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: extend_hierarchy(levels, strength, CF, keep) ml = multilevel_solver(levels, **kwargs) change_smoothers(ml, presmoother, postsmoother) return ml # internal function
def extend_hierarchy(levels, strength, CF, keep): """ helper function for local methods """ def unpack_arg(v): if isinstance(v,tuple): return v[0],v[1] else: return v,{} A = levels[-1].A # strength of connection 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 == 'ode') or (fn == 'evolution'): raise NotImplementedError('evolution method not supported for Classical AMG') elif fn is None: C = A else: raise ValueError('unrecognized strength of connection method: %s' % fn) if CF in [ 'RS', 'PMIS', 'PMISc', 'CLJP', 'CLJPc']: import split splitting = getattr(split, CF)(C) else: raise ValueError('unknown C/F splitting method (%s)' % CF) P = direct_interpolation(A, C, splitting) R = P.T.tocsr() if keep: levels[-1].C = C # strength of connection matrix levels[-1].splitting = splitting # C/F splitting levels[-1].P = P # prolongation operator levels[-1].R = R # restriction operator levels.append( multilevel_solver.level() ) A = R * A * P #galerkin operator levels[-1].A = A