"""Relaxation methods for linear systems"""
__docformat__ = "restructuredtext en"
from warnings import warn
import numpy
from scipy import sparse
from pyamg.util.utils import type_prep, get_diagonal, get_block_diag
from pyamg.relaxation.smoothing import schwarz_parameters
from pyamg import amg_core
__all__ = ['sor', 'gauss_seidel', 'jacobi', 'polynomial', 'schwarz']
__all__ += ['jacobi_ne', 'gauss_seidel_ne', 'gauss_seidel_nr']
__all__ += ['gauss_seidel_indexed', 'block_jacobi', 'block_gauss_seidel']
def make_system(A, x, b, formats=None):
"""
Return A,x,b suitable for relaxation or raise an exception
Parameters
----------
A : {sparse-matrix}
n x n system
x : {array}
n-vector, initial guess
b : {array}
n-vector, right-hand side
formats: {'csr', 'csc', 'bsr', 'lil', 'dok',...}
desired sparse matrix format
default is no change to A's format
Returns
-------
(A,x,b), where A is in the desired sparse-matrix format
and x and b are "raveled", i.e. (n,) vectors.
Notes
-----
Does some rudimentary error checking on the system,
such as checking for compatible dimensions and checking
for compatible type, i.e. float or complex.
Examples
--------
>>> from pyamg.relaxation.relaxation import make_system
>>> from pyamg.gallery import poisson
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> (A,x,b) = make_system(A,x,b,formats=['csc'])
>>> print str(x.shape)
(100,)
>>> print str(b.shape)
(100,)
>>> print A.format
csc
"""
if formats is None:
pass
elif formats == ['csr']:
if sparse.isspmatrix_csr(A):
pass
elif sparse.isspmatrix_bsr(A):
A = A.tocsr()
else:
warn('implicit conversion to CSR', sparse.SparseEfficiencyWarning)
A = sparse.csr_matrix(A)
else:
if sparse.isspmatrix(A) and A.format in formats:
pass
else:
A = sparse.csr_matrix(A).asformat(formats[0])
if not isinstance(x, numpy.ndarray):
raise ValueError('expected numpy array for argument x')
if not isinstance(b, numpy.ndarray):
raise ValueError('expected numpy array for argument b')
M, N = A.shape
if M != N:
raise ValueError('expected square matrix')
if x.shape not in [(M,), (M, 1)]:
raise ValueError('x has invalid dimensions')
if b.shape not in [(M,), (M, 1)]:
raise ValueError('b has invalid dimensions')
if A.dtype != x.dtype or A.dtype != b.dtype:
raise TypeError('arguments A, x, and b must have the same dtype')
if not x.flags.carray:
raise ValueError('x must be contiguous in memory')
x = numpy.ravel(x)
b = numpy.ravel(b)
return A, x, b
[docs]def sor(A, x, b, omega, iterations=1, sweep='forward'):
"""Perform SOR iteration on the linear system Ax=b
Parameters
----------
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
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)
"""
A, x, b = make_system(A, x, b, formats=['csr', 'bsr'])
x_old = numpy.empty_like(x)
for i in range(iterations):
x_old[:] = x
gauss_seidel(A, x, b, iterations=1, sweep=sweep)
x *= omega
x_old *= (1-omega)
x += x_old
[docs]def 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
Parameters
----------
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
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 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)
"""
A, x, b = make_system(A, x, b, formats=['csr'])
if subdomain is None and inv_subblock is not None:
raise ValueError("inv_subblock must be None if subdomain is None")
# If no subdomains are defined, default is to use the sparsity pattern of A
# to define the overlapping regions
(subdomain, subdomain_ptr, inv_subblock, inv_subblock_ptr) = \
schwarz_parameters(A, subdomain, subdomain_ptr,
inv_subblock, inv_subblock_ptr)
if sweep == 'forward':
row_start, row_stop, row_step = 0, subdomain_ptr.shape[0]-1, 1
elif sweep == 'backward':
row_start, row_stop, row_step = subdomain_ptr.shape[0]-2, -1, -1
elif sweep == 'symmetric':
for iter in xrange(iterations):
schwarz(A, x, b, iterations=1, subdomain=subdomain,
subdomain_ptr=subdomain_ptr, inv_subblock=inv_subblock,
inv_subblock_ptr=inv_subblock_ptr, sweep='forward')
schwarz(A, x, b, iterations=1, subdomain=subdomain,
subdomain_ptr=subdomain_ptr, inv_subblock=inv_subblock,
inv_subblock_ptr=inv_subblock_ptr, sweep='backward')
return
else:
raise ValueError("valid sweep directions are 'forward',\
'backward', and 'symmetric'")
# Call C code, need to make sure that subdomains are sorted and unique
for iter in xrange(iterations):
amg_core.overlapping_schwarz_csr(A.indptr, A.indices, A.data,
x, b, inv_subblock, inv_subblock_ptr,
subdomain, subdomain_ptr,
subdomain_ptr.shape[0]-1, A.shape[0],
row_start, row_stop, row_step)
[docs]def gauss_seidel(A, x, b, iterations=1, sweep='forward'):
"""Perform Gauss-Seidel iteration on the linear system Ax=b
Parameters
----------
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
Returns
-------
Nothing, x will be modified in place.
Examples
--------
>>> # 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)
"""
A, x, b = make_system(A, x, b, formats=['csr', 'bsr'])
if sweep == 'forward':
row_start, row_stop, row_step = 0, len(x), 1
elif sweep == 'backward':
row_start, row_stop, row_step = len(x)-1, -1, -1
elif sweep == 'symmetric':
for iter in xrange(iterations):
gauss_seidel(A, x, b, iterations=1, sweep='forward')
gauss_seidel(A, x, b, iterations=1, sweep='backward')
return
else:
raise ValueError("valid sweep directions are 'forward',\
'backward', and 'symmetric'")
if sparse.isspmatrix_csr(A):
for iter in xrange(iterations):
amg_core.gauss_seidel(A.indptr, A.indices, A.data, x, b,
row_start, row_stop, row_step)
else:
R, C = A.blocksize
if R != C:
raise ValueError('BSR blocks must be square')
row_start = row_start / R
row_stop = row_stop / R
for iter in xrange(iterations):
amg_core.bsr_gauss_seidel(A.indptr, A.indices, numpy.ravel(A.data),
x, b, row_start, row_stop, row_step, R)
[docs]def jacobi(A, x, b, iterations=1, omega=1.0):
"""Perform Jacobi iteration on the linear system Ax=b
Parameters
----------
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
Returns
-------
Nothing, x will be modified in place.
Examples
--------
>>> # 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)
"""
A, x, b = make_system(A, x, b, formats=['csr', 'bsr'])
sweep = slice(None)
(row_start, row_stop, row_step) = sweep.indices(A.shape[0])
if (row_stop - row_start) * row_step <= 0: # no work to do
return
temp = numpy.empty_like(x)
# Create uniform type, convert possibly complex scalars to length 1 arrays
[omega] = type_prep(A.dtype, [omega])
if sparse.isspmatrix_csr(A):
for iter in xrange(iterations):
amg_core.jacobi(A.indptr, A.indices, A.data, x, b, temp,
row_start, row_stop, row_step, omega)
else:
R, C = A.blocksize
if R != C:
raise ValueError('BSR blocks must be square')
row_start = row_start / R
row_stop = row_stop / R
for iter in xrange(iterations):
amg_core.bsr_jacobi(A.indptr, A.indices, numpy.ravel(A.data),
x, b, temp, row_start, row_stop,
row_step, R, omega)
[docs]def block_jacobi(A, x, b, Dinv=None, blocksize=1, iterations=1, omega=1.0):
"""Perform block Jacobi iteration on the linear system Ax=b
Parameters
----------
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
Returns
-------
Nothing, x will be modified in place.
Examples
--------
>>> # 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)
"""
A, x, b = make_system(A, x, b, formats=['csr', 'bsr'])
A = A.tobsr(blocksize=(blocksize, blocksize))
if Dinv is None:
Dinv = get_block_diag(A, blocksize=blocksize, inv_flag=True)
elif Dinv.shape[0] != A.shape[0]/blocksize:
raise ValueError('Dinv and A have incompatible dimensions')
elif (Dinv.shape[1] != blocksize) or (Dinv.shape[2] != blocksize):
raise ValueError('Dinv and blocksize are incompatible')
sweep = slice(None)
(row_start, row_stop, row_step) = sweep.indices(A.shape[0]/blocksize)
if (row_stop - row_start) * row_step <= 0: # no work to do
return
temp = numpy.empty_like(x)
# Create uniform type, convert possibly complex scalars to length 1 arrays
[omega] = type_prep(A.dtype, [omega])
for iter in xrange(iterations):
amg_core.block_jacobi(A.indptr, A.indices, numpy.ravel(A.data),
x, b, numpy.ravel(Dinv), temp,
row_start, row_stop, row_step,
omega, blocksize)
[docs]def 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
Parameters
----------
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
Returns
-------
Nothing, x will be modified in place.
Examples
--------
>>> # 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)
"""
A, x, b = make_system(A, x, b, formats=['csr', 'bsr'])
A = A.tobsr(blocksize=(blocksize, blocksize))
if Dinv is None:
Dinv = get_block_diag(A, blocksize=blocksize, inv_flag=True)
elif Dinv.shape[0] != A.shape[0]/blocksize:
raise ValueError('Dinv and A have incompatible dimensions')
elif (Dinv.shape[1] != blocksize) or (Dinv.shape[2] != blocksize):
raise ValueError('Dinv and blocksize are incompatible')
if sweep == 'forward':
row_start, row_stop, row_step = 0, len(x)/blocksize, 1
elif sweep == 'backward':
row_start, row_stop, row_step = len(x)/blocksize-1, -1, -1
elif sweep == 'symmetric':
for iter in xrange(iterations):
block_gauss_seidel(A, x, b, iterations=1, sweep='forward',
blocksize=blocksize, Dinv=Dinv)
block_gauss_seidel(A, x, b, iterations=1, sweep='backward',
blocksize=blocksize, Dinv=Dinv)
return
else:
raise ValueError("valid sweep directions are 'forward',\
'backward', and 'symmetric'")
for iter in xrange(iterations):
amg_core.block_gauss_seidel(A.indptr, A.indices, numpy.ravel(A.data),
x, b, numpy.ravel(Dinv),
row_start, row_stop, row_step, blocksize)
[docs]def polynomial(A, x, b, coefficients, iterations=1):
"""Apply a polynomial smoother to the system Ax=b
Parameters
----------
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
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'.
- 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).
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)
"""
A, x, b = make_system(A, x, b, formats=None)
for i in range(iterations):
from pyamg.util.linalg import norm
if norm(x) == 0:
residual = b
else:
residual = (b - A*x)
h = coefficients[0]*residual
for c in coefficients[1:]:
h = c*residual + A*h
x += h
[docs]def 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.
Parameters
----------
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
Returns
-------
Nothing, x will be modified in place.
Examples
--------
>>> 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
"""
A, x, b = make_system(A, x, b, formats=['csr'])
indices = numpy.asarray(indices, dtype='intc')
# if indices.min() < 0:
# raise ValueError('row index (%d) is invalid' % indices.min())
# if indices.max() >= A.shape[0]
# raise ValueError('row index (%d) is invalid' % indices.max())
if sweep == 'forward':
row_start, row_stop, row_step = 0, len(indices), 1
elif sweep == 'backward':
row_start, row_stop, row_step = len(indices)-1, -1, -1
elif sweep == 'symmetric':
for iter in xrange(iterations):
gauss_seidel_indexed(A, x, b, indices, iterations=1,
sweep='forward')
gauss_seidel_indexed(A, x, b, indices, iterations=1,
sweep='backward')
return
else:
raise ValueError("valid sweep directions are 'forward',\
'backward', and 'symmetric'")
for iter in xrange(iterations):
amg_core.gauss_seidel_indexed(A.indptr, A.indices, A.data,
x, b, indices,
row_start, row_stop, row_step)
[docs]def 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)
Parameters
----------
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
Returns
-------
Nothing, x will be modified in place.
References
----------
.. [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.
Examples
--------
>>> # 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)
"""
A, x, b = make_system(A, x, b, formats=['csr'])
sweep = slice(None)
(row_start, row_stop, row_step) = sweep.indices(A.shape[0])
temp = numpy.zeros_like(x)
# Dinv for A*A.H
Dinv = get_diagonal(A, norm_eq=2, inv=True)
# Create uniform type, convert possibly complex scalars to length 1 arrays
[omega] = type_prep(A.dtype, [omega])
for i in range(iterations):
delta = (numpy.ravel(b - A*x)*numpy.ravel(Dinv)).astype(A.dtype)
amg_core.jacobi_ne(A.indptr, A.indices, A.data,
x, b, delta, temp, row_start,
row_stop, row_step, omega)
[docs]def 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)
Parameters
----------
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)
Returns
-------
Nothing, x will be modified in place.
References
----------
.. [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
Examples
--------
>>> # 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)
"""
A, x, b = make_system(A, x, b, formats=['csr'])
# Dinv for A*A.H
if Dinv is None:
Dinv = numpy.ravel(get_diagonal(A, norm_eq=2, inv=True))
if sweep == 'forward':
row_start, row_stop, row_step = 0, len(x), 1
elif sweep == 'backward':
row_start, row_stop, row_step = len(x)-1, -1, -1
elif sweep == 'symmetric':
for iter in xrange(iterations):
gauss_seidel_ne(A, x, b, iterations=1, sweep='forward',
omega=omega, Dinv=Dinv)
gauss_seidel_ne(A, x, b, iterations=1, sweep='backward',
omega=omega, Dinv=Dinv)
return
else:
raise ValueError("valid sweep directions are 'forward',\
'backward', and 'symmetric'")
for i in xrange(iterations):
amg_core.gauss_seidel_ne(A.indptr, A.indices, A.data,
x, b, row_start,
row_stop, row_step, Dinv, omega)
[docs]def 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
Parameters
----------
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)
Returns
-------
Nothing, x will be modified in place.
References
----------
.. [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
Examples
--------
>>> # 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)
"""
A, x, b = make_system(A, x, b, formats=['csc'])
# Dinv for A.H*A
if Dinv is None:
Dinv = numpy.ravel(get_diagonal(A, norm_eq=1, inv=True))
if sweep == 'forward':
col_start, col_stop, col_step = 0, len(x), 1
elif sweep == 'backward':
col_start, col_stop, col_step = len(x)-1, -1, -1
elif sweep == 'symmetric':
for iter in xrange(iterations):
gauss_seidel_nr(A, x, b, iterations=1, sweep='forward',
omega=omega, Dinv=Dinv)
gauss_seidel_nr(A, x, b, iterations=1, sweep='backward',
omega=omega, Dinv=Dinv)
return
else:
raise ValueError("valid sweep directions are 'forward',\
'backward', and 'symmetric'")
# Calculate initial residual
r = b - A*x
for i in xrange(iterations):
amg_core.gauss_seidel_nr(A.indptr, A.indices, A.data,
x, r, col_start,
col_stop, col_step, Dinv, omega)
# from pyamg.utils import dispatcher
# dispatch = dispatcher( dict([ (fn,eval(fn)) for fn in __all__ ]) )