This page contains the Util Package documentation.
Utility Functions for reading and writing individual rows in BSR matrices
Return row i in BSR matrix A. Only nonzero entries are returned
Parameters : | A : bsr_matrix
i : int
|
---|---|
Returns : | z : array
|
Examples
>>> 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]
Write a scalar at each nonzero location in row i of BSR matrix A
Parameters : | A : bsr_matrix
i : int
x : float
|
---|---|
Returns : | A : bsr_matrix
|
Examples
>>> 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)
Overwrite the nonzeros in row i of BSR matrix A with the vector x. length(x) and nnz(A[i,:]) must be equivalent.
Parameters : | A : bsr_matrix
i : int
x : array
|
---|---|
Returns : | A : bsr_matrix
|
Examples
>>> 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]))
Linear Algebra Helper Routines
Approximate the spectral radius of a matrix
Parameters : | A : {dense or sparse matrix}
tol : {scalar}
maxiter : {integer}
restart : {integer}
symmetric : {boolean}
initial_guess : {array|None}
|
---|---|
Returns : | An approximation to the spectral radius of A : |
Notes
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).
References
[R34] | 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. |
Examples
>>> 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
Infinity norm of a matrix (maximum absolute row sum).
Parameters : | A : csr_matrix, csc_matrix, sparse, or numpy matrix
|
---|---|
Returns : | n : float
|
See also
Notes
Examples
>>> 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
2-norm of a vector
Parameters : | x : array_like
pnorm : string
|
---|---|
Returns : | n : float
|
See also
Notes
Estimates the condition number of A
Parameters : | A : {dense or sparse matrix}
tol : {float}
maxiter: {int} :
symmetric : {bool}
|
---|---|
Returns : | Estimate of cond(A) with |lambda_max| / |lambda_min| : through the use of Arnoldi or Lanczos iterations, depending on : the symmetric flag : |
Notes
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.
Examples
>>> import numpy
>>> from pyamg.util.linalg import condest
>>> c = condest(numpy.array([[1.,0.],[0.,2.]]))
>>> print c
2.0
Returns condition number of A
Parameters : | A : {dense or sparse matrix}
|
---|---|
Returns : | 2-norm condition number through use of the SVD : Use for small to moderate sized dense matrices. : For large sparse matrices, use condest. : |
Notes
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.
Examples
>>> import numpy
>>> from pyamg.util.linalg import cond
>>> c = condest(numpy.array([[1.0,0.],[0.,2.0]]))
>>> print c
2.0
Returns True if A is Hermitian to within tol
Parameters : | A : {dense or sparse matrix}
fast_check : {bool}
tol : {float}
verbose: {bool} :
|
---|---|
Returns : | True if hermitian : False if nonhermitian : |
Notes
This function applies a simple test of conjugate symmetry
Examples
>>> 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
Calculate the Moore-Penrose pseudo inverse of each block of the three dimensional array a.
Parameters : | a : {dense array}
cond : {float}
|
---|---|
Returns : | Nothing, a is modified in place so that a[k] holds the pseudoinverse : of that block. : |
Notes
By using lapack wrappers, this can be much faster for large n, than directly calling pinv2
Examples
>>> 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)
General utility functions for pyamg
Parameters : | A : sparse matrix or rank 1 array
|
---|---|
Returns : | B : array or sparse matrix
|
Examples
>>> 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.]]
A quick solver to profile a particular multilevel object
Parameters : | ml : multilevel
accel : function pointer
|
---|---|
Returns : | residuals : array
|
See also
multilevel.psolve, multilevel.solve
Examples
>>> 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)
Loop over all elements of varlist and convert them to upcasttype
Parameters : | upcast_type : data type
varlist : list
|
---|---|
Returns : | Returns upcast-ed varlist to upcast_type : |
Notes
Useful when harmonizing the types of variables, such as if A and b are complex, but x,y and z are not.
Examples
>>> 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])
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
Parameters : | upcast_type : data type
varlist : list
|
---|---|
Returns : | Returns upcast-ed varlist to upcast_type : |
Notes
Useful when harmonizing the types of variables, such as if A and b are complex, but x,y and z are not.
Examples
>>> 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])
Return the diagonal or inverse of diagonal for A, (A.H A) or (A A.H)
Parameters : | A : {dense or sparse matrix}
norm_eq : {0, 1, 2}
inv : {True, False}
|
---|---|
Returns : | diagonal, D, of appropriate system : |
Notes
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
Examples
>>> 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 ]
Unamalgamate a CSR A with blocks of 1’s.
Equivalent to Kronecker_Product(A, ones(RowsPerBlock, ColsPerBlock))
Parameters : | A : csr_matrix
RowsPerBlock : int
ColsPerBlock : int
|
---|---|
Returns : | A_UnAmal : bsr_matrix
|
Examples
>>> 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.]])
Convert 2D or 3D coordinates into Rigid body modes for use as near nullspace modes in elasticity AMG solvers
Parameters : | numNodes : int
numPDEs : :
x,y,z : array_like
|
---|---|
Returns : | rbm : matrix
|
Examples
>>> 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.]])
Examine a multilevel hierarchy’s spectrum
Parameters : | mg { pyamg multilevel hierarchy } :
|
---|---|
Returns : | (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 :
|
Notes
This can be useful for troubleshooting and when examining how your problem’s nature changes from level to level
Examples
>>> 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
Print a table from a list of lists representing the rows of a table
Parameters : | table : list
title : string
delim : string
centering : {‘left’, ‘right’, ‘center’}
col_padding : int
header : {True, False}
headerchar : {string}
|
---|---|
Returns : | string representing table that’s ready to be printed : |
Notes
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.
Examples
>>> 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')
Return the block diagonal of A, in array form
Parameters : | A : csr_matrix
blocksize : int
inv_flag : bool
|
---|---|
Returns : | block_diag : array
|
Examples
>>> 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)
Amalgamate matrix A
Parameters : | A : csr_matrix
blocksize : int
|
---|---|
Returns : | A_amal : csr_matrix
|
Notes
inverse operation of UnAmal for square matrices
Examples
>>> 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.]])
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.
Parameters : | A : sparse matrix
copy : {True,False}
|
---|---|
Returns : | D_sqrt : array
D_sqrt_inv : array
DAD : csr_matrix
|
Notes
Examples
>>> 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. ]]
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
Parameters : | A : {sparse matrix}
B : {array}
BH : {None, array}
|
---|---|
Returns : | 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 : |
Notes
Examples
>>> 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]]
Create a linear operator that applies a relaxation method for the given right-hand-side
Parameters : | methods : {tuple or string}
|
---|---|
Returns : | linear operator that applies the relaxation method to a vector for a : fixed right-hand-side, b. : |
Notes
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.
Examples
>>> 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
Filter the matrix A according to the matrix graph of C, while ensuring that the new, filtered A satisfies: A_new*B = Bf.
Returns : | A : sparse matrix updated such that sparsity structure of A now matches that of C, and that the relationship A*B = Bf holds. : |
---|
Notes
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
Examples
>>> 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. ]])
Helper function that scales T with a right multiplication by a block diagonal inverse, so that T is the identity at C-node rows.
Parameters : | T : {bsr_matrix}
P_I : {bsr_matrix}
I_F : {bsr_matrix}
|
---|---|
Returns : | T : {bsr_matrix}
|
Notes
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.
Examples
>>> 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. ]])
Helper function that returns a dictionary of sparse matrices and arrays which allow us to easily operate on Cpts and Fpts separately.
Parameters : | A : {csr_matrix, bsr_matrix}
Cnodes : {array}
AggOp : {csr_matrix}
T : {bsr_matrix}
|
---|---|
Returns : | Dictionary containing these parameters: : P_I : {bsr_matrix}
I_F : {bsr_matrix}
I_C : {bsr_matrix}
Cpts : {array}
Fpts : {array}
|
Notes
The principal calling routine is aggregation.smooth.energy_prolongation_smoother, which uses the Cpt_param dictionary for root-node style prolongation smoothing
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.
Parameters : | B : {array}
C : {csr_matrix, bsr_matrix}
|
---|---|
Returns : | BtBinv : {array}
|
Notes
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.