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(...)
- table to standard out detailing the spectrum of each level in mg
- 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.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)
utils Module¶
General utility functions for pyamg
- 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(...)
- table to standard out detailing the spectrum of each level in mg
- 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]