Util Documentation

This page contains the Util Package documentation.

The 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

Parameters :

A : bsr_matrix

Input matrix

i : int

Row number

Returns :

z : array

Actual nonzero values for row i colindx Array of column indices for the nonzeros of row i

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]
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

Parameters :

A : bsr_matrix

Input matrix

i : int

Row number

x : float

Scalar to overwrite nonzeros of row i in A

Returns :

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

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)
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.

Parameters :

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

Returns :

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

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]))

The 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)[source]

Approximate the spectral radius of a matrix

Parameters :

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.

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
pyamg.util.linalg.infinity_norm(A)[source]

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

Parameters :

A : csr_matrix, csc_matrix, sparse, or numpy matrix

Sparse or dense matrix

Returns :

n : float

Infinity norm of the matrix

See also

scipy.linalg.norm
dense matrix norms

Notes

  • This serves as an upper bound on spectral radius.
  • csr and csc avoid a deep copy
  • dense calls scipy.linalg.norm

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
pyamg.util.linalg.norm(x, pnorm='2')[source]

2-norm of a vector

Parameters :

x : array_like

Vector of complex or real values

pnorm : string

‘2’ calculates the 2-norm ‘inf’ calculates the infinity-norm

Returns :

n : float

2-norm of a vector

See also

scipy.linalg.norm
scipy general matrix or vector norm

Notes

  • 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
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

Parameters :

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

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
pyamg.util.linalg.cond(A)[source]

Returns condition number of A

Parameters :

A : {dense or sparse matrix}

e.g. array, matrix, csr_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
pyamg.util.linalg.ishermitian(A, fast_check=True, tol=1e-06, verbose=False)[source]

Returns True if A is Hermitian to within tol

Parameters :

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

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
pyamg.util.linalg.pinv_array(a, cond=None)[source]

Calculate the Moore-Penrose pseudo inverse of each block of the three dimensional array a.

Parameters :

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.

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)

The setup Module

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

The 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
Parameters :

A : sparse matrix or rank 1 array

General sparse matrix or array of diagonal entries

Returns :

B : array or sparse matrix

Diagonal sparse is returned as csr if A is dense otherwise return an array of the diagonal

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.]]
pyamg.util.utils.profile_solver(ml, accel=None, **kwargs)[source]

A quick solver to profile a particular multilevel object

Parameters :

ml : multilevel

Fully constructed multilevel object

accel : function pointer

Pointer to a valid Krylov solver (e.g. gmres, cg)

Returns :

residuals : array

Array of residuals for each iteration

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)
pyamg.util.utils.to_type(upcast_type, varlist)[source]

Loop over all elements of varlist and convert them to upcasttype

Parameters :

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 :

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])
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

Parameters :

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 :

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])
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)

Parameters :

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

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       ]
pyamg.util.utils.UnAmal(A, RowsPerBlock, ColsPerBlock)[source]

Unamalgamate a CSR A with blocks of 1’s.

Equivalent to Kronecker_Product(A, ones(RowsPerBlock, ColsPerBlock))

Parameters :

A : csr_matrix

Amalgamted matrix

RowsPerBlock : int

Give A blocks of size (RowsPerBlock, ColsPerBlock)

ColsPerBlock : int

Give A blocks of size (RowsPerBlock, ColsPerBlock)

Returns :

A_UnAmal : bsr_matrix

Similar to a Kronecker product of A and ones(RowsPerBlock, ColsPerBlock)

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.]])
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

Parameters :

numNodes : int

Number of nodes

numPDEs : :

Number of dofs per node

x,y,z : array_like

Coordinate vectors

Returns :

rbm : matrix

A matrix of size (numNodes*numPDEs) x (1 | 6) containing the 6 rigid body modes

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.]])
pyamg.util.utils.hierarchy_spectrum(mg, filter=True, plot=False)[source]

Examine a multilevel hierarchy’s spectrum

Parameters :

mg { pyamg multilevel hierarchy } :

e.g. generated with smoothed_aggregation_solver(...) or ruge_stuben_solver(...)

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 :

spectrum at each level

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 
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

Parameters :

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

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')
pyamg.util.utils.get_block_diag(A, blocksize, inv_flag=True)[source]

Return the block diagonal of A, in array form

Parameters :

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

Returns :

block_diag : array

block diagonal of A in array form, array size is (A.shape[0]/blocksize, blocksize, blocksize)

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)
pyamg.util.utils.amalgamate(A, blocksize)[source]

Amalgamate matrix A

Parameters :

A : csr_matrix

Matrix to amalgamate

blocksize : int

blocksize to use while amalgamating

Returns :

A_amal : csr_matrix

Amalgamated matrix A, first, convert A to BSR with square blocksize and then return CSR matrix using the resulting BSR indptr and indices

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.]])
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.

Parameters :

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)
Returns :

D_sqrt : array

Array of sqrt(diag(A))

D_sqrt_inv : array

Array of 1/sqrt(diag(A))

DAD : csr_matrix

Symmetrically scaled A

Notes

  • if A is not csr, it is converted to csr and sent to scale_rows

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. ]]
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
Parameters :

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.

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

  • if A is not csr, it is converted to csr and sent to scale_rows

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]]
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

Parameters :

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.

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
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.
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. ]])
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.

Parameters :

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

Returns :

T : {bsr_matrix}

Tentative prolongator scaled to be identity at C-pt nodes

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. ]])
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.

Parameters :

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

Returns :

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

Notes

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.

Parameters :

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.

Returns :

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.

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.

Table Of Contents

Previous topic

Testing Documentation

Next topic

Vis Documentation

This Page