util Package

util Package

linalg.py provides some linear algebra functionality not yet found in scipy

utils.py provides some utility functions for use with pyamg

BSR_utils.py provides utility functions for accessing and writing individual rows of BSR matrices

pyamg.util.Coord2RBM(numNodes, numPDEs, x, y, z)

Convert 2D or 3D coordinates into Rigid body modes for use as near nullspace modes in elasticity AMG solvers

numNodes : int
Number of nodes
numPDEs :
Number of dofs per node
x,y,z : array_like
Coordinate vectors
rbm : matrix
A matrix of size (numNodes*numPDEs) x (1 | 6) containing the 6 rigid body modes
>>> import numpy
>>> from pyamg.util.utils import Coord2RBM
>>> a = numpy.array([0,1,2])
>>> Coord2RBM(3,6,a,a,a)
matrix([[ 1.,  0.,  0.,  0.,  0., -0.],
        [ 0.,  1.,  0., -0.,  0.,  0.],
        [ 0.,  0.,  1.,  0., -0.,  0.],
        [ 0.,  0.,  0.,  1.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  1.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  1.],
        [ 1.,  0.,  0.,  0.,  1., -1.],
        [ 0.,  1.,  0., -1.,  0.,  1.],
        [ 0.,  0.,  1.,  1., -1.,  0.],
        [ 0.,  0.,  0.,  1.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  1.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  1.],
        [ 1.,  0.,  0.,  0.,  2., -2.],
        [ 0.,  1.,  0., -2.,  0.,  2.],
        [ 0.,  0.,  1.,  2., -2.,  0.],
        [ 0.,  0.,  0.,  1.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  1.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  1.]])
pyamg.util.UnAmal(A, RowsPerBlock, ColsPerBlock)

Unamalgamate a CSR A with blocks of 1’s. This operation is equivalent to replacing each entry of A with ones(RowsPerBlock, ColsPerBlock), i.e., this is equivalent to setting all of A’s nonzeros to 1 and then doing a Kronecker product between A and ones(RowsPerBlock, ColsPerBlock).

A : csr_matrix
Amalgamted matrix
RowsPerBlock : int
Give A blocks of size (RowsPerBlock, ColsPerBlock)
ColsPerBlock : int
Give A blocks of size (RowsPerBlock, ColsPerBlock)
A : bsr_matrix
Returns A.data[:] = 1, followed by a Kronecker product of A and ones(RowsPerBlock, ColsPerBlock)
>>> from numpy import array
>>> from scipy.sparse import csr_matrix
>>> from pyamg.util.utils import UnAmal
>>> row = array([0,0,1,2,2,2])
>>> col = array([0,2,2,0,1,2])
>>> data = array([1,2,3,4,5,6])
>>> A = csr_matrix( (data,(row,col)), shape=(3,3) )
>>> A.todense()
matrix([[1, 0, 2],
        [0, 0, 3],
        [4, 5, 6]])
>>> UnAmal(A,2,2).todense()
matrix([[ 1.,  1.,  0.,  0.,  1.,  1.],
        [ 1.,  1.,  0.,  0.,  1.,  1.],
        [ 0.,  0.,  0.,  0.,  1.,  1.],
        [ 0.,  0.,  0.,  0.,  1.,  1.],
        [ 1.,  1.,  1.,  1.,  1.,  1.],
        [ 1.,  1.,  1.,  1.,  1.,  1.]])
pyamg.util.amalgamate(A, blocksize)

Amalgamate matrix A

A : csr_matrix
Matrix to amalgamate
blocksize : int
blocksize to use while amalgamating
A_amal : csr_matrix
Amalgamated matrix A, first, convert A to BSR with square blocksize and then return a CSR matrix of ones using the resulting BSR indptr and indices

inverse operation of UnAmal for square matrices

>>> from numpy import array
>>> from scipy.sparse import csr_matrix
>>> from pyamg.util.utils import amalgamate
>>> row = array([0,0,1])
>>> col = array([0,2,1])
>>> data = array([1,2,3])
>>> A = csr_matrix( (data,(row,col)), shape=(4,4) )
>>> A.todense()
matrix([[1, 0, 2, 0],
        [0, 3, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]])
>>> amalgamate(A,2).todense()
matrix([[ 1.,  1.],
        [ 0.,  0.]])
pyamg.util.approximate_spectral_radius(A, tol=0.01, maxiter=15, restart=5, symmetric=None, initial_guess=None, return_vector=False)

Approximate the spectral radius of a matrix

A : {dense or sparse matrix}
E.g. csr_matrix, csc_matrix, ndarray, etc.
tol : {scalar}
Relative tolerance of approximation, i.e., the error divided by the approximate spectral radius is compared to tol.
maxiter : {integer}
Maximum number of iterations to perform
restart : {integer}
Number of restarted Arnoldi processes. For example, a value of 0 will run Arnoldi once, for maxiter iterations, and a value of 1 will restart Arnoldi once, using the maximal eigenvector from the first Arnoldi process as the initial guess.
symmetric : {boolean}
True - if A is symmetric
Lanczos iteration is used (more efficient)
False - if A is non-symmetric (default
Arnoldi iteration is used (less efficient)
initial_guess : {array|None}
If n x 1 array, then use as initial guess for Arnoldi/Lanczos. If None, then use a random initial guess.
return_vector : {boolean}
True - return an approximate dominant eigenvector, in addition to the
spectral radius.

False - Do not return the approximate dominant eigenvector

An approximation to the spectral radius of A, and if return_vector=True, then also return the approximate dominant eigenvector

The spectral radius is approximated by looking at the Ritz eigenvalues. Arnoldi iteration (or Lanczos) is used to project the matrix A onto a Krylov subspace: H = Q* A Q. The eigenvalues of H (i.e. the Ritz eigenvalues) should represent the eigenvalues of A in the sense that the minimum and maximum values are usually well matched (for the symmetric case it is true since the eigenvalues are real).

[1]Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, editors. “Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide”, SIAM, Philadelphia, 2000.
>>> from pyamg.util.linalg import approximate_spectral_radius
>>> import numpy
>>> from scipy.linalg import eigvals, norm
>>> A = numpy.array([[1.,0.],[0.,1.]])
>>> print approximate_spectral_radius(A,maxiter=3)
1.0
>>> print max([norm(x) for x in eigvals(A)])
1.0
pyamg.util.blocksize(A)
pyamg.util.compute_BtBinv(B, C)
Helper function that creates inv(B_i.T B_i) for each block row i in C,
where B_i is B restricted to the sparsity pattern of block row i.
B : {array}
(M,k) array, typically near-nullspace modes for coarse grid, i.e., B_c.
C : {csr_matrix, bsr_matrix}
Sparse NxM matrix, whose sparsity structure (i.e., matrix graph) is used to determine BtBinv.
BtBinv : {array}
BtBinv[i] = inv(B_i.T B_i), where B_i is B restricted to the nonzero pattern of block row i in C.
>>> from numpy import array
>>> from scipy.sparse import bsr_matrix
>>> from pyamg.util.utils import compute_BtBinv
>>> T = array([[ 1.,  0.],
...            [ 1.,  0.],
...            [ 0.,  .5],
...            [ 0.,  .25]])
>>> T = bsr_matrix(T)
>>> B = array([[1.],[2.]])
>>> compute_BtBinv(B, T)
array([[[ 1.  ]],

       [[ 1.  ]],

       [[ 0.25]],

       [[ 0.25]]])

The principal calling routines are aggregation.smooth.energy_prolongation_smoother, and util.utils.filter_operator.

BtBinv is used in the prolongation smoothing process that incorporates B into the span of prolongation with row-wise projection operators. It is these projection operators that BtBinv is part of.

pyamg.util.cond(A)

Returns condition number of A

A : {dense or sparse matrix}
e.g. array, matrix, csr_matrix, ...

2-norm condition number through use of the SVD Use for small to moderate sized dense matrices. For large sparse matrices, use condest.

The condition number measures how large of a change in the problems solution is caused by a change in problem’s input. Large condition numbers indicate that small perturbations and numerical errors are magnified greatly when solving the system.

>>> import numpy
>>> from pyamg.util.linalg import condest
>>> c = condest(numpy.array([[1.0,0.],[0.,2.0]]))
>>> print c
2.0
pyamg.util.condest(A, tol=0.1, maxiter=25, symmetric=False)

Estimates the condition number of A

A : {dense or sparse matrix}
e.g. array, matrix, csr_matrix, ...
tol : {float}
Approximation tolerance, currently not used
maxiter: {int}
Max number of Arnoldi/Lanczos iterations
symmetric : {bool}
If symmetric use the far more efficient Lanczos algorithm, Else use Arnoldi

Estimate of cond(A) with |lambda_max| / |lambda_min| through the use of Arnoldi or Lanczos iterations, depending on the symmetric flag

The condition number measures how large of a change in the the problems solution is caused by a change in problem’s input. Large condition numbers indicate that small perturbations and numerical errors are magnified greatly when solving the system.

>>> import numpy
>>> from pyamg.util.linalg import condest
>>> c = condest(numpy.array([[1.,0.],[0.,2.]]))
>>> print c
2.0
pyamg.util.diag_sparse(A)
If A is a sparse matrix (e.g. csr_matrix or csc_matrix)
  • return the diagonal of A as an array
Otherwise
  • return a csr_matrix with A on the diagonal
A : sparse matrix or 1d array
General sparse matrix or array of diagonal entries
B : array or sparse matrix
Diagonal sparse is returned as csr if A is dense otherwise return an array of the diagonal
>>> import numpy
>>> from pyamg.util.utils import diag_sparse
>>> d = 2.0*numpy.ones((3,)).ravel()
>>> print diag_sparse(d).todense()
[[ 2.  0.  0.]
 [ 0.  2.  0.]
 [ 0.  0.  2.]]
pyamg.util.eliminate_diag_dom_nodes(A, C, theta=1.02)

Helper function that eliminates diagonally dominant rows and cols from A in the separate matrix C. This is useful because it eliminates nodes in C which we don’t want coarsened. These eliminated nodes in C just become the rows and columns of the identity.

A : {csr_matrix, bsr_matrix}
Sparse NxN matrix
C : {csr_matrix}
Sparse MxM matrix, where M is the number of nodes in A. M=N if A is CSR or is BSR with blocksize 1. Otherwise M = N/blocksize.
theta : {float}
determines diagonal dominance threshhold
C : {csr_matrix}
C updated such that the rows and columns corresponding to diagonally dominant rows in A have been eliminated and replaced with rows and columns of the identity.
Diagonal dominance is defined as
|| (e_i, A) - a_ii ||_1 < theta a_ii

that is, the 1-norm of the off diagonal elements in row i must be less than theta times the diagonal element.

>>> from pyamg.gallery import poisson
>>> from pyamg.util.utils import eliminate_diag_dom_nodes
>>> A = poisson( (4,), format='csr' )
>>> C = eliminate_diag_dom_nodes(A, A.copy(), 1.1)
>>> C.todense()
matrix([[ 1.,  0.,  0.,  0.],
        [ 0.,  2., -1.,  0.],
        [ 0., -1.,  2.,  0.],
        [ 0.,  0.,  0.,  1.]])
pyamg.util.filter_operator(A, C, B, Bf, BtBinv=None)

Filter the matrix A according to the matrix graph of C, while ensuring that the new, filtered A satisfies: A_new*B = Bf.

A : {csr_matrix, bsr_matrix}
n x m matrix to filter
C : {csr_matrix, bsr_matrix}
n x m matrix representing the couplings in A to keep
B : {array}
m x k array of near nullspace vectors
Bf : {array}
n x k array of near nullspace vectors to place in span(A)
BtBinv : {None, array}
3 dimensional array such that, BtBinv[i] = pinv(B_i.H Bi), and B_i is B restricted to the neighborhood (with respect to the matrix graph of C) of dof of i. If None is passed in, this array is computed internally.

A : sparse matrix updated such that sparsity structure of A now matches that of C, and that the relationship A*B = Bf holds.

This procedure allows for certain important modes (i.e., Bf) to be placed in span(A) by way of row-wise l2-projections that enforce the relationship A*B = Bf. This is useful for maintaining certain modes (e.g., the constant) in the span of prolongation.

>>> from numpy import ones, array
>>> from scipy.sparse import csr_matrix
>>> from pyamg.util.utils import filter_operator
>>> A = array([ [1.,1,1],[1,1,1],[0,1,0],[0,1,0],[0,0,1],[0,0,1]])
>>> C = array([ [1.,1,0],[1,1,0],[0,1,0],[0,1,0],[0,0,1],[0,0,1]])
>>> B = ones((3,1))
>>> Bf = ones((6,1))
>>> filter_operator(csr_matrix(A), csr_matrix(C), B, Bf).todense()
matrix([[ 0.5,  0.5,  0. ],
        [ 0.5,  0.5,  0. ],
        [ 0. ,  1. ,  0. ],
        [ 0. ,  1. ,  0. ],
        [ 0. ,  0. ,  1. ],
        [ 0. ,  0. ,  1. ]])

This routine is primarily used in pyamg.aggregation.smooth.energy_prolongation_smoother, where it is used to generate a suitable initial guess for the energy-minimization process, when root-node style SA is used. Essentially, the tentative prolongator, T, is processed by this routine to produce fine-grid nullspace vectors when multiplying coarse-grid nullspace vectors, i.e., T*B = Bf. This is possible for any arbitrary vectors B and Bf, so long as the sparsity structure of T is rich enough.

When generating initial guesses for root-node style prolongation operators, this function is usually called before pyamg.uti.utils.scale_T

pyamg.util.get_Cpt_params(A, Cnodes, AggOp, T)
Helper function that returns a dictionary of sparse matrices and arrays
which allow us to easily operate on Cpts and Fpts separately.
A : {csr_matrix, bsr_matrix}
Operator
Cnodes : {array}
Array of all root node indices. This is an array of nodal indices, not degree-of-freedom indices. If the blocksize of T is 1, then nodal indices and degree-of-freedom indices coincide.
AggOp : {csr_matrix}
Aggregation operator corresponding to A
T : {bsr_matrix}
Tentative prolongator based on AggOp

Dictionary containing these parameters:

P_I : {bsr_matrix}
Interpolation operator that carries out only simple injection from the coarse grid to fine grid Cpts nodes
I_F : {bsr_matrix}
Identity operator on Fpts, i.e., the action of this matrix zeros out entries in a vector at all Cpts, leaving Fpts untouched
I_C : {bsr_matrix}
Identity operator on Cpts nodes, i.e., the action of this matrix zeros out entries in a vector at all Fpts, leaving Cpts untouched
Cpts : {array}
An array of all root node dofs, corresponding to the F/C splitting
Fpts : {array}
An array of all non root node dofs, corresponding to the F/C splitting
>>> from numpy import array
>>> from pyamg.util.utils import get_Cpt_params
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import csr_matrix, bsr_matrix
>>> A = poisson((10,), format='csr')
>>> Cpts = array([3, 7])
>>> AggOp = ([[ 1., 0.], [ 1., 0.],
...           [ 1., 0.], [ 1., 0.],
...           [ 1., 0.], [ 0., 1.],
...           [ 0., 1.], [ 0., 1.],
...           [ 0., 1.], [ 0., 1.]])
>>> AggOp = csr_matrix(AggOp)
>>> T = AggOp.copy().tobsr()
>>> params = get_Cpt_params(A, Cpts, AggOp, T)
>>> params['P_I'].todense()
matrix([[ 0.,  0.],
        [ 0.,  0.],
        [ 0.,  0.],
        [ 1.,  0.],
        [ 0.,  0.],
        [ 0.,  0.],
        [ 0.,  0.],
        [ 0.,  1.],
        [ 0.,  0.],
        [ 0.,  0.]])

The principal calling routine is aggregation.smooth.energy_prolongation_smoother, which uses the Cpt_param dictionary for root-node style prolongation smoothing

pyamg.util.get_block_diag(A, blocksize, inv_flag=True)

Return the block diagonal of A, in array form

A : csr_matrix
assumed to be square
blocksize : int
square block size for the diagonal
inv_flag : bool
if True, return the inverse of the block diagonal
block_diag : array
block diagonal of A in array form, array size is (A.shape[0]/blocksize, blocksize, blocksize)
>>> from scipy import arange
>>> from scipy.sparse import csr_matrix
>>> from pyamg.util import get_block_diag
>>> A = csr_matrix(arange(36).reshape(6,6))
>>> block_diag_inv = get_block_diag(A, blocksize=2, inv_flag=False)
>>> print block_diag_inv
[[[  0.   1.]
  [  6.   7.]]

 [[ 14.  15.]
  [ 20.  21.]]

 [[ 28.  29.]
  [ 34.  35.]]]
>>> block_diag_inv = get_block_diag(A, blocksize=2, inv_flag=True)
pyamg.util.get_diagonal(A, norm_eq=False, inv=False)
Return the diagonal or inverse of diagonal for
A, (A.H A) or (A A.H)
A : {dense or sparse matrix}
e.g. array, matrix, csr_matrix, ...
norm_eq : {0, 1, 2}
0 ==> D = diag(A) 1 ==> D = diag(A.H A) 2 ==> D = diag(A A.H)
inv : {True, False}
If True, D = 1.0/D

diagonal, D, of appropriate system

This function is especially useful for its fast methods of obtaining diag(A A.H) and diag(A.H A). Dinv is zero wherever D is zero

>>> from pyamg.util.utils import get_diagonal
>>> from pyamg.gallery import poisson
>>> A = poisson( (5,), format='csr' )
>>> D = get_diagonal(A)
>>> print D
[ 2.  2.  2.  2.  2.]
>>> D = get_diagonal(A, norm_eq=1, inv=True)
>>> print D
[ 0.2         0.16666667  0.16666667  0.16666667  0.2       ]
pyamg.util.hierarchy_spectrum(mg, filter=True, plot=False)

Examine a multilevel hierarchy’s spectrum

mg { pyamg multilevel hierarchy }
e.g. generated with smoothed_aggregation_solver(...) or ruge_stuben_solver(...)
  1. table to standard out detailing the spectrum of each level in mg
  2. if plot==True, a sequence of plots in the complex plane of the spectrum at each level

This can be useful for troubleshooting and when examining how your problem’s nature changes from level to level

>>> from pyamg import smoothed_aggregation_solver
>>> from pyamg.gallery import poisson
>>> from pyamg.util.utils import hierarchy_spectrum
>>> A = poisson( (1,), format='csr' )
>>> ml = smoothed_aggregation_solver(A)
>>> hierarchy_spectrum(ml)

 Level min(re(eig)) max(re(eig)) num re(eig) < 0 num re(eig) > 0 cond_2(A)
---------------------------------------------------------------------------
   0      2.000        2.000            0               1         1.00e+00


 Level min(im(eig)) max(im(eig)) num im(eig) < 0 num im(eig) > 0 cond_2(A)
---------------------------------------------------------------------------
   0      0.000        0.000            0               0         1.00e+00
pyamg.util.infinity_norm(A)

Infinity norm of a matrix (maximum absolute row sum).

A : csr_matrix, csc_matrix, sparse, or numpy matrix
Sparse or dense matrix
n : float
Infinity norm of the matrix
  • This serves as an upper bound on spectral radius.
  • csr and csc avoid a deep copy
  • dense calls scipy.linalg.norm

scipy.linalg.norm : dense matrix norms

>>> import numpy
>>> from scipy.sparse import spdiags
>>> from pyamg.util.linalg import infinity_norm
>>> n=10
>>> e = numpy.ones((n,1)).ravel()
>>> data = [ -1*e, 2*e, -1*e ]
>>> A = spdiags(data,[-1,0,1],n,n)
>>> print infinity_norm(A)
4.0
pyamg.util.ishermitian(A, fast_check=True, tol=1e-06, verbose=False)

Returns True if A is Hermitian to within tol

A : {dense or sparse matrix}
e.g. array, matrix, csr_matrix, ...
fast_check : {bool}
If True, use the heuristic < Ax, y> = < x, Ay> for random vectors x and y to check for conjugate symmetry. If False, compute A - A.H.
tol : {float}
Symmetry tolerance
verbose: {bool}
prints max( |A - A.H| ) if nonhermitian and fast_check=False abs( <Ax, y> - <x, Ay> ) if nonhermitian and fast_check=False

True if hermitian False if nonhermitian

This function applies a simple test of conjugate symmetry

>>> import numpy
>>> from pyamg.util.linalg import ishermitian
>>> ishermitian(numpy.array([[1,2],[1,1]]))
False
>>> from pyamg.gallery import poisson
>>> ishermitian(poisson((10,10)))
True
pyamg.util.levelize_smooth_or_improve_candidates(to_levelize, max_levels)

Helper function to preprocess the smooth and improve_candidates parameters passed to smoothed_aggregation_solver and rootnode_solver.

to_levelize : {string, tuple, list}
Parameter to preprocess, i.e., levelize and convert to a level-by-level list such that entry i specifies the parameter at level i
max_levels : int
Defines the maximum number of levels considered
to_levelize : list
The parameter list such that entry i specifies the parameter choice at level i.

This routine is needed because the user will pass in a parameter option such as smooth=’jacobi’, or smooth=[‘jacobi’, None], and this option must be “levelized”, or converted to a list of length max_levels such that entry [i] in that list is the parameter choice for level i.

The parameter choice in to_levelize can be a string, tuple or list. If it is a string or tuple, then that option is assumed to be the parameter setting at every level. If to_levelize is inititally a list, if the length of the list is less than max_levels, the last entry in the list defines that parameter for all subsequent levels.

>>> from pyamg.util.utils import levelize_smooth_or_improve_candidates
>>> improve_candidates = ['gauss_seidel', None]
>>> levelize_smooth_or_improve_candidates(improve_candidates, 4)
['gauss_seidel', None, None, None]
pyamg.util.levelize_strength_or_aggregation(to_levelize, max_levels, max_coarse)

Helper function to preprocess the strength and aggregation parameters passed to smoothed_aggregation_solver and rootnode_solver.

to_levelize : {string, tuple, list}
Parameter to preprocess, i.e., levelize and convert to a level-by-level list such that entry i specifies the parameter at level i
max_levels : int
Defines the maximum number of levels considered
max_coarse : int
Defines the maximum coarse grid size allowed
(max_levels, max_coarse, to_levelize) : tuple
New max_levels and max_coarse values and then the parameter list to_levelize, such that entry i specifies the parameter choice at level i. max_levels and max_coarse are returned, because they may be updated if strength or aggregation set a predefined coarsening and possibly change these values.

This routine is needed because the user will pass in a parameter option such as smooth=’jacobi’, or smooth=[‘jacobi’, None], and this option must be “levelized”, or converted to a list of length max_levels such that entry [i] in that list is the parameter choice for level i.

The parameter choice in to_levelize can be a string, tuple or list. If it is a string or tuple, then that option is assumed to be the parameter setting at every level. If to_levelize is inititally a list, if the length of the list is less than max_levels, the last entry in the list defines that parameter for all subsequent levels.

>>> from pyamg.util.utils import levelize_strength_or_aggregation
>>> strength = ['evolution', 'classical']
>>> levelize_strength_or_aggregation(strength, 4, 10)
(4, 10, ['evolution', 'classical', 'classical'])
pyamg.util.norm(x, pnorm='2')

2-norm of a vector

x : array_like
Vector of complex or real values
pnorm : string
‘2’ calculates the 2-norm ‘inf’ calculates the infinity-norm
n : float
2-norm of a vector
  • currently 1+ order of magnitude faster than scipy.linalg.norm(x), which calls sqrt(numpy.sum(real((conjugate(x)*x)),axis=0)) resulting in an extra copy
  • only handles the 2-norm and infinity-norm for vectors

scipy.linalg.norm : scipy general matrix or vector norm

pyamg.util.pinv_array(a, cond=None)
Calculate the Moore-Penrose pseudo inverse of each block of
the three dimensional array a.
a : {dense array}
Is of size (n, m, m)
cond : {float}
Used by *gelss to filter numerically zeros singular values. If None, a suitable value is chosen for you.

Nothing, a is modified in place so that a[k] holds the pseudoinverse of that block.

By using lapack wrappers, this can be much faster for large n, than directly calling pinv2

>>> import numpy
>>> from pyamg.util.linalg import pinv_array
>>> a = numpy.array([[[1.,2.],[1.,1.]], [[1.,1.],[3.,3.]]])
>>> ac = a.copy()
>>> # each block of a is inverted in-place
>>> pinv_array(a)
pyamg.util.print_table(table, title='', delim='|', centering='center', col_padding=2, header=True, headerchar='-')

Print a table from a list of lists representing the rows of a table

table : list
list of lists, e.g. a table with 3 columns and 2 rows could be [ [‘0,0’, ‘0,1’, ‘0,2’], [‘1,0’, ‘1,1’, ‘1,2’] ]
title : string
Printed centered above the table
delim : string
character to delimit columns
centering : {‘left’, ‘right’, ‘center’}
chooses justification for columns
col_padding : int
number of blank spaces to add to each column
header : {True, False}
Does the first entry of table contain column headers?
headerchar : {string}
character to separate column headers from rest of table

string representing table that’s ready to be printed

The string for the table will have correctly justified columns with extra padding added into each column entry to ensure columns align. The characters to delimit the columns can be user defined. This should be useful for printing convergence data from tests.

>>> from pyamg.util.utils import print_table
>>> table = [ ['cos(0)', 'cos(pi/2)', 'cos(pi)'], ['0.0', '1.0', '0.0'] ]
>>> table1 = print_table(table)                 # string to print
>>> table2 = print_table(table, delim='||')
>>> table3 = print_table(table, headerchar='*')
>>> table4 = print_table(table, col_padding=6, centering='left')
pyamg.util.profile_solver(ml, accel=None, **kwargs)

A quick solver to profile a particular multilevel object

ml : multilevel
Fully constructed multilevel object
accel : function pointer
Pointer to a valid Krylov solver (e.g. gmres, cg)
residuals : array
Array of residuals for each iteration

multilevel.psolve, multilevel.solve

>>> import numpy
>>> from scipy.sparse import spdiags, csr_matrix
>>> from scipy.sparse.linalg import cg
>>> from pyamg.classical import ruge_stuben_solver
>>> from pyamg.util.utils import profile_solver
>>> n=100
>>> e = numpy.ones((n,1)).ravel()
>>> data = [ -1*e, 2*e, -1*e ]
>>> A = csr_matrix(spdiags(data,[-1,0,1],n,n))
>>> b = A*numpy.ones(A.shape[0])
>>> ml = ruge_stuben_solver(A, max_coarse=10)
>>> res = profile_solver(ml,accel=cg)
pyamg.util.relaxation_as_linear_operator(method, A, b)

Create a linear operator that applies a relaxation method for the given right-hand-side

methods : {tuple or string}
Relaxation descriptor: Each tuple must be of the form (‘method’,’opts’) where ‘method’ is the name of a supported smoother, e.g., gauss_seidel, and ‘opts’ a dict of keyword arguments to the smoother, e.g., opts = {‘sweep’:symmetric}. If string, must be that of a supported smoother, e.g., gauss_seidel.

linear operator that applies the relaxation method to a vector for a fixed right-hand-side, b.

This method is primarily used to improve B during the aggregation setup phase. Here b = 0, and each relaxation call can improve the quality of B, especially near the boundaries.

>>> from pyamg.gallery import poisson
>>> from pyamg.util.utils import relaxation_as_linear_operator
>>> import numpy
>>> A = poisson((100,100), format='csr')           # matrix
>>> B = numpy.ones((A.shape[0],1))                 # Candidate vector
>>> b = numpy.zeros((A.shape[0]))                  # RHS
>>> relax = relaxation_as_linear_operator('gauss_seidel', A, b)
>>> B = relax*B
pyamg.util.residual_norm(A, x, b)

Compute ||b - A*x||

pyamg.util.scale_T(T, P_I, I_F)

Helper function that scales T with a right multiplication by a block diagonal inverse, so that T is the identity at C-node rows.

T : {bsr_matrix}
Tentative prolongator, with square blocks in the BSR data structure, and a non-overlapping block-diagonal structure
P_I : {bsr_matrix}
Interpolation operator that carries out only simple injection from the coarse grid to fine grid Cpts nodes
I_F : {bsr_matrix}
Identity operator on Fpts, i.e., the action of this matrix zeros out entries in a vector at all Cpts, leaving Fpts untouched
T : {bsr_matrix}
Tentative prolongator scaled to be identity at C-pt nodes
>>> from scipy.sparse import csr_matrix, bsr_matrix
>>> from scipy import matrix, array
>>> from pyamg.util.utils import scale_T
>>> T = matrix([[ 1.0,  0.,   0. ],
...             [ 0.5,  0.,   0. ],
...             [ 0. ,  1.,   0. ],
...             [ 0. ,  0.5,  0. ],
...             [ 0. ,  0.,   1. ],
...             [ 0. ,  0.,   0.25 ]])
>>> P_I = matrix([[ 0.,  0.,   0. ],
...               [ 1.,  0.,   0. ],
...               [ 0.,  1.,   0. ],
...               [ 0.,  0.,   0. ],
...               [ 0.,  0.,   0. ],
...               [ 0.,  0.,   1. ]])
>>> I_F = matrix([[ 1.,  0.,  0.,  0.,  0.,  0.],
...               [ 0.,  0.,  0.,  0.,  0.,  0.],
...               [ 0.,  0.,  0.,  0.,  0.,  0.],
...               [ 0.,  0.,  0.,  1.,  0.,  0.],
...               [ 0.,  0.,  0.,  0.,  1.,  0.],
...               [ 0.,  0.,  0.,  0.,  0.,  0.]])
>>> scale_T(bsr_matrix(T), bsr_matrix(P_I), bsr_matrix(I_F)).todense()
matrix([[ 2. ,  0. ,  0. ],
        [ 1. ,  0. ,  0. ],
        [ 0. ,  1. ,  0. ],
        [ 0. ,  0.5,  0. ],
        [ 0. ,  0. ,  4. ],
        [ 0. ,  0. ,  1. ]])

This routine is primarily used in pyamg.aggregation.smooth.energy_prolongation_smoother, where it is used to generate a suitable initial guess for the energy-minimization process, when root-node style SA is used. This function, scale_T, takes an existing tentative prolongator and ensures that it injects from the coarse-grid to fine-grid root-nodes.

When generating initial guesses for root-node style prolongation operators, this function is usually called after pyamg.uti.utils.filter_operator

This function assumes that the eventual coarse-grid nullspace vectors equal coarse-grid injection applied to the fine-grid nullspace vectors.

pyamg.util.symmetric_rescaling(A, copy=True)

Scale the matrix symmetrically:

A = D^{-1/2} A D^{-1/2}

where D=diag(A).

The left multiplication is accomplished through scale_rows and the right multiplication is done through scale columns.

A : sparse matrix
Sparse matrix with N rows
copy : {True,False}
  • If copy=True, then the matrix is copied to a new and different return matrix (e.g. B=symmetric_rescaling(A))
  • If copy=False, then the matrix is overwritten deeply (e.g. symmetric_rescaling(A,copy=False) overwrites A)
D_sqrt : array
Array of sqrt(diag(A))
D_sqrt_inv : array
Array of 1/sqrt(diag(A))
DAD : csr_matrix
Symmetrically scaled A
  • if A is not csr, it is converted to csr and sent to scale_rows
>>> import numpy
>>> from scipy.sparse import spdiags
>>> from pyamg.util.utils import symmetric_rescaling
>>> n=5
>>> e = numpy.ones((n,1)).ravel()
>>> data = [ -1*e, 2*e, -1*e ]
>>> A = spdiags(data,[-1,0,1],n,n).tocsr()
>>> Ds, Dsi, DAD = symmetric_rescaling(A)
>>> print DAD.todense()
[[ 1.  -0.5  0.   0.   0. ]
 [-0.5  1.  -0.5  0.   0. ]
 [ 0.  -0.5  1.  -0.5  0. ]
 [ 0.   0.  -0.5  1.  -0.5]
 [ 0.   0.   0.  -0.5  1. ]]
pyamg.util.symmetric_rescaling_sa(A, B, BH=None)

Scale the matrix symmetrically:

A = D^{-1/2} A D^{-1/2}

where D=diag(A). The left multiplication is accomplished through scale_rows and the right multiplication is done through scale columns.

The candidates B and BH are scaled accordingly:

B = D^{1/2} B
BH = D^{1/2} BH
A : {sparse matrix}
Sparse matrix with N rows
B : {array}
N x m array
BH : {None, array}
If A.symmetry == ‘nonsymmetric, then BH must be an N x m array. Otherwise, BH is ignored.

Appropriately scaled A, B and BH, i.e., A = D^{-1/2} A D^{-1/2}, B = D^{1/2} B, and BH = D^{1/2} BH

  • if A is not csr, it is converted to csr and sent to scale_rows
>>> import numpy
>>> from scipy.sparse import spdiags
>>> from pyamg.util.utils import symmetric_rescaling_sa
>>> n=5
>>> e = numpy.ones((n,1)).ravel()
>>> data = [ -1*e, 2*e, -1*e ]
>>> A = spdiags(data,[-1,0,1],n,n).tocsr()
>>> B = e.copy().reshape(-1,1)
>>> [DAD, DB, DBH] = symmetric_rescaling_sa(A,B,BH=None)
>>> print DAD.todense()
[[ 1.  -0.5  0.   0.   0. ]
 [-0.5  1.  -0.5  0.   0. ]
 [ 0.  -0.5  1.  -0.5  0. ]
 [ 0.   0.  -0.5  1.  -0.5]
 [ 0.   0.   0.  -0.5  1. ]]
>>> print DB
[[ 1.41421356]
 [ 1.41421356]
 [ 1.41421356]
 [ 1.41421356]
 [ 1.41421356]]
pyamg.util.to_type(upcast_type, varlist)

Loop over all elements of varlist and convert them to upcasttype

upcast_type : data type
e.g. complex, float64 or complex128
varlist : list
list may contain arrays, mat’s, sparse matrices, or scalars the elements may be float, int or complex

Returns upcast-ed varlist to upcast_type

Useful when harmonizing the types of variables, such as if A and b are complex, but x,y and z are not.

>>> import numpy
>>> from pyamg.util.utils import to_type
>>> from scipy.sparse.sputils import upcast
>>> x = numpy.ones((5,1))
>>> y = 2.0j*numpy.ones((5,1))
>>> varlist = to_type(upcast(x.dtype, y.dtype), [x, y])
pyamg.util.type_prep(upcast_type, varlist)

Loop over all elements of varlist and convert them to upcasttype The only difference with pyamg.util.utils.to_type(...), is that scalars are wrapped into (1,0) arrays. This is desirable when passing the numpy complex data type to C routines and complex scalars aren’t handled correctly

upcast_type : data type
e.g. complex, float64 or complex128
varlist : list
list may contain arrays, mat’s, sparse matrices, or scalars the elements may be float, int or complex

Returns upcast-ed varlist to upcast_type

Useful when harmonizing the types of variables, such as if A and b are complex, but x,y and z are not.

>>> import numpy
>>> from pyamg.util.utils import type_prep
>>> from scipy.sparse.sputils import upcast
>>> x = numpy.ones((5,1))
>>> y = 2.0j*numpy.ones((5,1))
>>> z = 2.3
>>> varlist = type_prep(upcast(x.dtype, y.dtype), [x, y, z])

BSR_utils Module

Utility Functions for reading and writing individual rows in BSR matrices

pyamg.util.BSR_utils.BSR_Get_Row(A, i)[source]

Return row i in BSR matrix A. Only nonzero entries are returned

A : bsr_matrix
Input matrix
i : int
Row number
z : array
Actual nonzero values for row i colindx Array of column indices for the nonzeros of row i
>>> from numpy import array
>>> from scipy.sparse import bsr_matrix
>>> from pyamg.util.BSR_utils import BSR_Get_Row
>>> indptr  = array([0,2,3,6])
>>> indices = array([0,2,2,0,1,2])
>>> data    = array([1,2,3,4,5,6]).repeat(4).reshape(6,2,2)
>>> B = bsr_matrix( (data,indices,indptr), shape=(6,6) )
>>> Brow = BSR_Get_Row(B,2)
>>> print Brow[1]
[4 5]
pyamg.util.BSR_utils.BSR_Row_WriteScalar(A, i, x)[source]

Write a scalar at each nonzero location in row i of BSR matrix A

A : bsr_matrix
Input matrix
i : int
Row number
x : float
Scalar to overwrite nonzeros of row i in A
A : bsr_matrix
All nonzeros in row i of A have been overwritten with x. If x is a vector, the first length(x) nonzeros in row i of A have been overwritten with entries from x
>>> from numpy import array
>>> from scipy.sparse import bsr_matrix
>>> from pyamg.util.BSR_utils import BSR_Row_WriteScalar
>>> indptr  = array([0,2,3,6])
>>> indices = array([0,2,2,0,1,2])
>>> data    = array([1,2,3,4,5,6]).repeat(4).reshape(6,2,2)
>>> B = bsr_matrix( (data,indices,indptr), shape=(6,6) )
>>> BSR_Row_WriteScalar(B,5,22)
pyamg.util.BSR_utils.BSR_Row_WriteVect(A, i, x)[source]

Overwrite the nonzeros in row i of BSR matrix A with the vector x. length(x) and nnz(A[i,:]) must be equivalent.

A : bsr_matrix
Matrix assumed to be in BSR format
i : int
Row number
x : array
Array of values to overwrite nonzeros in row i of A
A : bsr_matrix
The nonzeros in row i of A have been overwritten with entries from x. x must be same length as nonzeros of row i. This is guaranteed when this routine is used with vectors derived form Get_BSR_Row
>>> from numpy import array
>>> from scipy.sparse import bsr_matrix
>>> from pyamg.util.BSR_utils import BSR_Row_WriteVect
>>> indptr  = array([0,2,3,6])
>>> indices = array([0,2,2,0,1,2])
>>> data    = array([1,2,3,4,5,6]).repeat(4).reshape(6,2,2)
>>> B = bsr_matrix( (data,indices,indptr), shape=(6,6) )
>>> BSR_Row_WriteVect(B,5,array([11,22,33,44,55,66]))

info Module

linalg.py provides some linear algebra functionality not yet found in scipy

utils.py provides some utility functions for use with pyamg

BSR_utils.py provides utility functions for accessing and writing individual rows of BSR matrices

linalg Module

Linear Algebra Helper Routines

pyamg.util.linalg.approximate_spectral_radius(A, tol=0.01, maxiter=15, restart=5, symmetric=None, initial_guess=None, return_vector=False)[source]

Approximate the spectral radius of a matrix

A : {dense or sparse matrix}
E.g. csr_matrix, csc_matrix, ndarray, etc.
tol : {scalar}
Relative tolerance of approximation, i.e., the error divided by the approximate spectral radius is compared to tol.
maxiter : {integer}
Maximum number of iterations to perform
restart : {integer}
Number of restarted Arnoldi processes. For example, a value of 0 will run Arnoldi once, for maxiter iterations, and a value of 1 will restart Arnoldi once, using the maximal eigenvector from the first Arnoldi process as the initial guess.
symmetric : {boolean}
True - if A is symmetric
Lanczos iteration is used (more efficient)
False - if A is non-symmetric (default
Arnoldi iteration is used (less efficient)
initial_guess : {array|None}
If n x 1 array, then use as initial guess for Arnoldi/Lanczos. If None, then use a random initial guess.
return_vector : {boolean}
True - return an approximate dominant eigenvector, in addition to the
spectral radius.

False - Do not return the approximate dominant eigenvector

An approximation to the spectral radius of A, and if return_vector=True, then also return the approximate dominant eigenvector

The spectral radius is approximated by looking at the Ritz eigenvalues. Arnoldi iteration (or Lanczos) is used to project the matrix A onto a Krylov subspace: H = Q* A Q. The eigenvalues of H (i.e. the Ritz eigenvalues) should represent the eigenvalues of A in the sense that the minimum and maximum values are usually well matched (for the symmetric case it is true since the eigenvalues are real).

[1]Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, editors. “Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide”, SIAM, Philadelphia, 2000.
>>> from pyamg.util.linalg import approximate_spectral_radius
>>> import numpy
>>> from scipy.linalg import eigvals, norm
>>> A = numpy.array([[1.,0.],[0.,1.]])
>>> print approximate_spectral_radius(A,maxiter=3)
1.0
>>> print max([norm(x) for x in eigvals(A)])
1.0
pyamg.util.linalg.infinity_norm(A)[source]

Infinity norm of a matrix (maximum absolute row sum).

A : csr_matrix, csc_matrix, sparse, or numpy matrix
Sparse or dense matrix
n : float
Infinity norm of the matrix
  • This serves as an upper bound on spectral radius.
  • csr and csc avoid a deep copy
  • dense calls scipy.linalg.norm

scipy.linalg.norm : dense matrix norms

>>> import numpy
>>> from scipy.sparse import spdiags
>>> from pyamg.util.linalg import infinity_norm
>>> n=10
>>> e = numpy.ones((n,1)).ravel()
>>> data = [ -1*e, 2*e, -1*e ]
>>> A = spdiags(data,[-1,0,1],n,n)
>>> print infinity_norm(A)
4.0
pyamg.util.linalg.norm(x, pnorm='2')[source]

2-norm of a vector

x : array_like
Vector of complex or real values
pnorm : string
‘2’ calculates the 2-norm ‘inf’ calculates the infinity-norm
n : float
2-norm of a vector
  • currently 1+ order of magnitude faster than scipy.linalg.norm(x), which calls sqrt(numpy.sum(real((conjugate(x)*x)),axis=0)) resulting in an extra copy
  • only handles the 2-norm and infinity-norm for vectors

scipy.linalg.norm : scipy general matrix or vector norm

pyamg.util.linalg.residual_norm(A, x, b)[source]

Compute ||b - A*x||

pyamg.util.linalg.condest(A, tol=0.1, maxiter=25, symmetric=False)[source]

Estimates the condition number of A

A : {dense or sparse matrix}
e.g. array, matrix, csr_matrix, ...
tol : {float}
Approximation tolerance, currently not used
maxiter: {int}
Max number of Arnoldi/Lanczos iterations
symmetric : {bool}
If symmetric use the far more efficient Lanczos algorithm, Else use Arnoldi

Estimate of cond(A) with |lambda_max| / |lambda_min| through the use of Arnoldi or Lanczos iterations, depending on the symmetric flag

The condition number measures how large of a change in the the problems solution is caused by a change in problem’s input. Large condition numbers indicate that small perturbations and numerical errors are magnified greatly when solving the system.

>>> import numpy
>>> from pyamg.util.linalg import condest
>>> c = condest(numpy.array([[1.,0.],[0.,2.]]))
>>> print c
2.0
pyamg.util.linalg.cond(A)[source]

Returns condition number of A

A : {dense or sparse matrix}
e.g. array, matrix, csr_matrix, ...

2-norm condition number through use of the SVD Use for small to moderate sized dense matrices. For large sparse matrices, use condest.

The condition number measures how large of a change in the problems solution is caused by a change in problem’s input. Large condition numbers indicate that small perturbations and numerical errors are magnified greatly when solving the system.

>>> import numpy
>>> from pyamg.util.linalg import condest
>>> c = condest(numpy.array([[1.0,0.],[0.,2.0]]))
>>> print c
2.0
pyamg.util.linalg.ishermitian(A, fast_check=True, tol=1e-06, verbose=False)[source]

Returns True if A is Hermitian to within tol

A : {dense or sparse matrix}
e.g. array, matrix, csr_matrix, ...
fast_check : {bool}
If True, use the heuristic < Ax, y> = < x, Ay> for random vectors x and y to check for conjugate symmetry. If False, compute A - A.H.
tol : {float}
Symmetry tolerance
verbose: {bool}
prints max( |A - A.H| ) if nonhermitian and fast_check=False abs( <Ax, y> - <x, Ay> ) if nonhermitian and fast_check=False

True if hermitian False if nonhermitian

This function applies a simple test of conjugate symmetry

>>> import numpy
>>> from pyamg.util.linalg import ishermitian
>>> ishermitian(numpy.array([[1,2],[1,1]]))
False
>>> from pyamg.gallery import poisson
>>> ishermitian(poisson((10,10)))
True
pyamg.util.linalg.pinv_array(a, cond=None)[source]
Calculate the Moore-Penrose pseudo inverse of each block of
the three dimensional array a.
a : {dense array}
Is of size (n, m, m)
cond : {float}
Used by *gelss to filter numerically zeros singular values. If None, a suitable value is chosen for you.

Nothing, a is modified in place so that a[k] holds the pseudoinverse of that block.

By using lapack wrappers, this can be much faster for large n, than directly calling pinv2

>>> import numpy
>>> from pyamg.util.linalg import pinv_array
>>> a = numpy.array([[[1.,2.],[1.,1.]], [[1.,1.],[3.,3.]]])
>>> ac = a.copy()
>>> # each block of a is inverted in-place
>>> pinv_array(a)

setup Module

pyamg.util.setup.configuration(parent_package='', top_path=None)[source]

utils Module

General utility functions for pyamg

pyamg.util.utils.blocksize(A)[source]
pyamg.util.utils.diag_sparse(A)[source]
If A is a sparse matrix (e.g. csr_matrix or csc_matrix)
  • return the diagonal of A as an array
Otherwise
  • return a csr_matrix with A on the diagonal
A : sparse matrix or 1d array
General sparse matrix or array of diagonal entries
B : array or sparse matrix
Diagonal sparse is returned as csr if A is dense otherwise return an array of the diagonal
>>> import numpy
>>> from pyamg.util.utils import diag_sparse
>>> d = 2.0*numpy.ones((3,)).ravel()
>>> print diag_sparse(d).todense()
[[ 2.  0.  0.]
 [ 0.  2.  0.]
 [ 0.  0.  2.]]
pyamg.util.utils.profile_solver(ml, accel=None, **kwargs)[source]

A quick solver to profile a particular multilevel object

ml : multilevel
Fully constructed multilevel object
accel : function pointer
Pointer to a valid Krylov solver (e.g. gmres, cg)
residuals : array
Array of residuals for each iteration

multilevel.psolve, multilevel.solve

>>> import numpy
>>> from scipy.sparse import spdiags, csr_matrix
>>> from scipy.sparse.linalg import cg
>>> from pyamg.classical import ruge_stuben_solver
>>> from pyamg.util.utils import profile_solver
>>> n=100
>>> e = numpy.ones((n,1)).ravel()
>>> data = [ -1*e, 2*e, -1*e ]
>>> A = csr_matrix(spdiags(data,[-1,0,1],n,n))
>>> b = A*numpy.ones(A.shape[0])
>>> ml = ruge_stuben_solver(A, max_coarse=10)
>>> res = profile_solver(ml,accel=cg)
pyamg.util.utils.to_type(upcast_type, varlist)[source]

Loop over all elements of varlist and convert them to upcasttype

upcast_type : data type
e.g. complex, float64 or complex128
varlist : list
list may contain arrays, mat’s, sparse matrices, or scalars the elements may be float, int or complex

Returns upcast-ed varlist to upcast_type

Useful when harmonizing the types of variables, such as if A and b are complex, but x,y and z are not.

>>> import numpy
>>> from pyamg.util.utils import to_type
>>> from scipy.sparse.sputils import upcast
>>> x = numpy.ones((5,1))
>>> y = 2.0j*numpy.ones((5,1))
>>> varlist = to_type(upcast(x.dtype, y.dtype), [x, y])
pyamg.util.utils.type_prep(upcast_type, varlist)[source]

Loop over all elements of varlist and convert them to upcasttype The only difference with pyamg.util.utils.to_type(...), is that scalars are wrapped into (1,0) arrays. This is desirable when passing the numpy complex data type to C routines and complex scalars aren’t handled correctly

upcast_type : data type
e.g. complex, float64 or complex128
varlist : list
list may contain arrays, mat’s, sparse matrices, or scalars the elements may be float, int or complex

Returns upcast-ed varlist to upcast_type

Useful when harmonizing the types of variables, such as if A and b are complex, but x,y and z are not.

>>> import numpy
>>> from pyamg.util.utils import type_prep
>>> from scipy.sparse.sputils import upcast
>>> x = numpy.ones((5,1))
>>> y = 2.0j*numpy.ones((5,1))
>>> z = 2.3
>>> varlist = type_prep(upcast(x.dtype, y.dtype), [x, y, z])
pyamg.util.utils.get_diagonal(A, norm_eq=False, inv=False)[source]
Return the diagonal or inverse of diagonal for
A, (A.H A) or (A A.H)
A : {dense or sparse matrix}
e.g. array, matrix, csr_matrix, ...
norm_eq : {0, 1, 2}
0 ==> D = diag(A) 1 ==> D = diag(A.H A) 2 ==> D = diag(A A.H)
inv : {True, False}
If True, D = 1.0/D

diagonal, D, of appropriate system

This function is especially useful for its fast methods of obtaining diag(A A.H) and diag(A.H A). Dinv is zero wherever D is zero

>>> from pyamg.util.utils import get_diagonal
>>> from pyamg.gallery import poisson
>>> A = poisson( (5,), format='csr' )
>>> D = get_diagonal(A)
>>> print D
[ 2.  2.  2.  2.  2.]
>>> D = get_diagonal(A, norm_eq=1, inv=True)
>>> print D
[ 0.2         0.16666667  0.16666667  0.16666667  0.2       ]
pyamg.util.utils.UnAmal(A, RowsPerBlock, ColsPerBlock)[source]

Unamalgamate a CSR A with blocks of 1’s. This operation is equivalent to replacing each entry of A with ones(RowsPerBlock, ColsPerBlock), i.e., this is equivalent to setting all of A’s nonzeros to 1 and then doing a Kronecker product between A and ones(RowsPerBlock, ColsPerBlock).

A : csr_matrix
Amalgamted matrix
RowsPerBlock : int
Give A blocks of size (RowsPerBlock, ColsPerBlock)
ColsPerBlock : int
Give A blocks of size (RowsPerBlock, ColsPerBlock)
A : bsr_matrix
Returns A.data[:] = 1, followed by a Kronecker product of A and ones(RowsPerBlock, ColsPerBlock)
>>> from numpy import array
>>> from scipy.sparse import csr_matrix
>>> from pyamg.util.utils import UnAmal
>>> row = array([0,0,1,2,2,2])
>>> col = array([0,2,2,0,1,2])
>>> data = array([1,2,3,4,5,6])
>>> A = csr_matrix( (data,(row,col)), shape=(3,3) )
>>> A.todense()
matrix([[1, 0, 2],
        [0, 0, 3],
        [4, 5, 6]])
>>> UnAmal(A,2,2).todense()
matrix([[ 1.,  1.,  0.,  0.,  1.,  1.],
        [ 1.,  1.,  0.,  0.,  1.,  1.],
        [ 0.,  0.,  0.,  0.,  1.,  1.],
        [ 0.,  0.,  0.,  0.,  1.,  1.],
        [ 1.,  1.,  1.,  1.,  1.,  1.],
        [ 1.,  1.,  1.,  1.,  1.,  1.]])
pyamg.util.utils.Coord2RBM(numNodes, numPDEs, x, y, z)[source]

Convert 2D or 3D coordinates into Rigid body modes for use as near nullspace modes in elasticity AMG solvers

numNodes : int
Number of nodes
numPDEs :
Number of dofs per node
x,y,z : array_like
Coordinate vectors
rbm : matrix
A matrix of size (numNodes*numPDEs) x (1 | 6) containing the 6 rigid body modes
>>> import numpy
>>> from pyamg.util.utils import Coord2RBM
>>> a = numpy.array([0,1,2])
>>> Coord2RBM(3,6,a,a,a)
matrix([[ 1.,  0.,  0.,  0.,  0., -0.],
        [ 0.,  1.,  0., -0.,  0.,  0.],
        [ 0.,  0.,  1.,  0., -0.,  0.],
        [ 0.,  0.,  0.,  1.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  1.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  1.],
        [ 1.,  0.,  0.,  0.,  1., -1.],
        [ 0.,  1.,  0., -1.,  0.,  1.],
        [ 0.,  0.,  1.,  1., -1.,  0.],
        [ 0.,  0.,  0.,  1.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  1.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  1.],
        [ 1.,  0.,  0.,  0.,  2., -2.],
        [ 0.,  1.,  0., -2.,  0.,  2.],
        [ 0.,  0.,  1.,  2., -2.,  0.],
        [ 0.,  0.,  0.,  1.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  1.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  1.]])
pyamg.util.utils.hierarchy_spectrum(mg, filter=True, plot=False)[source]

Examine a multilevel hierarchy’s spectrum

mg { pyamg multilevel hierarchy }
e.g. generated with smoothed_aggregation_solver(...) or ruge_stuben_solver(...)
  1. table to standard out detailing the spectrum of each level in mg
  2. if plot==True, a sequence of plots in the complex plane of the spectrum at each level

This can be useful for troubleshooting and when examining how your problem’s nature changes from level to level

>>> from pyamg import smoothed_aggregation_solver
>>> from pyamg.gallery import poisson
>>> from pyamg.util.utils import hierarchy_spectrum
>>> A = poisson( (1,), format='csr' )
>>> ml = smoothed_aggregation_solver(A)
>>> hierarchy_spectrum(ml)

 Level min(re(eig)) max(re(eig)) num re(eig) < 0 num re(eig) > 0 cond_2(A)
---------------------------------------------------------------------------
   0      2.000        2.000            0               1         1.00e+00


 Level min(im(eig)) max(im(eig)) num im(eig) < 0 num im(eig) > 0 cond_2(A)
---------------------------------------------------------------------------
   0      0.000        0.000            0               0         1.00e+00
pyamg.util.utils.print_table(table, title='', delim='|', centering='center', col_padding=2, header=True, headerchar='-')[source]

Print a table from a list of lists representing the rows of a table

table : list
list of lists, e.g. a table with 3 columns and 2 rows could be [ [‘0,0’, ‘0,1’, ‘0,2’], [‘1,0’, ‘1,1’, ‘1,2’] ]
title : string
Printed centered above the table
delim : string
character to delimit columns
centering : {‘left’, ‘right’, ‘center’}
chooses justification for columns
col_padding : int
number of blank spaces to add to each column
header : {True, False}
Does the first entry of table contain column headers?
headerchar : {string}
character to separate column headers from rest of table

string representing table that’s ready to be printed

The string for the table will have correctly justified columns with extra padding added into each column entry to ensure columns align. The characters to delimit the columns can be user defined. This should be useful for printing convergence data from tests.

>>> from pyamg.util.utils import print_table
>>> table = [ ['cos(0)', 'cos(pi/2)', 'cos(pi)'], ['0.0', '1.0', '0.0'] ]
>>> table1 = print_table(table)                 # string to print
>>> table2 = print_table(table, delim='||')
>>> table3 = print_table(table, headerchar='*')
>>> table4 = print_table(table, col_padding=6, centering='left')
pyamg.util.utils.get_block_diag(A, blocksize, inv_flag=True)[source]

Return the block diagonal of A, in array form

A : csr_matrix
assumed to be square
blocksize : int
square block size for the diagonal
inv_flag : bool
if True, return the inverse of the block diagonal
block_diag : array
block diagonal of A in array form, array size is (A.shape[0]/blocksize, blocksize, blocksize)
>>> from scipy import arange
>>> from scipy.sparse import csr_matrix
>>> from pyamg.util import get_block_diag
>>> A = csr_matrix(arange(36).reshape(6,6))
>>> block_diag_inv = get_block_diag(A, blocksize=2, inv_flag=False)
>>> print block_diag_inv
[[[  0.   1.]
  [  6.   7.]]

 [[ 14.  15.]
  [ 20.  21.]]

 [[ 28.  29.]
  [ 34.  35.]]]
>>> block_diag_inv = get_block_diag(A, blocksize=2, inv_flag=True)
pyamg.util.utils.amalgamate(A, blocksize)[source]

Amalgamate matrix A

A : csr_matrix
Matrix to amalgamate
blocksize : int
blocksize to use while amalgamating
A_amal : csr_matrix
Amalgamated matrix A, first, convert A to BSR with square blocksize and then return a CSR matrix of ones using the resulting BSR indptr and indices

inverse operation of UnAmal for square matrices

>>> from numpy import array
>>> from scipy.sparse import csr_matrix
>>> from pyamg.util.utils import amalgamate
>>> row = array([0,0,1])
>>> col = array([0,2,1])
>>> data = array([1,2,3])
>>> A = csr_matrix( (data,(row,col)), shape=(4,4) )
>>> A.todense()
matrix([[1, 0, 2, 0],
        [0, 3, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]])
>>> amalgamate(A,2).todense()
matrix([[ 1.,  1.],
        [ 0.,  0.]])
pyamg.util.utils.symmetric_rescaling(A, copy=True)[source]

Scale the matrix symmetrically:

A = D^{-1/2} A D^{-1/2}

where D=diag(A).

The left multiplication is accomplished through scale_rows and the right multiplication is done through scale columns.

A : sparse matrix
Sparse matrix with N rows
copy : {True,False}
  • If copy=True, then the matrix is copied to a new and different return matrix (e.g. B=symmetric_rescaling(A))
  • If copy=False, then the matrix is overwritten deeply (e.g. symmetric_rescaling(A,copy=False) overwrites A)
D_sqrt : array
Array of sqrt(diag(A))
D_sqrt_inv : array
Array of 1/sqrt(diag(A))
DAD : csr_matrix
Symmetrically scaled A
  • if A is not csr, it is converted to csr and sent to scale_rows
>>> import numpy
>>> from scipy.sparse import spdiags
>>> from pyamg.util.utils import symmetric_rescaling
>>> n=5
>>> e = numpy.ones((n,1)).ravel()
>>> data = [ -1*e, 2*e, -1*e ]
>>> A = spdiags(data,[-1,0,1],n,n).tocsr()
>>> Ds, Dsi, DAD = symmetric_rescaling(A)
>>> print DAD.todense()
[[ 1.  -0.5  0.   0.   0. ]
 [-0.5  1.  -0.5  0.   0. ]
 [ 0.  -0.5  1.  -0.5  0. ]
 [ 0.   0.  -0.5  1.  -0.5]
 [ 0.   0.   0.  -0.5  1. ]]
pyamg.util.utils.symmetric_rescaling_sa(A, B, BH=None)[source]

Scale the matrix symmetrically:

A = D^{-1/2} A D^{-1/2}

where D=diag(A). The left multiplication is accomplished through scale_rows and the right multiplication is done through scale columns.

The candidates B and BH are scaled accordingly:

B = D^{1/2} B
BH = D^{1/2} BH
A : {sparse matrix}
Sparse matrix with N rows
B : {array}
N x m array
BH : {None, array}
If A.symmetry == ‘nonsymmetric, then BH must be an N x m array. Otherwise, BH is ignored.

Appropriately scaled A, B and BH, i.e., A = D^{-1/2} A D^{-1/2}, B = D^{1/2} B, and BH = D^{1/2} BH

  • if A is not csr, it is converted to csr and sent to scale_rows
>>> import numpy
>>> from scipy.sparse import spdiags
>>> from pyamg.util.utils import symmetric_rescaling_sa
>>> n=5
>>> e = numpy.ones((n,1)).ravel()
>>> data = [ -1*e, 2*e, -1*e ]
>>> A = spdiags(data,[-1,0,1],n,n).tocsr()
>>> B = e.copy().reshape(-1,1)
>>> [DAD, DB, DBH] = symmetric_rescaling_sa(A,B,BH=None)
>>> print DAD.todense()
[[ 1.  -0.5  0.   0.   0. ]
 [-0.5  1.  -0.5  0.   0. ]
 [ 0.  -0.5  1.  -0.5  0. ]
 [ 0.   0.  -0.5  1.  -0.5]
 [ 0.   0.   0.  -0.5  1. ]]
>>> print DB
[[ 1.41421356]
 [ 1.41421356]
 [ 1.41421356]
 [ 1.41421356]
 [ 1.41421356]]
pyamg.util.utils.relaxation_as_linear_operator(method, A, b)[source]

Create a linear operator that applies a relaxation method for the given right-hand-side

methods : {tuple or string}
Relaxation descriptor: Each tuple must be of the form (‘method’,’opts’) where ‘method’ is the name of a supported smoother, e.g., gauss_seidel, and ‘opts’ a dict of keyword arguments to the smoother, e.g., opts = {‘sweep’:symmetric}. If string, must be that of a supported smoother, e.g., gauss_seidel.

linear operator that applies the relaxation method to a vector for a fixed right-hand-side, b.

This method is primarily used to improve B during the aggregation setup phase. Here b = 0, and each relaxation call can improve the quality of B, especially near the boundaries.

>>> from pyamg.gallery import poisson
>>> from pyamg.util.utils import relaxation_as_linear_operator
>>> import numpy
>>> A = poisson((100,100), format='csr')           # matrix
>>> B = numpy.ones((A.shape[0],1))                 # Candidate vector
>>> b = numpy.zeros((A.shape[0]))                  # RHS
>>> relax = relaxation_as_linear_operator('gauss_seidel', A, b)
>>> B = relax*B
pyamg.util.utils.filter_operator(A, C, B, Bf, BtBinv=None)[source]

Filter the matrix A according to the matrix graph of C, while ensuring that the new, filtered A satisfies: A_new*B = Bf.

A : {csr_matrix, bsr_matrix}
n x m matrix to filter
C : {csr_matrix, bsr_matrix}
n x m matrix representing the couplings in A to keep
B : {array}
m x k array of near nullspace vectors
Bf : {array}
n x k array of near nullspace vectors to place in span(A)
BtBinv : {None, array}
3 dimensional array such that, BtBinv[i] = pinv(B_i.H Bi), and B_i is B restricted to the neighborhood (with respect to the matrix graph of C) of dof of i. If None is passed in, this array is computed internally.

A : sparse matrix updated such that sparsity structure of A now matches that of C, and that the relationship A*B = Bf holds.

This procedure allows for certain important modes (i.e., Bf) to be placed in span(A) by way of row-wise l2-projections that enforce the relationship A*B = Bf. This is useful for maintaining certain modes (e.g., the constant) in the span of prolongation.

>>> from numpy import ones, array
>>> from scipy.sparse import csr_matrix
>>> from pyamg.util.utils import filter_operator
>>> A = array([ [1.,1,1],[1,1,1],[0,1,0],[0,1,0],[0,0,1],[0,0,1]])
>>> C = array([ [1.,1,0],[1,1,0],[0,1,0],[0,1,0],[0,0,1],[0,0,1]])
>>> B = ones((3,1))
>>> Bf = ones((6,1))
>>> filter_operator(csr_matrix(A), csr_matrix(C), B, Bf).todense()
matrix([[ 0.5,  0.5,  0. ],
        [ 0.5,  0.5,  0. ],
        [ 0. ,  1. ,  0. ],
        [ 0. ,  1. ,  0. ],
        [ 0. ,  0. ,  1. ],
        [ 0. ,  0. ,  1. ]])

This routine is primarily used in pyamg.aggregation.smooth.energy_prolongation_smoother, where it is used to generate a suitable initial guess for the energy-minimization process, when root-node style SA is used. Essentially, the tentative prolongator, T, is processed by this routine to produce fine-grid nullspace vectors when multiplying coarse-grid nullspace vectors, i.e., T*B = Bf. This is possible for any arbitrary vectors B and Bf, so long as the sparsity structure of T is rich enough.

When generating initial guesses for root-node style prolongation operators, this function is usually called before pyamg.uti.utils.scale_T

pyamg.util.utils.scale_T(T, P_I, I_F)[source]

Helper function that scales T with a right multiplication by a block diagonal inverse, so that T is the identity at C-node rows.

T : {bsr_matrix}
Tentative prolongator, with square blocks in the BSR data structure, and a non-overlapping block-diagonal structure
P_I : {bsr_matrix}
Interpolation operator that carries out only simple injection from the coarse grid to fine grid Cpts nodes
I_F : {bsr_matrix}
Identity operator on Fpts, i.e., the action of this matrix zeros out entries in a vector at all Cpts, leaving Fpts untouched
T : {bsr_matrix}
Tentative prolongator scaled to be identity at C-pt nodes
>>> from scipy.sparse import csr_matrix, bsr_matrix
>>> from scipy import matrix, array
>>> from pyamg.util.utils import scale_T
>>> T = matrix([[ 1.0,  0.,   0. ],
...             [ 0.5,  0.,   0. ],
...             [ 0. ,  1.,   0. ],
...             [ 0. ,  0.5,  0. ],
...             [ 0. ,  0.,   1. ],
...             [ 0. ,  0.,   0.25 ]])
>>> P_I = matrix([[ 0.,  0.,   0. ],
...               [ 1.,  0.,   0. ],
...               [ 0.,  1.,   0. ],
...               [ 0.,  0.,   0. ],
...               [ 0.,  0.,   0. ],
...               [ 0.,  0.,   1. ]])
>>> I_F = matrix([[ 1.,  0.,  0.,  0.,  0.,  0.],
...               [ 0.,  0.,  0.,  0.,  0.,  0.],
...               [ 0.,  0.,  0.,  0.,  0.,  0.],
...               [ 0.,  0.,  0.,  1.,  0.,  0.],
...               [ 0.,  0.,  0.,  0.,  1.,  0.],
...               [ 0.,  0.,  0.,  0.,  0.,  0.]])
>>> scale_T(bsr_matrix(T), bsr_matrix(P_I), bsr_matrix(I_F)).todense()
matrix([[ 2. ,  0. ,  0. ],
        [ 1. ,  0. ,  0. ],
        [ 0. ,  1. ,  0. ],
        [ 0. ,  0.5,  0. ],
        [ 0. ,  0. ,  4. ],
        [ 0. ,  0. ,  1. ]])

This routine is primarily used in pyamg.aggregation.smooth.energy_prolongation_smoother, where it is used to generate a suitable initial guess for the energy-minimization process, when root-node style SA is used. This function, scale_T, takes an existing tentative prolongator and ensures that it injects from the coarse-grid to fine-grid root-nodes.

When generating initial guesses for root-node style prolongation operators, this function is usually called after pyamg.uti.utils.filter_operator

This function assumes that the eventual coarse-grid nullspace vectors equal coarse-grid injection applied to the fine-grid nullspace vectors.

pyamg.util.utils.get_Cpt_params(A, Cnodes, AggOp, T)[source]
Helper function that returns a dictionary of sparse matrices and arrays
which allow us to easily operate on Cpts and Fpts separately.
A : {csr_matrix, bsr_matrix}
Operator
Cnodes : {array}
Array of all root node indices. This is an array of nodal indices, not degree-of-freedom indices. If the blocksize of T is 1, then nodal indices and degree-of-freedom indices coincide.
AggOp : {csr_matrix}
Aggregation operator corresponding to A
T : {bsr_matrix}
Tentative prolongator based on AggOp

Dictionary containing these parameters:

P_I : {bsr_matrix}
Interpolation operator that carries out only simple injection from the coarse grid to fine grid Cpts nodes
I_F : {bsr_matrix}
Identity operator on Fpts, i.e., the action of this matrix zeros out entries in a vector at all Cpts, leaving Fpts untouched
I_C : {bsr_matrix}
Identity operator on Cpts nodes, i.e., the action of this matrix zeros out entries in a vector at all Fpts, leaving Cpts untouched
Cpts : {array}
An array of all root node dofs, corresponding to the F/C splitting
Fpts : {array}
An array of all non root node dofs, corresponding to the F/C splitting
>>> from numpy import array
>>> from pyamg.util.utils import get_Cpt_params
>>> from pyamg.gallery import poisson
>>> from scipy.sparse import csr_matrix, bsr_matrix
>>> A = poisson((10,), format='csr')
>>> Cpts = array([3, 7])
>>> AggOp = ([[ 1., 0.], [ 1., 0.],
...           [ 1., 0.], [ 1., 0.],
...           [ 1., 0.], [ 0., 1.],
...           [ 0., 1.], [ 0., 1.],
...           [ 0., 1.], [ 0., 1.]])
>>> AggOp = csr_matrix(AggOp)
>>> T = AggOp.copy().tobsr()
>>> params = get_Cpt_params(A, Cpts, AggOp, T)
>>> params['P_I'].todense()
matrix([[ 0.,  0.],
        [ 0.,  0.],
        [ 0.,  0.],
        [ 1.,  0.],
        [ 0.,  0.],
        [ 0.,  0.],
        [ 0.,  0.],
        [ 0.,  1.],
        [ 0.,  0.],
        [ 0.,  0.]])

The principal calling routine is aggregation.smooth.energy_prolongation_smoother, which uses the Cpt_param dictionary for root-node style prolongation smoothing

pyamg.util.utils.compute_BtBinv(B, C)[source]
Helper function that creates inv(B_i.T B_i) for each block row i in C,
where B_i is B restricted to the sparsity pattern of block row i.
B : {array}
(M,k) array, typically near-nullspace modes for coarse grid, i.e., B_c.
C : {csr_matrix, bsr_matrix}
Sparse NxM matrix, whose sparsity structure (i.e., matrix graph) is used to determine BtBinv.
BtBinv : {array}
BtBinv[i] = inv(B_i.T B_i), where B_i is B restricted to the nonzero pattern of block row i in C.
>>> from numpy import array
>>> from scipy.sparse import bsr_matrix
>>> from pyamg.util.utils import compute_BtBinv
>>> T = array([[ 1.,  0.],
...            [ 1.,  0.],
...            [ 0.,  .5],
...            [ 0.,  .25]])
>>> T = bsr_matrix(T)
>>> B = array([[1.],[2.]])
>>> compute_BtBinv(B, T)
array([[[ 1.  ]],

       [[ 1.  ]],

       [[ 0.25]],

       [[ 0.25]]])

The principal calling routines are aggregation.smooth.energy_prolongation_smoother, and util.utils.filter_operator.

BtBinv is used in the prolongation smoothing process that incorporates B into the span of prolongation with row-wise projection operators. It is these projection operators that BtBinv is part of.

pyamg.util.utils.eliminate_diag_dom_nodes(A, C, theta=1.02)[source]

Helper function that eliminates diagonally dominant rows and cols from A in the separate matrix C. This is useful because it eliminates nodes in C which we don’t want coarsened. These eliminated nodes in C just become the rows and columns of the identity.

A : {csr_matrix, bsr_matrix}
Sparse NxN matrix
C : {csr_matrix}
Sparse MxM matrix, where M is the number of nodes in A. M=N if A is CSR or is BSR with blocksize 1. Otherwise M = N/blocksize.
theta : {float}
determines diagonal dominance threshhold
C : {csr_matrix}
C updated such that the rows and columns corresponding to diagonally dominant rows in A have been eliminated and replaced with rows and columns of the identity.
Diagonal dominance is defined as
|| (e_i, A) - a_ii ||_1 < theta a_ii

that is, the 1-norm of the off diagonal elements in row i must be less than theta times the diagonal element.

>>> from pyamg.gallery import poisson
>>> from pyamg.util.utils import eliminate_diag_dom_nodes
>>> A = poisson( (4,), format='csr' )
>>> C = eliminate_diag_dom_nodes(A, A.copy(), 1.1)
>>> C.todense()
matrix([[ 1.,  0.,  0.,  0.],
        [ 0.,  2., -1.,  0.],
        [ 0., -1.,  2.,  0.],
        [ 0.,  0.,  0.,  1.]])
pyamg.util.utils.levelize_strength_or_aggregation(to_levelize, max_levels, max_coarse)[source]

Helper function to preprocess the strength and aggregation parameters passed to smoothed_aggregation_solver and rootnode_solver.

to_levelize : {string, tuple, list}
Parameter to preprocess, i.e., levelize and convert to a level-by-level list such that entry i specifies the parameter at level i
max_levels : int
Defines the maximum number of levels considered
max_coarse : int
Defines the maximum coarse grid size allowed
(max_levels, max_coarse, to_levelize) : tuple
New max_levels and max_coarse values and then the parameter list to_levelize, such that entry i specifies the parameter choice at level i. max_levels and max_coarse are returned, because they may be updated if strength or aggregation set a predefined coarsening and possibly change these values.

This routine is needed because the user will pass in a parameter option such as smooth=’jacobi’, or smooth=[‘jacobi’, None], and this option must be “levelized”, or converted to a list of length max_levels such that entry [i] in that list is the parameter choice for level i.

The parameter choice in to_levelize can be a string, tuple or list. If it is a string or tuple, then that option is assumed to be the parameter setting at every level. If to_levelize is inititally a list, if the length of the list is less than max_levels, the last entry in the list defines that parameter for all subsequent levels.

>>> from pyamg.util.utils import levelize_strength_or_aggregation
>>> strength = ['evolution', 'classical']
>>> levelize_strength_or_aggregation(strength, 4, 10)
(4, 10, ['evolution', 'classical', 'classical'])
pyamg.util.utils.levelize_smooth_or_improve_candidates(to_levelize, max_levels)[source]

Helper function to preprocess the smooth and improve_candidates parameters passed to smoothed_aggregation_solver and rootnode_solver.

to_levelize : {string, tuple, list}
Parameter to preprocess, i.e., levelize and convert to a level-by-level list such that entry i specifies the parameter at level i
max_levels : int
Defines the maximum number of levels considered
to_levelize : list
The parameter list such that entry i specifies the parameter choice at level i.

This routine is needed because the user will pass in a parameter option such as smooth=’jacobi’, or smooth=[‘jacobi’, None], and this option must be “levelized”, or converted to a list of length max_levels such that entry [i] in that list is the parameter choice for level i.

The parameter choice in to_levelize can be a string, tuple or list. If it is a string or tuple, then that option is assumed to be the parameter setting at every level. If to_levelize is inititally a list, if the length of the list is less than max_levels, the last entry in the list defines that parameter for all subsequent levels.

>>> from pyamg.util.utils import levelize_smooth_or_improve_candidates
>>> improve_candidates = ['gauss_seidel', None]
>>> levelize_smooth_or_improve_candidates(improve_candidates, 4)
['gauss_seidel', None, None, None]