"""Method to create pre- and post-smoothers on the levels of a multilevel_solver
"""
import numpy
import scipy
import relaxation
from chebyshev import chebyshev_polynomial_coefficients
from pyamg.util.utils import scale_rows, get_block_diag, UnAmal, get_diagonal
from pyamg.util.linalg import approximate_spectral_radius
from pyamg.krylov import gmres, cgne, cgnr, cg
from pyamg import amg_core
from scipy.linalg import lapack as la
__docformat__ = "restructuredtext en"
__all__ = ['change_smoothers']
def unpack_arg(v):
if isinstance(v,tuple):
return v[0],v[1]
else:
return v,{}
[docs]def change_smoothers(ml, presmoother, postsmoother):
'''
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}
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
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)
'''
# interpret arguments into list
if isinstance(presmoother, str) or isinstance(presmoother, tuple) or (presmoother == None):
presmoother = [presmoother]
elif not isinstance(presmoother, list):
raise ValueError,'Unrecognized presmoother'
if isinstance(postsmoother, str) or isinstance(postsmoother, tuple) or (postsmoother == None):
postsmoother = [postsmoother]
elif not isinstance(postsmoother, list):
raise ValueError,'Unrecognized postsmoother'
# set ml.levels[i].presmoother = presmoother[i]
i = 0
for i in range( min(len(presmoother), len(ml.levels[:-1])) ):
# unpack presmoother[i]
fn,kwargs = unpack_arg(presmoother[i])
# get function handle
try:
setup_presmoother = eval('setup_' + str(fn))
except NameError, ne:
raise NameError("invalid presmoother method: ", fn)
ml.levels[i].presmoother = setup_presmoother(ml.levels[i], **kwargs)
# Fill in remaining levels
for j in range(i+1, len(ml.levels[:-1])):
ml.levels[j].presmoother = setup_presmoother(ml.levels[j], **kwargs)
# set ml.levels[i].postsmoother = postsmoother[i]
i = 0
for i in range( min(len(postsmoother), len(ml.levels[:-1])) ):
# unpack postsmoother[i]
fn,kwargs = unpack_arg(postsmoother[i])
# get function handle
try:
setup_postsmoother = eval('setup_' + str(fn))
except NameError, ne:
raise NameError("invalid postsmoother method: ", fn)
ml.levels[i].postsmoother = setup_postsmoother(ml.levels[i], **kwargs)
# Fill in remaining levels
for j in range(i+1, len(ml.levels[:-1])):
ml.levels[j].postsmoother = setup_postsmoother(ml.levels[j], **kwargs)
def rho_D_inv_A(A):
"""
Return the (approx.) spectral radius of D^-1 * A
Parameters
----------
A : {sparse-matrix}
Returns
-------
approximate spectral radius of diag(A)^{-1} A
Examples
--------
>>> from pyamg.gallery import poisson
>>> from pyamg.relaxation.smoothing import rho_D_inv_A
>>> from scipy.sparse import csr_matrix
>>> import numpy
>>> A = csr_matrix(numpy.array([[1,0,0],[0,2,0],[0,0,3]]))
>>> print rho_D_inv_A(A)
1.0
"""
if not hasattr(A, 'rho_D_inv'):
D_inv = get_diagonal(A, inv=True)
D_inv_A = scale_rows(A, D_inv, copy=True)
A.rho_D_inv = approximate_spectral_radius(D_inv_A)
return A.rho_D_inv
def rho_block_D_inv_A(A, Dinv):
"""
Return the (approx.) spectral radius of block D^-1 * A
Parameters
----------
A : {sparse-matrix}
size NxN
Dinv : {array}
Inverse of diagonal blocks of A
size (N/blocksize, blocksize, blocksize)
Returns
-------
approximate spectral radius of (Dinv A)
Examples
--------
>>> from pyamg.gallery import poisson
>>> from pyamg.relaxation.smoothing import rho_block_D_inv_A
>>> from pyamg.util.utils import get_block_diag
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> Dinv = get_block_diag(A, blocksize=4, inv_flag=True)
"""
if not hasattr(A, 'rho_block_D_inv'):
from scipy.sparse.linalg import LinearOperator
blocksize = Dinv.shape[1]
if Dinv.shape[1] != Dinv.shape[2]:
raise ValueError('Dinv has incorrect dimensions')
elif Dinv.shape[0] != A.shape[0]/blocksize:
raise ValueError('Dinv and A have incompatible dimensions')
Dinv = scipy.sparse.bsr_matrix( (Dinv, \
scipy.arange(Dinv.shape[0]), scipy.arange(Dinv.shape[0]+1)), shape=A.shape)
# Don't explicitly form Dinv*A
def matvec(x):
return Dinv*(A*x)
D_inv_A = LinearOperator(A.shape, matvec, dtype=A.dtype)
A.rho_block_D_inv = approximate_spectral_radius(D_inv_A)
return A.rho_block_D_inv
def matrix_asformat(lvl, name, format, blocksize=None):
'''
This routine looks for the matrix "name" in the specified format as a
member of the level instance, lvl. For example, if name='A', format='bsr'
and blocksize=(4,4), and if lvl.Absr44 exists with the correct blocksize,
then lvl.Absr is returned. If the matrix doesn't already exist, lvl.name
is converted to the desired format, and made a member of lvl.
Only create such persistent copies of a matrix for routines such as
presmoothing and postsmoothing, where the matrix conversion is done every
cycle.
Calling this function can _dramatically_ increase your memory costs.
Be careful with it's usage.
'''
desired_matrix = 'lvl.' + name + format
if format == 'bsr':
desired_matrix += str(blocksize[0])+str(blocksize[1])
base_matrix = 'lvl.' + name
if hasattr(lvl, desired_matrix[4:]):
# if lvl already contains lvl.name+format
pass
elif eval(base_matrix).format == format and format != 'bsr':
# is base_matrix already in the correct format?
exec desired_matrix +' = '+ base_matrix
elif eval(base_matrix).format == format and format == 'bsr':
# make sure blocksize is correct
if eval(base_matrix).blocksize != blocksize:
exec desired_matrix +' = '+ base_matrix+'.to'+format+'(blocksize='+str(blocksize)+')'
else:
exec desired_matrix +' = '+ base_matrix
elif format == 'bsr':
# convert
exec desired_matrix +' = '+ base_matrix+'.to'+format+'(blocksize='+str(blocksize)+')'
else:
# convert
exec desired_matrix + ' = lvl.' + name + '.to' + format + '()'
return eval(desired_matrix)
#Appends
#A.schwarz_tuple = (subdomain, subdomain_ptr, inv_subblock, inv_subblock_ptr) to
#the matrix passed in If this tuple already exists, return it, otherwise compute
#it Should work for passing in both A and the Strength matrix C Make sure to
#wrap an Acsr into the schwarz call
#make sure that it handles passing in preset subdomain stuff, and recomputing it
#If subdomain and subdomain_ptr are passed in, check to see that they are the
#same as any preexisting subdomain and subdomain_ptr?
def schwarz_parameters(A, subdomain=None, subdomain_ptr=None,
inv_subblock=None, inv_subblock_ptr=None):
'''
Helper function for setting up Schwarz relaxation. This function avoids
recomputing the subdomains and block inverses manytimes, e.g., it avoids
a costly double computation when setting up pre and post smoothing with Schwarz.
Parameters
----------
A {csr_matrix}
Returns
-------
A.schwarz_parameters[0] is subdomain
A.schwarz_parameters[1] is subdomain_ptr
A.schwarz_parameters[2] is inv_subblock
A.schwarz_parameters[3] is inv_subblock_ptr
'''
# Check if A has a pre-existing set of Schwarz parameters
if hasattr(A, 'schwarz_parameters'):
if subdomain != None and subdomain_ptr != None:
# check that the existing parameters correspond to the same subdomains
if numpy.array(A.schwarz_parameters[0] == subdomain).all() and \
numpy.array(A.schwarz_parameters[1] == subdomain_ptr).all():
return A.schwarz_parameters
else:
return A.schwarz_parameters
# Default is to use the overlapping regions defined by A's sparsity pattern
if subdomain is None or subdomain_ptr is None:
subdomain_ptr = A.indptr.copy()
subdomain = A.indices.copy()
##
# Extract each subdomain's block from the matrix
if inv_subblock is None or inv_subblock_ptr is None:
inv_subblock_ptr = numpy.zeros(subdomain_ptr.shape, dtype=A.indices.dtype)
blocksize = (subdomain_ptr[1:] - subdomain_ptr[:-1])
inv_subblock_ptr[1:] = numpy.cumsum(blocksize*blocksize)
##
# Extract each block column from A
inv_subblock = numpy.zeros((inv_subblock_ptr[-1],), dtype=A.dtype)
amg_core.extract_subblocks(A.indptr, A.indices, A.data, inv_subblock,
inv_subblock_ptr, subdomain, subdomain_ptr,
int(subdomain_ptr.shape[0]-1), A.shape[0])
##
# Choose tolerance for which singular values are zero in *gelss below
t = A.dtype.char
eps = numpy.finfo(numpy.float).eps
feps = numpy.finfo(numpy.single).eps
geps = numpy.finfo(numpy.longfloat).eps
_array_precision = {'f': 0, 'd': 1, 'g': 2, 'F': 0, 'D': 1, 'G':2}
cond = {0: feps*1e3, 1: eps*1e6, 2: geps*1e6}[_array_precision[t]]
##
# Invert each block column
my_pinv, = la.get_lapack_funcs(['gelss'], (numpy.ones((1,), dtype=A.dtype)) )
for i in xrange(subdomain_ptr.shape[0]-1):
m = blocksize[i]
rhs = scipy.eye(m,m, dtype=A.dtype)
gelssoutput = my_pinv(inv_subblock[inv_subblock_ptr[i]:inv_subblock_ptr[i+1]].reshape(m,m),
rhs, cond=cond, overwrite_a=True, overwrite_b=True)
inv_subblock[inv_subblock_ptr[i]:inv_subblock_ptr[i+1]] = numpy.ravel(gelssoutput[1])
A.schwarz_parameters = (subdomain, subdomain_ptr, inv_subblock, inv_subblock_ptr)
return A.schwarz_parameters
"""
The following setup_smoother_name functions are helper functions for
parsing user input and assigning each level the appropriate smoother for
the above functions "change_smoothers".
The standard interface is
Parameters
----------
lvl : {multilevel level}
the level in the hierarchy for which to assign a smoother
iterations : {int}
how many smoother iterations
optional_params : {}
optional params specific for each method such as omega or sweep
Returns
-------
Function pointer for the appropriate relaxation method for level=lvl
Examples
--------
See change_smoothers above
"""
def setup_gauss_seidel(lvl, iterations=1, sweep='forward'):
def smoother(A,x,b):
relaxation.gauss_seidel(A, x, b, iterations=iterations, sweep=sweep)
return smoother
def setup_jacobi(lvl, iterations=1, omega=1.0, withrho=True):
if withrho:
omega = omega/rho_D_inv_A(lvl.A)
def smoother(A,x,b):
relaxation.jacobi(A, x, b, iterations=iterations, omega=omega)
return smoother
def setup_schwarz(lvl, iterations=1, subdomain=None, subdomain_ptr=None, \
inv_subblock=None, inv_subblock_ptr=None, sweep='symmetric'):
Acsr = matrix_asformat(lvl, 'A', 'csr')
subdomain, subdomain_ptr, inv_subblock, inv_subblock_ptr = \
schwarz_parameters(Acsr, subdomain, subdomain_ptr, inv_subblock, inv_subblock_ptr)
def smoother(A,x,b):
relaxation.schwarz(A, x, b, iterations=iterations, subdomain=subdomain, \
subdomain_ptr=subdomain_ptr, inv_subblock=inv_subblock,\
inv_subblock_ptr=inv_subblock_ptr, sweep=sweep)
return smoother
def setup_strength_based_schwarz(lvl, iterations=1, sweep='symmetric'):
# Use the overlapping regions defined by the strength of connection matrix C
# for the overlapping Schwarz method
if not hasattr(lvl,'C'):
C = lvl.A.tocsr()
else:
C = lvl.C.tocsr()
subdomain_ptr = C.indptr.copy()
subdomain = C.indices.copy()
return setup_schwarz(lvl, iterations=iterations, subdomain=subdomain, \
subdomain_ptr=subdomain_ptr, sweep=sweep)
def setup_block_jacobi(lvl, iterations=1, omega=1.0, Dinv=None, blocksize=None, withrho=True):
##
# Determine Blocksize
if blocksize == None and Dinv == None:
if scipy.sparse.isspmatrix_csr(lvl.A):
blocksize = 1
elif scipy.sparse.isspmatrix_bsr(lvl.A):
blocksize = lvl.A.blocksize[0]
elif blocksize == None:
blocksize = Dinv.shape[1]
if blocksize == 1:
# Block Jacobi is equivalent to normal Jacobi
return setup_jacobi(lvl, iterations=iterations, omega=omega, withrho=withrho)
else:
# Use Block Jacobi
if Dinv == None:
Dinv = get_block_diag(lvl.A, blocksize=blocksize, inv_flag=True)
if withrho:
omega = omega/rho_block_D_inv_A(lvl.A, Dinv)
def smoother(A,x,b):
relaxation.block_jacobi(A, x, b, iterations=iterations, omega=omega, \
Dinv=Dinv, blocksize=blocksize)
return smoother
def setup_block_gauss_seidel(lvl, iterations=1, sweep='forward', Dinv=None, blocksize=None):
##
# Determine Blocksize
if blocksize == None and Dinv == None:
if scipy.sparse.isspmatrix_csr(lvl.A):
blocksize = 1
elif scipy.sparse.isspmatrix_bsr(lvl.A):
blocksize = lvl.A.blocksize[0]
elif blocksize == None:
blocksize = Dinv.shape[1]
if blocksize == 1:
# Block GS is equivalent to normal GS
return setup_gauss_seidel(lvl, iterations=iterations, sweep=sweep)
else:
# Use Block GS
if Dinv == None:
Dinv = get_block_diag(lvl.A, blocksize=blocksize, inv_flag=True)
def smoother(A,x,b):
relaxation.block_gauss_seidel(A, x, b, iterations=iterations, \
Dinv=Dinv, blocksize=blocksize, sweep=sweep)
return smoother
def setup_richardson(lvl, iterations=1, omega=1.0):
omega = omega/approximate_spectral_radius(lvl.A)
def smoother(A,x,b):
relaxation.polynomial(A, x, b, coefficients=[omega], iterations=iterations)
return smoother
def setup_sor(lvl, omega=0.5, iterations=1, sweep='forward'):
def smoother(A,x,b):
relaxation.sor(A, x, b, omega=omega, iterations=iterations, sweep=sweep)
return smoother
def setup_chebyshev(lvl, lower_bound=1.0/30.0, upper_bound=1.1, degree=3, iterations=1):
rho = approximate_spectral_radius(lvl.A)
a = rho * lower_bound
b = rho * upper_bound
coefficients = -chebyshev_polynomial_coefficients(a, b, degree)[:-1] # drop the constant coefficient
def smoother(A,x,b):
relaxation.polynomial(A, x, b, coefficients=coefficients, iterations=iterations)
return smoother
def setup_jacobi_ne(lvl, iterations=1, omega=1.0, withrho=True):
Acsr = matrix_asformat(lvl, 'A', 'csr')
if withrho:
omega = omega/rho_D_inv_A(Acsr)**2
def smoother(A,x,b):
relaxation.jacobi_ne(Acsr, x, b, iterations=iterations, omega=omega)
return smoother
def setup_gauss_seidel_ne(lvl, iterations=1, sweep='forward', omega=1.0):
Acsr = matrix_asformat(lvl, 'A', 'csr')
def smoother(A,x,b):
relaxation.gauss_seidel_ne(Acsr, x, b, iterations=iterations, sweep=sweep, omega=omega)
return smoother
def setup_gauss_seidel_nr(lvl, iterations=1, sweep='forward', omega=1.0):
Acsc = matrix_asformat(lvl, 'A', 'csc')
def smoother(A,x,b):
relaxation.gauss_seidel_nr(Acsc, x, b, iterations=iterations, sweep=sweep, omega=omega)
return smoother
def setup_gmres(lvl, tol=1e-12, maxiter=1, restrt=None, M=None, callback=None, residuals=None):
def smoother(A,x,b):
x[:] = (gmres(A, b, x0=x, tol=tol, maxiter=maxiter, restrt=restrt, M=M, callback=callback, residuals=residuals)[0]).reshape(x.shape)
return smoother
def setup_cg(lvl, tol=1e-12, maxiter=1, M=None, callback=None, residuals=None):
def smoother(A,x,b):
x[:] = (cg(A, b, x0=x, tol=tol, maxiter=maxiter, M=M, callback=callback, residuals=residuals)[0]).reshape(x.shape)
return smoother
def setup_cgne(lvl, tol=1e-12, maxiter=1, M=None, callback=None, residuals=None):
def smoother(A,x,b):
x[:] = (cgne(A, b, x0=x, tol=tol, maxiter=maxiter, M=M, callback=callback, residuals=residuals)[0]).reshape(x.shape)
return smoother
def setup_cgnr(lvl, tol=1e-12, maxiter=1, M=None, callback=None, residuals=None):
def smoother(A,x,b):
x[:] = (cgnr(A, b, x0=x, tol=tol, maxiter=maxiter, M=M, callback=callback, residuals=residuals)[0]).reshape(x.shape)
return smoother
def setup_None(lvl):
def smoother(A,x,b):
pass
return smoother