relaxation Package¶
relaxation Package¶
Relaxation methods¶
The multigrid cycle is formed by two complementary procedures: relaxation and coarse-grid correction. The role of relaxation is to rapidly damp oscillatory (high-frequency) errors out of the approximate solution. When the error is smooth, it can then be accurately represented on the coarser grid, where a solution, or approximate solution, can be computed.
Iterative methods for linear systems that have an error smoothing property are valid relaxation methods. Since the purpose of a relaxation method is to smooth oscillatory errors, its effectiveness on non-oscillatory errors is not important. This point explains why simple iterative methods like Gauss-Seidel iteration are effective relaxation methods while being very slow to converge to the solution of Ax=b.
- PyAMG implements relaxation methods of the following varieties:
- Jacobi iteration
- Gauss-Seidel iteration
- Successive Over-Relaxation
- Polynomial smoothing (e.g. Chebyshev)
- Jacobi and Gauss-Seidel on the normal equations (A.H A and A A.H)
- Krylov methods: gmres, cg, cgnr, cgne
- No pre- or postsmoother
Refer to the docstrings of the individual methods for additional information.
- pyamg.relaxation.block_gauss_seidel(A, x, b, iterations=1, sweep='forward', blocksize=1, Dinv=None)¶
Perform block Gauss-Seidel iteration on the linear system Ax=b
- A : {csr_matrix, bsr_matrix}
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
- Dinv : array
- Array holding block diagonal inverses of A size (N/blocksize, blocksize, blocksize)
- blocksize : int
- Desired dimension of blocks
Nothing, x will be modified in place.
>>> # Use Gauss-Seidel as a Stand-Alone Solver >>> from pyamg.relaxation import block_gauss_seidel >>> 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 >>> opts = {'sweep':'symmetric', 'blocksize' : 4} >>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)), ... coarse_solver='pinv2', max_coarse=50, ... presmoother=('block_gauss_seidel', opts), ... postsmoother=('block_gauss_seidel', opts)) >>> x0=numpy.zeros((A.shape[0],1)) >>> residuals=[] >>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
- pyamg.relaxation.block_jacobi(A, x, b, Dinv=None, blocksize=1, iterations=1, omega=1.0)¶
Perform block Jacobi iteration on the linear system Ax=b
- A : csr_matrix or bsr_matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- Dinv : array
- Array holding block diagonal inverses of A size (N/blocksize, blocksize, blocksize)
- blocksize : int
- Desired dimension of blocks
- iterations : int
- Number of iterations to perform
- omega : scalar
- Damping parameter
Nothing, x will be modified in place.
>>> # Use block Jacobi as a Stand-Alone Solver >>> from pyamg.relaxation import block_jacobi >>> 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 >>> opts = {'omega': 4.0/3.0, 'iterations' : 2, 'blocksize' : 4} >>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)), ... coarse_solver='pinv2', max_coarse=50, ... presmoother=('block_jacobi', opts), ... postsmoother=('block_jacobi', opts)) >>> x0=numpy.zeros((A.shape[0],1)) >>> residuals=[] >>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
- pyamg.relaxation.gauss_seidel(A, x, b, iterations=1, sweep='forward')¶
Perform Gauss-Seidel iteration on the linear system Ax=b
- A : {csr_matrix, bsr_matrix}
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
Nothing, x will be modified in place.
>>> # Use Gauss-Seidel as a Stand-Alone Solver >>> from pyamg.relaxation import gauss_seidel >>> 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)
- pyamg.relaxation.gauss_seidel_indexed(A, x, b, indices, iterations=1, sweep='forward')¶
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.
- A : csr_matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- indices : ndarray
- Row indices to relax.
- iterations : int
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
Nothing, x will be modified in place.
>>> from pyamg.gallery import poisson >>> from pyamg.relaxation import gauss_seidel_indexed >>> 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 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') # 3, then 2 >>> gauss_seidel_indexed(A, x, b, [2,0,2]) # relax row 2, 0, 2
- pyamg.relaxation.gauss_seidel_ne(A, x, b, iterations=1, sweep='forward', omega=1.0, Dinv=None)¶
- Perform Gauss-Seidel iterations on the linear system A A.H x = b
- (Also known as Kaczmarz relaxation)
- A : csr_matrix
- Sparse NxN matrix
- x : { ndarray }
- Approximate solution (length N)
- b : { ndarray }
- Right-hand side (length N)
- iterations : { int }
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
- omega : { float}
- Relaxation parameter typically in (0, 2) if omega != 1.0, then algorithm becomes SOR on A A.H
- Dinv : { ndarray}
- Inverse of diag(A A.H), (length N)
Nothing, x will be modified in place.
[1] 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; [2] Kaczmarz. Angenaeherte Aufloesung von Systemen Linearer Gleichungen. Bull. Acad. Polon. Sci. Lett. A 35, 355-57. 1937 >>> # Use NE Gauss-Seidel as a Stand-Alone Solver >>> from pyamg.relaxation import gauss_seidel_ne >>> 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)
- pyamg.relaxation.gauss_seidel_nr(A, x, b, iterations=1, sweep='forward', omega=1.0, Dinv=None)¶
Perform Gauss-Seidel iterations on the linear system A.H A x = A.H b
- A : csr_matrix
- Sparse NxN matrix
- x : { ndarray }
- Approximate solution (length N)
- b : { ndarray }
- Right-hand side (length N)
- iterations : { int }
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
- omega : { float}
- Relaxation parameter typically in (0, 2) if omega != 1.0, then algorithm becomes SOR on A.H A
- Dinv : { ndarray}
- Inverse of diag(A.H A), (length N)
Nothing, x will be modified in place.
[1] Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 247-9, 2003 http://www-users.cs.umn.edu/~saad/books.html >>> # Use NR Gauss-Seidel as a Stand-Alone Solver >>> from pyamg.relaxation import gauss_seidel_nr >>> 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)
- pyamg.relaxation.jacobi(A, x, b, iterations=1, omega=1.0)¶
Perform Jacobi iteration on the linear system Ax=b
- A : csr_matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- omega : scalar
- Damping parameter
Nothing, x will be modified in place.
>>> # Use Jacobi as a Stand-Alone Solver >>> from pyamg.relaxation import jacobi >>> 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)
- pyamg.relaxation.jacobi_ne(A, x, b, iterations=1, omega=1.0)¶
- Perform Jacobi iterations on the linear system A A.H x = A.H b
- (Also known as Cimmino relaxation)
- A : csr_matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- omega : scalar
- Damping parameter
Nothing, x will be modified in place.
[1] 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; [2] Kaczmarz. Angenaeherte Aufloesung von Systemen Linearer Gleichungen. Bull. Acad. Polon. Sci. Lett. A 35, 355-57. 1937 [3] Cimmino. La ricerca scientifica ser. II 1. Pubbliz. dell’Inst. pre le Appl. del Calculo 34, 326-333, 1938. >>> # Use NE Jacobi as a Stand-Alone Solver >>> from pyamg.relaxation import jacobi_ne >>> 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 >>> opts = {'iterations' : 2, 'omega' : 4.0/3.0} >>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)), ... coarse_solver='pinv2', max_coarse=50, ... presmoother=('jacobi_ne', opts), ... postsmoother=('jacobi_ne', opts)) >>> x0=numpy.zeros((A.shape[0],1)) >>> residuals=[] >>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
- pyamg.relaxation.polynomial(A, x, b, coefficients, iterations=1)¶
Apply a polynomial smoother to the system Ax=b
- A : sparse matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- coefficients : {array_like}
- Coefficients of the polynomial. See Notes section for details.
- iterations : int
- Number of iterations to perform
Nothing, x will be modified in place.
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’.
- Richardson iteration p(A) = c_0:
polynomial_smoother(A, x, b, [c_0])
- Linear smoother p(A) = c_1*A + c_0:
polynomial_smoother(A, x, b, [c_1, c_0])
- Quadratic smoother p(A) = c_2*A^2 + c_1*A + 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).
>>> # 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)
- pyamg.relaxation.schwarz(A, x, b, iterations=1, subdomain=None, subdomain_ptr=None, inv_subblock=None, inv_subblock_ptr=None, sweep='forward')¶
- Perform Overlapping multiplicative Schwarz on
- the linear system Ax=b
- A : {csr_matrix, bsr_matrix}
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- subdomain : {int array}
- Linear array containing each subdomain’s elements
- subdomain_ptr : {int array}
- Pointer in subdomain, such that subdomain[subdomain_ptr[i]:subdomain_ptr[i+1]]] contains the _sorted_ indices in subdomain i
- inv_subblock : {int_array}
- Linear array containing each subdomain’s inverted diagonal block of A
- inv_subblock_ptr : {int array}
- Pointer in inv_subblock, such that inv_subblock[inv_subblock_ptr[i]:inv_subblock_ptr[i+1]]] contains the inverted diagonal block of A for the i-th subdomain in _row_ major order
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
Nothing, x will be modified in place.
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
>>> # Use Overlapping Schwarz as a Stand-Alone Solver >>> from pyamg.relaxation import schwarz >>> 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)
- pyamg.relaxation.sor(A, x, b, omega, iterations=1, sweep='forward')¶
Perform SOR iteration on the linear system Ax=b
- A : {csr_matrix, bsr_matrix}
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- omega : scalar
- Damping parameter
- iterations : int
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
Nothing, x will be modified in place.
When omega=1.0, SOR is equivalent to Gauss-Seidel.
>>> # 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)
chebyshev Module¶
Compute coefficients for polynomial smoothers
- pyamg.relaxation.chebyshev.chebyshev_polynomial_coefficients(a, b, degree)[source]¶
Chebyshev polynomial coefficients for the interval [a,b]
- a,b : float
- The left and right endpoints of the interval.
- degree : int
- Degree of desired Chebyshev polynomial
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.
a,b typically represent the interval of the spectrum for some matrix that you wish to damp with a Chebyshev smoother.
>>> from pyamg.relaxation.chebyshev import chebyshev_polynomial_coefficients >>> print chebyshev_polynomial_coefficients(1.0,2.0, 3) [-0.32323232 1.45454545 -2.12121212 1. ]
info Module¶
Relaxation methods¶
The multigrid cycle is formed by two complementary procedures: relaxation and coarse-grid correction. The role of relaxation is to rapidly damp oscillatory (high-frequency) errors out of the approximate solution. When the error is smooth, it can then be accurately represented on the coarser grid, where a solution, or approximate solution, can be computed.
Iterative methods for linear systems that have an error smoothing property are valid relaxation methods. Since the purpose of a relaxation method is to smooth oscillatory errors, its effectiveness on non-oscillatory errors is not important. This point explains why simple iterative methods like Gauss-Seidel iteration are effective relaxation methods while being very slow to converge to the solution of Ax=b.
- PyAMG implements relaxation methods of the following varieties:
- Jacobi iteration
- Gauss-Seidel iteration
- Successive Over-Relaxation
- Polynomial smoothing (e.g. Chebyshev)
- Jacobi and Gauss-Seidel on the normal equations (A.H A and A A.H)
- Krylov methods: gmres, cg, cgnr, cgne
- No pre- or postsmoother
Refer to the docstrings of the individual methods for additional information.
relaxation Module¶
Relaxation methods for linear systems
- pyamg.relaxation.relaxation.sor(A, x, b, omega, iterations=1, sweep='forward')[source]¶
Perform SOR iteration on the linear system Ax=b
- A : {csr_matrix, bsr_matrix}
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- omega : scalar
- Damping parameter
- iterations : int
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
Nothing, x will be modified in place.
When omega=1.0, SOR is equivalent to Gauss-Seidel.
>>> # 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)
- pyamg.relaxation.relaxation.gauss_seidel(A, x, b, iterations=1, sweep='forward')[source]¶
Perform Gauss-Seidel iteration on the linear system Ax=b
- A : {csr_matrix, bsr_matrix}
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
Nothing, x will be modified in place.
>>> # Use Gauss-Seidel as a Stand-Alone Solver >>> from pyamg.relaxation import gauss_seidel >>> 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)
- pyamg.relaxation.relaxation.jacobi(A, x, b, iterations=1, omega=1.0)[source]¶
Perform Jacobi iteration on the linear system Ax=b
- A : csr_matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- omega : scalar
- Damping parameter
Nothing, x will be modified in place.
>>> # Use Jacobi as a Stand-Alone Solver >>> from pyamg.relaxation import jacobi >>> 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)
- pyamg.relaxation.relaxation.polynomial(A, x, b, coefficients, iterations=1)[source]¶
Apply a polynomial smoother to the system Ax=b
- A : sparse matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- coefficients : {array_like}
- Coefficients of the polynomial. See Notes section for details.
- iterations : int
- Number of iterations to perform
Nothing, x will be modified in place.
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’.
- Richardson iteration p(A) = c_0:
polynomial_smoother(A, x, b, [c_0])
- Linear smoother p(A) = c_1*A + c_0:
polynomial_smoother(A, x, b, [c_1, c_0])
- Quadratic smoother p(A) = c_2*A^2 + c_1*A + 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).
>>> # 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)
- pyamg.relaxation.relaxation.schwarz(A, x, b, iterations=1, subdomain=None, subdomain_ptr=None, inv_subblock=None, inv_subblock_ptr=None, sweep='forward')[source]¶
- Perform Overlapping multiplicative Schwarz on
- the linear system Ax=b
- A : {csr_matrix, bsr_matrix}
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- subdomain : {int array}
- Linear array containing each subdomain’s elements
- subdomain_ptr : {int array}
- Pointer in subdomain, such that subdomain[subdomain_ptr[i]:subdomain_ptr[i+1]]] contains the _sorted_ indices in subdomain i
- inv_subblock : {int_array}
- Linear array containing each subdomain’s inverted diagonal block of A
- inv_subblock_ptr : {int array}
- Pointer in inv_subblock, such that inv_subblock[inv_subblock_ptr[i]:inv_subblock_ptr[i+1]]] contains the inverted diagonal block of A for the i-th subdomain in _row_ major order
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
Nothing, x will be modified in place.
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
>>> # Use Overlapping Schwarz as a Stand-Alone Solver >>> from pyamg.relaxation import schwarz >>> 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)
- pyamg.relaxation.relaxation.jacobi_ne(A, x, b, iterations=1, omega=1.0)[source]¶
- Perform Jacobi iterations on the linear system A A.H x = A.H b
- (Also known as Cimmino relaxation)
- A : csr_matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- omega : scalar
- Damping parameter
Nothing, x will be modified in place.
[1] 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; [2] Kaczmarz. Angenaeherte Aufloesung von Systemen Linearer Gleichungen. Bull. Acad. Polon. Sci. Lett. A 35, 355-57. 1937 [3] Cimmino. La ricerca scientifica ser. II 1. Pubbliz. dell’Inst. pre le Appl. del Calculo 34, 326-333, 1938. >>> # Use NE Jacobi as a Stand-Alone Solver >>> from pyamg.relaxation import jacobi_ne >>> 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 >>> opts = {'iterations' : 2, 'omega' : 4.0/3.0} >>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)), ... coarse_solver='pinv2', max_coarse=50, ... presmoother=('jacobi_ne', opts), ... postsmoother=('jacobi_ne', opts)) >>> x0=numpy.zeros((A.shape[0],1)) >>> residuals=[] >>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
- pyamg.relaxation.relaxation.gauss_seidel_ne(A, x, b, iterations=1, sweep='forward', omega=1.0, Dinv=None)[source]¶
- Perform Gauss-Seidel iterations on the linear system A A.H x = b
- (Also known as Kaczmarz relaxation)
- A : csr_matrix
- Sparse NxN matrix
- x : { ndarray }
- Approximate solution (length N)
- b : { ndarray }
- Right-hand side (length N)
- iterations : { int }
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
- omega : { float}
- Relaxation parameter typically in (0, 2) if omega != 1.0, then algorithm becomes SOR on A A.H
- Dinv : { ndarray}
- Inverse of diag(A A.H), (length N)
Nothing, x will be modified in place.
[1] 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; [2] Kaczmarz. Angenaeherte Aufloesung von Systemen Linearer Gleichungen. Bull. Acad. Polon. Sci. Lett. A 35, 355-57. 1937 >>> # Use NE Gauss-Seidel as a Stand-Alone Solver >>> from pyamg.relaxation import gauss_seidel_ne >>> 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)
- pyamg.relaxation.relaxation.gauss_seidel_nr(A, x, b, iterations=1, sweep='forward', omega=1.0, Dinv=None)[source]¶
Perform Gauss-Seidel iterations on the linear system A.H A x = A.H b
- A : csr_matrix
- Sparse NxN matrix
- x : { ndarray }
- Approximate solution (length N)
- b : { ndarray }
- Right-hand side (length N)
- iterations : { int }
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
- omega : { float}
- Relaxation parameter typically in (0, 2) if omega != 1.0, then algorithm becomes SOR on A.H A
- Dinv : { ndarray}
- Inverse of diag(A.H A), (length N)
Nothing, x will be modified in place.
[1] Yousef Saad, “Iterative Methods for Sparse Linear Systems, Second Edition”, SIAM, pp. 247-9, 2003 http://www-users.cs.umn.edu/~saad/books.html >>> # Use NR Gauss-Seidel as a Stand-Alone Solver >>> from pyamg.relaxation import gauss_seidel_nr >>> 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)
- pyamg.relaxation.relaxation.gauss_seidel_indexed(A, x, b, indices, iterations=1, sweep='forward')[source]¶
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.
- A : csr_matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- indices : ndarray
- Row indices to relax.
- iterations : int
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
Nothing, x will be modified in place.
>>> from pyamg.gallery import poisson >>> from pyamg.relaxation import gauss_seidel_indexed >>> 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 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') # 3, then 2 >>> gauss_seidel_indexed(A, x, b, [2,0,2]) # relax row 2, 0, 2
- pyamg.relaxation.relaxation.block_jacobi(A, x, b, Dinv=None, blocksize=1, iterations=1, omega=1.0)[source]¶
Perform block Jacobi iteration on the linear system Ax=b
- A : csr_matrix or bsr_matrix
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- Dinv : array
- Array holding block diagonal inverses of A size (N/blocksize, blocksize, blocksize)
- blocksize : int
- Desired dimension of blocks
- iterations : int
- Number of iterations to perform
- omega : scalar
- Damping parameter
Nothing, x will be modified in place.
>>> # Use block Jacobi as a Stand-Alone Solver >>> from pyamg.relaxation import block_jacobi >>> 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 >>> opts = {'omega': 4.0/3.0, 'iterations' : 2, 'blocksize' : 4} >>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)), ... coarse_solver='pinv2', max_coarse=50, ... presmoother=('block_jacobi', opts), ... postsmoother=('block_jacobi', opts)) >>> x0=numpy.zeros((A.shape[0],1)) >>> residuals=[] >>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
- pyamg.relaxation.relaxation.block_gauss_seidel(A, x, b, iterations=1, sweep='forward', blocksize=1, Dinv=None)[source]¶
Perform block Gauss-Seidel iteration on the linear system Ax=b
- A : {csr_matrix, bsr_matrix}
- Sparse NxN matrix
- x : ndarray
- Approximate solution (length N)
- b : ndarray
- Right-hand side (length N)
- iterations : int
- Number of iterations to perform
- sweep : {‘forward’,’backward’,’symmetric’}
- Direction of sweep
- Dinv : array
- Array holding block diagonal inverses of A size (N/blocksize, blocksize, blocksize)
- blocksize : int
- Desired dimension of blocks
Nothing, x will be modified in place.
>>> # Use Gauss-Seidel as a Stand-Alone Solver >>> from pyamg.relaxation import block_gauss_seidel >>> 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 >>> opts = {'sweep':'symmetric', 'blocksize' : 4} >>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)), ... coarse_solver='pinv2', max_coarse=50, ... presmoother=('block_gauss_seidel', opts), ... postsmoother=('block_gauss_seidel', opts)) >>> x0=numpy.zeros((A.shape[0],1)) >>> residuals=[] >>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
smoothing Module¶
Method to create pre and post-smoothers on the levels of a multilevel_solver
- pyamg.relaxation.smoothing.change_smoothers(ml, presmoother, postsmoother)[source]¶
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.
- ml : {pyamg multilevel hierarchy}
- Data structure that stores the multigrid hierarchy.
- presmoother : {None, string, tuple, list}
presmoother can be (1) the name of a supported smoother, e.g. “gauss_seidel”, (2) a tuple of the form (‘method’,’opts’) where ‘method’ is the name of a supported smoother and ‘opts’ a dict of keyword arguments to the smoother, or (3) a list of instances of options 1 or 2. See the Examples section for illustrations of the format.
If presmoother is a list, presmoother[i] determines the smoothing strategy for level i. Else, presmoother defines the same strategy for all levels.
If len(presmoother) < len(ml.levels), then presmoother[-1] is used for all remaining levels
If len(presmoother) > len(ml.levels), then the remaining smoothing strategies are ignored
- postsmoother : {string, tuple, list}
- Defines postsmoother in identical fashion to presmoother
ml changed in place ml.level[i].presmoother <=== presmoother[i] ml.level[i].postsmoother <=== postsmoother[i]
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
>>> 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)