This page contains the Relaxation Package documentation.
Compute coefficients for polynomial smoothers
Chebyshev polynomial coefficients for the interval [a,b]
Parameters : | a,b : float
degree : int
|
---|---|
Returns : | Coefficients of the Chebyshev polynomial C(t) with minimum : magnitude on the interval [a,b] such that C(0) = 1.0. : The coefficients are returned in descending order. : |
Notes
a,b typically represent the interval of the spectrum for some matrix that you wish to damp with a Chebyshev smoother.
Examples
>>> from pyamg.relaxation.chebyshev import chebyshev_polynomial_coefficients
>>> print chebyshev_polynomial_coefficients(1.0,2.0, 3)
[-0.32323232 1.45454545 -2.12121212 1. ]
Relaxation methods for linear systems
Perform SOR iteration on the linear system Ax=b
Parameters : | A : {csr_matrix, bsr_matrix}
x : ndarray
b : ndarray
omega : scalar
iterations : int
sweep : {‘forward’,’backward’,’symmetric’}
|
---|---|
Returns : | Nothing, x will be modified in place. : |
Notes
When omega=1.0, SOR is equivalent to Gauss-Seidel.
Examples
>>> ## Use SOR as stand-along solver
>>> from pyamg.relaxation import sor
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> sor(A, x0, b, 1.33, iterations=10)
>>> print norm(b-A*x0)
3.03888724811
>>> #
>>> ## Use SOR as the multigrid smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('sor', {'sweep':'symmetric', 'omega' : 1.33}),
... postsmoother=('sor', {'sweep':'symmetric', 'omega' : 1.33}))
>>> x0 = numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Perform Gauss-Seidel iteration on the linear system Ax=b
Parameters : | A : {csr_matrix, bsr_matrix}
x : ndarray
b : ndarray
iterations : int
sweep : {‘forward’,’backward’,’symmetric’}
|
---|---|
Returns : | Nothing, x will be modified in place. : |
Examples
>>> ## Use Gauss-Seidel as a Stand-Alone Solver
>>> from pyamg.relaxation import *
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> gauss_seidel(A, x0, b, iterations=10)
>>> print norm(b-A*x0)
4.00733716236
>>> #
>>> ## Use Gauss-Seidel as the Multigrid Smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('gauss_seidel', {'sweep':'symmetric'}),
... postsmoother=('gauss_seidel', {'sweep':'symmetric'}))
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Perform Jacobi iteration on the linear system Ax=b
Parameters : | A : csr_matrix
x : ndarray
b : ndarray
iterations : int
omega : scalar
|
---|---|
Returns : | Nothing, x will be modified in place. : |
Examples
>>> ## Use Jacobi as a Stand-Alone Solver
>>> from pyamg.relaxation import *
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> jacobi(A, x0, b, iterations=10, omega=1.0)
>>> print norm(b-A*x0)
5.83475132751
>>> #
>>> ## Use Jacobi as the Multigrid Smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('jacobi', {'omega': 4.0/3.0, 'iterations' : 2}),
... postsmoother=('jacobi', {'omega': 4.0/3.0, 'iterations' : 2}))
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Apply a polynomial smoother to the system Ax=b
Parameters : | A : sparse matrix
x : ndarray
b : ndarray
coefficients : {array_like}
iterations : int
|
---|---|
Returns : | Nothing, x will be modified in place. : |
Notes
The smoother has the form x[:] = x + p(A) (b - A*x) where p(A) is a polynomial in A whose scalar coefficients are specified (in descending order) by argument ‘coefficients’.
polynomial_smoother(A, x, b, [c_0])
polynomial_smoother(A, x, b, [c_1, c_0])
polynomial_smoother(A, x, b, [c_2, c_1, c_0])
Here, Horner’s Rule is applied to avoid computing A^k directly.
For efficience, the method detects the case x = 0 one matrix-vector product is avoided (since (b - A*x) is b).
Examples
>>> ## The polynomial smoother is not currently used directly
>>> ## in PyAMG. It is only used by the chebyshev smoothing option,
>>> ## which automatically calculates the correct coefficients.
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> from pyamg.aggregation import smoothed_aggregation_solver
>>> A = poisson((10,10), format='csr')
>>> b = numpy.ones((A.shape[0],1))
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('chebyshev', {'degree':3, 'iterations':1}),
... postsmoother=('chebyshev', {'degree':3, 'iterations':1}))
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Perform Overlapping multiplicative Schwarz on the linear system Ax=b
Parameters : | A : {csr_matrix, bsr_matrix}
x : ndarray
b : ndarray
iterations : int
subdomain : {int array}
subdomain_ptr : {int array}
inv_subblock : {int_array}
inv_subblock_ptr : {int array}
sweep : {‘forward’,’backward’,’symmetric’}
|
---|---|
Returns : | Nothing, x will be modified in place. : |
Notes
If subdomains is None, then a point-wise iteration takes place, with the overlapping region defined by each degree-of-freedom’s neighbors in the matrix graph.
If subdomains is not None, but subblocks is, then the subblocks are formed internally.
Currently only supports CSR matrices
Examples
>>> ## Use Overlapping Schwarz as a Stand-Alone Solver
>>> from pyamg.relaxation import *
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> schwarz(A, x0, b, iterations=10)
>>> print norm(b-A*x0)
0.126326160522
>>> #
>>> ## Schwarz as the Multigrid Smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother='schwarz',
... postsmoother='schwarz')
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Perform Jacobi iterations on the linear system A A.H x = A.H b (Also known as Cimmino relaxation)
Parameters : | A : csr_matrix
x : ndarray
b : ndarray
iterations : int
omega : scalar
|
---|---|
Returns : | Nothing, x will be modified in place. : |
References
[R28] | Brandt, Ta’asan. “Multigrid Method For Nearly Singular And Slightly Indefinite Problems.” 1985. NASA Technical Report Numbers: ICASE-85-57; NAS 1.26:178026; NASA-CR-178026; |
[R29] | Kaczmarz. Angenaeherte Aufloesung von Systemen Linearer Gleichungen. Bull. Acad. Polon. Sci. Lett. A 35, 355-57. 1937 |
[R30] | Cimmino. La ricerca scientifica ser. II 1. Pubbliz. dell’Inst. pre le Appl. del Calculo 34, 326-333, 1938. |
Examples
>>> ## Use NE Jacobi as a Stand-Alone Solver
>>> from pyamg.relaxation import *
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((50,50), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> jacobi_ne(A, x0, b, iterations=10, omega=2.0/3.0)
>>> print norm(b-A*x0)
49.3886046066
>>> #
>>> ## Use NE Jacobi as the Multigrid Smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('jacobi_ne', {'iterations' : 2, 'omega' : 4.0/3.0}),
... postsmoother=('jacobi_ne', {'iterations' : 2, 'omega' : 4.0/3.0}))
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Perform Gauss-Seidel iterations on the linear system A A.H x = b (Also known as Kaczmarz relaxation)
Parameters : | A : csr_matrix
x : { ndarray }
b : { ndarray }
iterations : { int }
sweep : {‘forward’,’backward’,’symmetric’}
omega : { float}
Dinv : { ndarray}
|
---|---|
Returns : | Nothing, x will be modified in place. : |
References
[R31] | Brandt, Ta’asan. “Multigrid Method For Nearly Singular And Slightly Indefinite Problems.” 1985. NASA Technical Report Numbers: ICASE-85-57; NAS 1.26:178026; NASA-CR-178026; |
[R32] | Kaczmarz. Angenaeherte Aufloesung von Systemen Linearer Gleichungen. Bull. Acad. Polon. Sci. Lett. A 35, 355-57. 1937 |
Examples
>>> ## Use NE Gauss-Seidel as a Stand-Alone Solver
>>> from pyamg.relaxation import *
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> gauss_seidel_ne(A, x0, b, iterations=10, sweep='symmetric')
>>> print norm(b-A*x0)
8.47576806771
>>> #
>>> ## Use NE Gauss-Seidel as the Multigrid Smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('gauss_seidel_ne', {'sweep' : 'symmetric'}),
... postsmoother=('gauss_seidel_ne', {'sweep' : 'symmetric'}))
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Perform Gauss-Seidel iterations on the linear system A.H A x = A.H b
Parameters : | A : csr_matrix
x : { ndarray }
b : { ndarray }
iterations : { int }
sweep : {‘forward’,’backward’,’symmetric’}
omega : { float}
Dinv : { ndarray}
|
---|---|
Returns : | Nothing, x will be modified in place. : |
References
[R33] | Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 247-9, 2003 http://www-users.cs.umn.edu/~saad/books.html |
Examples
>>> ## Use NR Gauss-Seidel as a Stand-Alone Solver
>>> from pyamg.relaxation import *
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> gauss_seidel_nr(A, x0, b, iterations=10, sweep='symmetric')
>>> print norm(b-A*x0)
8.45044864352
>>> #
>>> ## Use NR Gauss-Seidel as the Multigrid Smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('gauss_seidel_nr', {'sweep' : 'symmetric'}),
... postsmoother=('gauss_seidel_nr', {'sweep' : 'symmetric'}))
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Perform indexed Gauss-Seidel iteration on the linear system Ax=b
In indexed Gauss-Seidel, the sequence in which unknowns are relaxed is specified explicitly. In contrast, the standard Gauss-Seidel method always performs complete sweeps of all variables in increasing or decreasing order. The indexed method may be used to implement specialized smoothers, like F-smoothing in Classical AMG.
Parameters : | A : csr_matrix
x : ndarray
b : ndarray
indices : ndarray
iterations : int
sweep : {‘forward’,’backward’,’symmetric’}
|
---|---|
Returns : | Nothing, x will be modified in place. : |
Examples
>>> from pyamg.gallery import poisson
>>> import numpy
>>> A = poisson((4,), format='csr')
>>> x = numpy.array([0.0, 0.0, 0.0, 0.0])
>>> b = numpy.array([0.0, 1.0, 2.0, 3.0])
>>> gauss_seidel_indexed(A, x, b, [0,1,2,3]) #relax all four rows, in order
>>> gauss_seidel_indexed(A, x, b, [0,1]) #relax first two rows
>>> gauss_seidel_indexed(A, x, b, [2,0]) #relax row 2, then row 0
>>> gauss_seidel_indexed(A, x, b, [2,3], sweep='backward') #relax row 3, then row 2
>>> gauss_seidel_indexed(A, x, b, [2,0,2]) #relax row 2, then 0, then 2 again
Perform block Jacobi iteration on the linear system Ax=b
Parameters : | A : csr_matrix or bsr_matrix
x : ndarray
b : ndarray
Dinv : array
blocksize : int
iterations : int
omega : scalar
|
---|---|
Returns : | Nothing, x will be modified in place. : |
Examples
>>> ## Use block Jacobi as a Stand-Alone Solver
>>> from pyamg.relaxation import *
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> block_jacobi(A, x0, b, blocksize=4, iterations=10, omega=1.0)
>>> print norm(b-A*x0)
4.66474230129
>>> #
>>> ## Use block Jacobi as the Multigrid Smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('block_jacobi', {'omega': 4.0/3.0, 'iterations' : 2, 'blocksize' : 4}),
... postsmoother=('block_jacobi', {'omega': 4.0/3.0, 'iterations' : 2, 'blocksize' : 4}))
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Perform block Gauss-Seidel iteration on the linear system Ax=b
Parameters : | A : {csr_matrix, bsr_matrix}
x : ndarray
b : ndarray
iterations : int
sweep : {‘forward’,’backward’,’symmetric’}
Dinv : array
blocksize : int
|
---|---|
Returns : | Nothing, x will be modified in place. : |
Examples
>>> ## Use Gauss-Seidel as a Stand-Alone Solver
>>> from pyamg.relaxation import *
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> block_gauss_seidel(A, x0, b, iterations=10, blocksize=4, sweep='symmetric')
>>> print norm(b-A*x0)
0.958333817624
>>> #
>>> ## Use Gauss-Seidel as the Multigrid Smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('block_gauss_seidel', {'sweep':'symmetric', 'blocksize' : 4}),
... postsmoother=('block_gauss_seidel', {'sweep':'symmetric', 'blocksize' : 4}))
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Method to create pre- and post-smoothers on the levels of a multilevel_solver
Initialize pre- and post- smoothers throughout a multilevel_solver, with the option of having different smoothers at different levels
For each level of the multilevel_solver ‘ml’ (except the coarsest level), initialize the .presmoother() and .postsmoother() methods used in the multigrid cycle.
Parameters : | ml : {pyamg multilevel hierarchy}
presmoother : {None, string, tuple, list}
postsmoother : {string, tuple, list}
|
---|---|
Returns : | ml changed in place : ml.level[i].presmoother <=== presmoother[i] : ml.level[i].postsmoother <=== postsmoother[i] : |
Notes
Parameter ‘omega’ of the Jacobi, Richardson, and jacobi_ne methods is scaled by the spectral radius of the matrix on each level. Therefore ‘omega’ should be in the interval (0,2).
Parameter ‘withrho’ (default: True) controls whether the omega is rescaled by the spectral radius in jacobi, block_jacobi, and jacobi_ne
By initializing the smoothers after the hierarchy has been setup, allows for “algebraically” directed relaxation, such as strength_based_schwarz, which uses only the strong connections of a degree-of-freedom to define overlapping regions
Available smoother methods:
gauss_seidel
block_gauss_seidel
jacobi
block_jacobi
richardson
sor
chebyshev
gauss_seidel_nr
gauss_seidel_ne
jacobi_ne
cg
gmres
cgne
cgnr
schwarz
strength_based_schwarz
None
Examples
>>> from pyamg.gallery import poisson
>>> from pyamg.aggregation import smoothed_aggregation_solver
>>> from pyamg.relaxation.smoothing import change_smoothers
>>> from pyamg.util.linalg import norm
>>> from scipy import rand, array, mean
>>> A = poisson((10,10), format='csr')
>>> b = rand(A.shape[0],)
>>> ml = smoothed_aggregation_solver(A, max_coarse=10)
>>> #
>>> ## Set all levels to use gauss_seidel's defaults
>>> smoothers = 'gauss_seidel'
>>> change_smoothers(ml, presmoother=smoothers, postsmoother=smoothers)
>>> residuals=[]
>>> x = ml.solve(b, tol=1e-8, residuals=residuals)
>>> #
>>> ## Set all levels to use three iterations of gauss_seidel's defaults
>>> smoothers = ('gauss_seidel', {'iterations' : 3})
>>> change_smoothers(ml, presmoother=smoothers, postsmoother=None)
>>> residuals=[]
>>> x = ml.solve(b, tol=1e-8, residuals=residuals)
>>> #
>>> ## Set level 0 to use gauss_seidel's defaults, and all
>>> ## subsequent levels to use 5 iterations of cgnr
>>> smoothers = ['gauss_seidel', ('cgnr', {'maxiter' : 5})]
>>> change_smoothers(ml, presmoother=smoothers, postsmoother=smoothers)
>>> residuals=[]
>>> x = ml.solve(b, tol=1e-8, residuals=residuals)