This page contains the Aggregation Package documentation.
Adaptive Smoothed Aggregation
Create a multilevel solver using Adaptive Smoothed Aggregation (aSA)
Parameters : | A : {csr_matrix, bsr_matrix}
initial_candidates : {None, n x m dense matrix}
symmetry : {string}
pdef : {bool}
num_candidates : {integer}
candidate_iters : {integer}
improvement_iters : {integer}
epsilon : {float}
max_levels : {integer}
max_coarse : {integer}
prepostsmoother : {string or dict}
strength : [‘symmetric’, ‘classical’, ‘evolution’, (‘predefined’, {‘C’
aggregate : [‘standard’, ‘lloyd’, ‘naive’, (‘predefined’, {‘AggOp’
smooth : [‘jacobi’, ‘richardson’, ‘energy’, None]
coarse_solver : [‘splu’, ‘lu’, ‘cholesky, ‘pinv’, ‘gauss_seidel’, ... ]
eliminate_local : {tuple}
keep: {bool} : default False
|
---|---|
Returns : | multilevel_solver : multilevel_solver
|
Notes
References
[R5] | Brezina, Falgout, MacLachlan, Manteuffel, McCormick, and Ruge “Adaptive Smoothed Aggregation ($lpha$SA) Multigrid” SIAM Review Volume 47, Issue 2 (2005) http://www.cs.umn.edu/~maclach/research/aSA2.pdf |
Examples
>>> from pyamg.gallery import stencil_grid
>>> from pyamg.aggregation import adaptive_sa_solver
>>> import numpy
>>> A=stencil_grid([[-1,-1,-1],[-1,8.0,-1],[-1,-1,-1]], (31,31),format='csr')
>>> [asa,work] = adaptive_sa_solver(A,num_candidates=1)
>>> residuals=[]
>>> x=asa.solve(b=numpy.ones((A.shape[0],)),x0=numpy.ones((A.shape[0],)),residuals=residuals)
Aggregation methods
Compute the sparsity pattern of the tentative prolongator
Parameters : | C : csr_matrix
|
---|---|
Returns : | AggOp : csr_matrix
Cpts : array
|
See also
amg_core.standard_aggregation
Examples
>>> from scipy.sparse import csr_matrix
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation.aggregate import standard_aggregation
>>> A = poisson((4,), format='csr') # 1D mesh with 4 vertices
>>> A.todense()
matrix([[ 2., -1., 0., 0.],
[-1., 2., -1., 0.],
[ 0., -1., 2., -1.],
[ 0., 0., -1., 2.]])
>>> standard_aggregation(A)[0].todense() # two aggregates
matrix([[1, 0],
[1, 0],
[0, 1],
[0, 1]], dtype=int8)
>>> A = csr_matrix([[1,0,0],[0,1,1],[0,1,1]])
>>> A.todense() # first vertex is isolated
matrix([[1, 0, 0],
[0, 1, 1],
[0, 1, 1]])
>>> standard_aggregation(A)[0].todense() # one aggregate
matrix([[0],
[1],
[1]], dtype=int8)
Compute the sparsity pattern of the tentative prolongator
Parameters : | C : csr_matrix
|
---|---|
Returns : | AggOp : csr_matrix
Cpts : array
|
See also
amg_core.naive_aggregation
Notes
Differs from standard aggregation. Each dof is considered. If it has been aggregated, skip over. Otherwise, put dof and any unaggregated neighbors in an aggregate. Results in possibly much higher complexities than standard aggregation.
Examples
>>> from scipy.sparse import csr_matrix
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation.aggregate import naive_aggregation
>>> A = poisson((4,), format='csr') # 1D mesh with 4 vertices
>>> A.todense()
matrix([[ 2., -1., 0., 0.],
[-1., 2., -1., 0.],
[ 0., -1., 2., -1.],
[ 0., 0., -1., 2.]])
>>> naive_aggregation(A)[0].todense() # two aggregates
matrix([[1, 0],
[1, 0],
[0, 1],
[0, 1]], dtype=int8)
>>> A = csr_matrix([[1,0,0],[0,1,1],[0,1,1]])
>>> A.todense() # first vertex is isolated
matrix([[1, 0, 0],
[0, 1, 1],
[0, 1, 1]])
>>> naive_aggregation(A)[0].todense() # two aggregates
matrix([[1, 0],
[0, 1],
[0, 1]], dtype=int8)
Aggregated nodes using Lloyd Clustering
Parameters : | C : csr_matrix
ratio : scalar
distance : [‘unit’,’abs’,’inv’,None]
maxiter : int
|
||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|
Returns : | AggOp : csr_matrix
seeds : array
|
See also
amg_core.standard_aggregation
Examples
>>> from scipy.sparse import csr_matrix
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation.aggregate import lloyd_aggregation
>>> A = poisson((4,), format='csr') # 1D mesh with 4 vertices
>>> A.todense()
matrix([[ 2., -1., 0., 0.],
[-1., 2., -1., 0.],
[ 0., -1., 2., -1.],
[ 0., 0., -1., 2.]])
>>> lloyd_aggregation(A)[0].todense() # one aggregate
matrix([[1],
[1],
[1],
[1]], dtype=int8)
>>> Agg = lloyd_aggregation(A,ratio=0.5).todense() # more seeding for two aggregates
Support for aggregation-based AMG
Create a multilevel solver using classical-style Smoothed Aggregation (SA)
Parameters : | A : {csr_matrix, bsr_matrix}
B : {None, array_like}
BH : {None, array_like}
symmetry : {string}
strength : {list}
aggregate : {list}
smooth : {list}
presmoother : {tuple, string, list}
postsmoother : {tuple, string, list}
Bimprove : {list}
max_levels : {integer}
max_coarse : {integer}
keep : {bool}
|
---|---|
Returns : | ml : multilevel_solver
|
Other Parameters: | |
cycle_type : [‘V’,’W’,’F’]
coarse_solver : [‘splu’, ‘lu’, ‘cholesky, ‘pinv’, ‘gauss_seidel’, ... ]
|
See also
multilevel_solver, classical.ruge_stuben_solver, aggregation.smoothed_aggregation_solver
Notes
This method implements classical-style SA, not root-node style SA (see aggregation.rootnode_solver).
The additional parameters are passed through as arguments to multilevel_solver. Refer to pyamg.multilevel_solver for additional documentation.
At each level, four steps are executed in order to define the coarser level operator.
The parameters smooth, strength, aggregate, presmoother, postsmoother can be varied on a per level basis. For different methods on different levels, use a list as input so that the i-th entry defines the method at the i-th level. If there are more levels in the hierarchy than list entries, the last entry will define the method for all levels lower.
Examples are: smooth=[(‘jacobi’, {‘omega’:1.0}), None, ‘jacobi’] presmoother=[(‘block_gauss_seidel’, {‘sweep’:symmetric}), ‘sor’] aggregate=[‘standard’, ‘naive’] strength=[(‘symmetric’, {‘theta’:0.25}), (‘symmetric’,{‘theta’:0.08})]
Predefined strength of connection and aggregation schemes can be specified. These options are best used together, but aggregation can be predefined while strength of connection is not.
For predefined strength of connection, use a list consisting of tuples of the form (‘predefined’, {‘C’ : C0}), where C0 is a csr_matrix and each degree-of-freedom in C0 represents a supernode. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘C’ : C0}), (‘predefined’, {‘C’ : C1}) ].
Similarly for predefined aggregation, use a list of tuples. For instance to predefine a three-level hierarchy, use [(‘predefined’, {‘AggOp’ : Agg0}), (‘predefined’, {‘AggOp’ : Agg1}) ], where the dimensions of A, Agg0 and Agg1 are compatible, i.e. Agg0.shape[1] == A.shape[0] and Agg1.shape[1] == Agg0.shape[0]. Each AggOp is a csr_matrix.
References
[R6] | Vanek, P. and Mandel, J. and Brezina, M., “Algebraic Multigrid by Smoothed Aggregation for Second and Fourth Order Elliptic Problems”, Computing, vol. 56, no. 3, pp. 179–196, 1996. http://citeseer.ist.psu.edu/vanek96algebraic.html |
Examples
>>> from pyamg import smoothed_aggregation_solver
>>> from pyamg.gallery import poisson
>>> from scipy.sparse.linalg import cg
>>> import numpy
>>> A = poisson((100,100), format='csr') # matrix
>>> b = numpy.ones((A.shape[0])) # RHS
>>> ml = smoothed_aggregation_solver(A) # AMG solver
>>> M = ml.aspreconditioner(cycle='V') # preconditioner
>>> x,info = cg(A, b, tol=1e-8, maxiter=30, M=M) # solve with CG
Methods to smooth tentative prolongation operators
Jacobi prolongation smoother
Parameters : | S : {csr_matrix, bsr_matrix}
T : {csr_matrix, bsr_matrix}
C : {csr_matrix, bsr_matrix}
B : {array}
omega : {scalar}
filter : {boolean}
weighting : {string}
|
---|---|
Returns : | P : {csr_matrix, bsr_matrix}
|
Notes
If weighting is not ‘local’, then results using Jacobi prolongation smoother are not precisely reproducible due to a random initial guess used for the spectral radius approximation. For precise reproducibility, set numpy.random.seed(..) to the same value before each test.
Examples
>>> from pyamg.aggregation import jacobi_prolongation_smoother
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import coo_matrix
>>> import numpy
>>> data = numpy.ones((6,))
>>> row = numpy.arange(0,6)
>>> col = numpy.kron([0,1],numpy.ones((3,)))
>>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
>>> T.todense()
matrix([[ 1., 0.],
[ 1., 0.],
[ 1., 0.],
[ 0., 1.],
[ 0., 1.],
[ 0., 1.]])
>>> A = poisson((6,),format='csr')
>>> P = jacobi_prolongation_smoother(A,T,A,numpy.ones((2,1)))
>>> P.todense()
matrix([[ 0.64930164, 0. ],
[ 1. , 0. ],
[ 0.64930164, 0.35069836],
[ 0.35069836, 0.64930164],
[ 0. , 1. ],
[ 0. , 0.64930164]])
Richardson prolongation smoother
Parameters : | S : {csr_matrix, bsr_matrix}
T : {csr_matrix, bsr_matrix}
omega : {scalar}
|
---|---|
Returns : | P : {csr_matrix, bsr_matrix}
|
Notes
Results using Richardson prolongation smoother are not precisely reproducible due to a random initial guess used for the spectral radius approximation. For precise reproducibility, set numpy.random.seed(..) to the same value before each test.
Examples
>>> from pyamg.aggregation import richardson_prolongation_smoother
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import coo_matrix
>>> import numpy
>>> data = numpy.ones((6,))
>>> row = numpy.arange(0,6)
>>> col = numpy.kron([0,1],numpy.ones((3,)))
>>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
>>> T.todense()
matrix([[ 1., 0.],
[ 1., 0.],
[ 1., 0.],
[ 0., 1.],
[ 0., 1.],
[ 0., 1.]])
>>> A = poisson((6,),format='csr')
>>> P = richardson_prolongation_smoother(A,T)
>>> P.todense()
matrix([[ 0.64930164, 0. ],
[ 1. , 0. ],
[ 0.64930164, 0.35069836],
[ 0.35069836, 0.64930164],
[ 0. , 1. ],
[ 0. , 0.64930164]])
Minimize the energy of the coarse basis functions (columns of T). Both root-node and non-root-node style prolongation smoothing is available, see Cpt_params description below.
Parameters : | A : {csr_matrix, bsr_matrix}
T : {bsr_matrix}
Atilde : {csr_matrix}
B : {array}
Bf : {array}
Cpt_params : {tuple}
krylov : {string}
maxiter : integer
tol : {scalar}
degree : {int}
weighting : {string}
|
---|---|
Returns : | T : {bsr_matrix}
|
Notes
Only ‘diagonal’ weighting is supported for the CGNR method, because we are working with A^* A and not A.
When Cpt_params[0] == True, root-node style prolongation smoothing is used to minimize the energy of columns of T. Essentially, an identity block is maintained in T, corresponding to injection from the coarse-grid to the fine-grid root-nodes. See [2] for more details, and see util.utils.get_Cpt_params for the helper function to generate Cpt_params.
If Cpt_params[0] == False, the energy of columns of T are still minimized, but without maintaining the identity block.
References
[R7] | Jan Mandel, Marian Brezina, and Petr Vanek “Energy Optimization of Algebraic Multigrid Bases” Computing 62, 205-228, 1999 http://dx.doi.org/10.1007/s006070050022 |
[R8] | Olson, L. and Schroder, J. and Tuminaro, R., “A general interpolation strategy for algebraic multigrid using energy minimization”, SIAM Journal on Scientific Computing (SISC), vol. 33, pp. 966–991, 2011. |
Examples
>>> from pyamg.aggregation import energy_prolongation_smoother
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import coo_matrix
>>> import numpy
>>> data = numpy.ones((6,))
>>> row = numpy.arange(0,6)
>>> col = numpy.kron([0,1],numpy.ones((3,)))
>>> T = coo_matrix((data,(row,col)),shape=(6,2)).tocsr()
>>> print T.todense()
[[ 1. 0.]
[ 1. 0.]
[ 1. 0.]
[ 0. 1.]
[ 0. 1.]
[ 0. 1.]]
>>> A = poisson((6,),format='csr')
>>> P = energy_prolongation_smoother(A,T,A,numpy.ones((2,1),dtype=float), None, (False,{}) )
>>> print P.todense()
[[ 1. 0. ]
[ 1. 0. ]
[ 0.66666667 0.33333333]
[ 0.33333333 0.66666667]
[ 0. 1. ]
[ 0. 1. ]]
Tentative prolongator
Fit near-nullspace candidates to form the tentative prolongator
Parameters : | AggOp : csr_matrix
B : array
tol : scalar
|
---|---|
Returns : | (Q,R) : (bsr_matrix,array)
|
See also
amg_core.fit_candidates
Notes
Assuming that each row of AggOp contains exactly one non-zero entry, i.e. all unknowns belong to an aggregate, then Q and R satisfy the relationship B = Q*R. In other words, the near-nullspace candidates are represented exactly by the tentative prolongator.
If AggOp contains rows with no non-zero entries, then the range of the tentative prolongator will not include those degrees of freedom. This situation is illustrated in the examples below.
References
[R9] | Vanek, P. and Mandel, J. and Brezina, M., “Algebraic Multigrid by Smoothed Aggregation for Second and Fourth Order Elliptic Problems”, Computing, vol. 56, no. 3, pp. 179–196, 1996. http://citeseer.ist.psu.edu/vanek96algebraic.html |
Examples
>>> from scipy.sparse import csr_matrix
>>> from pyamg.aggregation.tentative import fit_candidates
>>> # four nodes divided into two aggregates
... AggOp = csr_matrix( [[1,0],
... [1,0],
... [0,1],
... [0,1]] )
>>> # B contains one candidate, the constant vector
... B = [[1],
... [1],
... [1],
... [1]]
>>> Q,R = fit_candidates(AggOp, B)
>>> Q.todense()
matrix([[ 0.70710678, 0. ],
[ 0.70710678, 0. ],
[ 0. , 0.70710678],
[ 0. , 0.70710678]])
>>> R
array([[ 1.41421356],
[ 1.41421356]])
>>> # Two candidates, the constant vector and a linear function
... B = [[1,0],
... [1,1],
... [1,2],
... [1,3]]
>>> Q,R = fit_candidates(AggOp, B)
>>> Q.todense()
matrix([[ 0.70710678, -0.70710678, 0. , 0. ],
[ 0.70710678, 0.70710678, 0. , 0. ],
[ 0. , 0. , 0.70710678, -0.70710678],
[ 0. , 0. , 0.70710678, 0.70710678]])
>>> R
array([[ 1.41421356, 0.70710678],
[ 0. , 0.70710678],
[ 1.41421356, 3.53553391],
[ 0. , 0.70710678]])
>>> # aggregation excludes the third node
... AggOp = csr_matrix( [[1,0],
... [1,0],
... [0,0],
... [0,1]] )
>>> B = [[1],
... [1],
... [1],
... [1]]
>>> Q,R = fit_candidates(AggOp, B)
>>> Q.todense()
matrix([[ 0.70710678, 0. ],
[ 0.70710678, 0. ],
[ 0. , 0. ],
[ 0. , 1. ]])
>>> R
array([[ 1.41421356],
[ 1. ]])