pyamg.amg_core#

amg_core - a C++ implementation of AMG-related routines.

pyamg.amg_core.apply_absolute_distance_filter()#

Return a filtered strength-of-connection matrix by applying a drop tolerance.

Strength values are assumed to be “distance”-like, i.e. the smaller the value the stronger the connection. Strength values are _Not_ evaluated relatively, i.e. an off-diagonal entry A[i,j] is a strong connection iff:

S[i,j] <= epsilon,   where k != i

Also, set the diagonal to 1.0, as each node is perfectly close to itself.

Parameters:

n_row{int}: Dimension of matrix, S
epsilon{float}: Drop tolerance
Sp{int array}: Row pointer array for CSR matrix S
Sj{int array}: Col index array for CSR matrix S
Sx{float|complex array}: Value array for CSR matrix S

Returns:

Sx{float|complex array}: Modified in place such that the above dropping strategy has been applied There will be explicit zero entries for each weak connection

Notes

Principle calling routines are strength of connection routines, e.g., distance_strength_of_connection

Examples

>>> from scipy.sparse import csr_matrix
>>> from pyamg.amg_core import apply_absolute_distance_filter
>>> from scipy import array
>>> # Graph in CSR where entries in row i represent distances from dof i
>>> indptr = array([0,3,6,9])
>>> indices = array([0,1,2,0,1,2,0,1,2])
>>> data = array([1.,2.,3.,4.,1.,2.,3.,9.,1.])
>>> S = csr_matrix( (data,indices,indptr), shape=(3,3) )
>>> print "Matrix Before Applying Filter\n" + str(S.todense())
>>> apply_absolute_distance_filter(3, 1.9, S.indptr, S.indices, S.data)
>>> print "Matrix After Applying Filter\n" + str(S.todense())

pyamg.amg_core.apply_distance_filter()#

Return a filtered strength-of-connection matrix by applying a drop tolerance: Strength values are assumed to be “distance”-like, i.e. the smaller the value the stronger the connection

An off-diagonal entry A[i,j] is a strong connection iff

S[i,j] <= epsilon * min( S[i,k] ) where k != i

Also, set the diagonal to 1.0, as each node is perfectly close to itself

Parameters:

n_row{int}: Dimension of matrix, S
epsilon{float}: Drop tolerance
Sp{int array}: Row pointer array for CSR matrix S
Sj{int array}: Col index array for CSR matrix S
Sx{float|complex array}: Value array for CSR matrix S

Returns:

Sx{float|complex array}: Modified in place such that the above dropping strategy has been applied There will be explicit zero entries for each weak connection

Notes

Principle calling routines are strength of connection routines, e.g. evolution_strength_of_connection(…) in strength.py. It is used to apply a drop tolerance.

Examples

>>> from scipy.sparse import csr_matrix
>>> from pyamg.amg_core import apply_distance_filter
>>> from scipy import array
>>> # Graph in CSR where entries in row i represent distances from dof i
>>> indptr = array([0,3,6,9])
>>> indices = array([0,1,2,0,1,2,0,1,2])
>>> data = array([1.,2.,3.,4.,1.,2.,3.,9.,1.])
>>> S = csr_matrix( (data,indices,indptr), shape=(3,3) )
>>> print "Matrix before Applying Filter\n" + str(S.todense())
>>> apply_distance_filter(3, 1.9, S.indptr, S.indices, S.data)
>>> print "Matrix after Applying Filter\n" + str(S.todense())

pyamg.amg_core.apply_givens()#

Apply the first nrot Givens rotations in B to x

Parameters:

xarray: n-vector to be operated on
Barray: Each 4 entries represent a Givens rotation length nrot*4
nint: dimensionality of x
nrotint: number of rotations in B

Returns:

x is modified in place to reflect the application of the nrot
rotations in B. It is assumed that the first rotation operates on
degrees of freedom 0 and 1. The second rotation operates on dof’s 1 and 2,
and so on

Notes

Principle calling routine is gmres(…) and fgmres(…) in krylov.py

pyamg.amg_core.apply_householders()#

Apply Householder reflectors in B to z

Implements the below python

for j in range(start,stop,step):
  z = z - 2.0*dot(conjugate(B[j,:]), v)*B[j,:]

Parameters:

zarray: length n vector to be operated on
Barray: n x m matrix of householder reflectors must be in row major form
nint: dimensionality of z
start, stop, stepint: control the choice of vectors in B to use

Returns:

z is modified in place to reflect the application of
the Householder reflectors, B[:,range(start,stop,step)]

Notes

Principle calling routine is gmres(…) and fgmres(…) in krylov.py

pyamg.amg_core.approx_ideal_restriction_pass1()#

Build row_pointer for approximate ideal restriction in CSR or BSR form.

Parameters:

Rparray<int>: Empty row-pointer for R
Cpconst array<int>: Row pointer for SOC matrix, C
Cjconst array<int>: Column indices for SOC matrix, C
Cptsarray<int>: List of global C-point indices
splittingconst array<int>: Boolean array with 1 denoting C-points and 0 F-points
distanceint, default 2: Distance of F-point neighborhood to consider, options are 1 and 2.

Returns:

Nothing, Rp[] modified in place.

pyamg.amg_core.approx_ideal_restriction_pass2()#

Build column indices and data array for approximate ideal restriction in CSR format.

Parameters:

Rpconst array<int>: Pre-determined row-pointer for R in CSR format
Rjarray<int>: Empty array for column indices for R in CSR format
Rxarray<float>: Empty array for data for R in CSR format
Apconst array<int>: Row pointer for matrix A
Ajconst array<int>: Column indices for matrix A
Axconst array<float>: Data array for matrix A
Cpconst array<int>: Row pointer for SOC matrix, C
Cjconst array<int>: Column indices for SOC matrix, C
Cxconst array<float>: Data array for SOC matrix, C
Cptsarray<int>: List of global C-point indices
splittingconst array<int>: Boolean array with 1 denoting C-points and 0 F-points
distanceint, default 2: Distance of F-point neighborhood to consider, options are 1 and 2.
use_gmresbool, default 0: Use GMRES for local dense solve
maxiterint, default 10: Maximum GMRES iterations
preconditionbool, default True: Diagonally precondition GMRES

Returns:

Nothing, Rj[] and Rx[] modified in place.

Notes

Rx[] must be passed in initialized to zero.

pyamg.amg_core.bellman_ford()#

Apply one iteration of Bellman-Ford iteration on a distance graph stored in CSR format.

Parameters:

num_nodesint: number of nodes (number of rows in A)
Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array (edge lengths)
darray, inplace: distance to nearest center
cmarray, inplace: cluster index for each node

References

[1]

Bellman-Ford Wikipedia: http://en.wikipedia.org/wiki/Bellman-Ford_algorithm

pyamg.amg_core.block_approx_ideal_restriction_pass2()#

Build column indices and data array for approximate ideal restriction in BSR format.

Parameters:

Rpconst array<int>: Pre-determined row-pointer for R in CSR format
Rjarray<int>: Empty array for column indices for R in CSR format
Rxarray<float>: Empty array for data for R in CSR format
Apconst array<int>: Row pointer for matrix A
Ajconst array<int>: Column indices for matrix A
Axconst array<float>: Data array for matrix A
Cpconst array<int>: Row pointer for SOC matrix, C
Cjconst array<int>: Column indices for SOC matrix, C
Cxconst array<float>: Data array for SOC matrix, C
Cptsarray<int>: List of global C-point indices
splittingconst array<int>: Boolean array with 1 denoting C-points and 0 F-points
blocksizeint: Blocksize of matrix (assume square blocks)
distanceint, default 2: Distance of F-point neighborhood to consider, options are 1 and 2.
use_gmresbool, default 0: Use GMRES for local dense solve
maxiterint, default 10: Maximum GMRES iterations
preconditionbool, default True: Diagonally precondition GMRES

Returns:

Nothing, Rj[] and Rx[] modified in place.

Notes

Rx[] must be passed in initialized to zero.

pyamg.amg_core.block_gauss_seidel()#

Perform one iteration of block Gauss-Seidel relaxation on the linear system Ax = b, where A is stored in BSR format and x and b are column vectors.

Refer to gauss_seidel for additional information regarding row_start, row_stop, and row_step.

Parameters:

Aparray: BSR row pointer
Ajarray: BSR index array
Axarray: BSR data array, blocks assumed square
xarray, inplace: approximate solution
barray: right hand side
Txarray: Inverse of each diagonal block of A stored as a (n/blocksize, blocksize, blocksize) array
row_startint: beginning of the sweep
row_stopint: end of the sweep (i.e. one past the last unknown)
row_stepint: stride used during the sweep (may be negative)
blocksizeint: dimension of square blocks in BSR matrix A

pyamg.amg_core.block_jacobi()#

Perform one iteration of block Jacobi relaxation on the linear system Ax = b, where A is stored in BSR format and x and b are column vectors. Damping is controlled by the omega parameter.

Refer to gauss_seidel for additional information regarding row_start, row_stop, and row_step.

Parameters:

Aparray: BSR row pointer
Ajarray: BSR index array
Axarray: BSR data array, blocks assumed square
xarray, inplace: approximate solution
barray: right hand side
Txarray: Inverse of each diagonal block of A stored as a (n/blocksize, blocksize, blocksize) array
temparray: temporary vector the same size as x
row_startint: beginning of the sweep
row_stopint: end of the sweep (i.e. one past the last unknown)
row_stepint: stride used during the sweep (may be negative)
omegafloat: damping parameter
blocksize int: dimension of square blocks in BSR matrix A

pyamg.amg_core.block_jacobi_indexed()#

Perform one iteration of block Jacobi relaxation on the linear: system Ax = b for a given set of (block) row indices. A is stored in BSR format and x and b are column vectors. Damping is controlled by the parameter omega.

Parameters:

Aparray: BSR row pointer
Ajarray: BSR index array
Axarray: BSR data array, blocks assumed square
xarray: approximate solution
barray: right hand side
Txarray: Inverse of each diagonal block of A stored as a (n/blocksize, blocksize, blocksize) array
indicesarray: Indices
omegafloat: damping parameter
blocksizeint: dimension of square blocks in BSR matrix A

Returns:

Nothing, x will be modified in place

pyamg.amg_core.breadth_first_search()#

Compute a breadth first search of a graph in CSR format beginning at a given seed vertex.

Parameters:

num_rowsint: number of rows in A (number of vertices)
Aparray: CSR row pointer
Ajarray: CSR index array
orderarray, num_rows, inplace: records the order in which vertices were searched
levelarray, num_rows, inplace: records the level set of the searched vertices (i.e. the minimum distance to the seed)

Notes

The values of the level must be initialized to -1

pyamg.amg_core.bsr_gauss_seidel()#

Perform one iteration of Gauss-Seidel relaxation on the linear system Ax = b, where A is stored in Block CSR format and x and b are column vectors. This method applies point-wise relaxation to the BSR as opposed to "block relaxation".

Refer to gauss_seidel for additional information regarding row_start, row_stop, and row_step.

Parameters:

Aparray: BSR row pointer
Ajarray: BSR index array
Axarray: BSR data array
xarray, inplace: approximate solution
barray: right hand side
row_startint: beginning of the sweep (block row index)
row_stopint: end of the sweep (i.e. one past the last unknown)
row_stepint: stride used during the sweep (may be negative)
blocksizeint: BSR blocksize (blocks must be square)

Returns:

Nothing, x will be modified inplace

pyamg.amg_core.bsr_jacobi()#

Perform one iteration of Jacobi relaxation on the linear system Ax = b, where A is stored in Block CSR format and x and b are column vectors. This method applies point-wise relaxation to the BSR as opposed to "block relaxation".

Refer to jacobi for additional information regarding row_start, row_stop, and row_step.

Parameters:

Aparray: BSR row pointer
Ajarray: BSR index array
Axarray: BSR data array
xarray, inplace: approximate solution
barray: right hand side
temparray, inplace: temporary vector the same size as x
row_startint: beginning of the sweep (block row index)
row_stopint: end of the sweep (i.e. one past the last unknown)
row_stepint: stride used during the sweep (may be negative)
blocksizeint: BSR blocksize (blocks must be square)
omegafloat: damping parameter

Returns:

Nothing, x will be modified inplace

pyamg.amg_core.bsr_jacobi_indexed()#

Perform one iteration of Jacobi relaxation on the linear: system Ax = b for a given set of row indices, where A is stored in Block CSR format and x and b are column vectors. This method applies point-wise relaxation to the BSR matrix for a given set of row block indices, as opposed to “block relaxation”.

Parameters:

Aparray: BSR row pointer
Ajarray: BSR index array
Axarray: BSR data array
xarray: approximate solution
barray: right hand side
indicesarray: list of row indices to perform Jacobi on, e.g., F-points. Note, it is assumed that indices correspond to blocks in A.
blocksizeint: BSR blocksize (blocks must be square)
omegafloat: damping parameter

Returns:

Nothing, x will be modified in place

pyamg.amg_core.calc_BtB()#

Helper routine for energy_prolongation_smoother

Parameters:

NullDimint

Number of near nullspace vectors

Nnodesint

Number of nodes, i.e. number of block rows in BSR matrix, S

cols_per_blockint

Columns per block in S

barray

Nnodes x BsqCols array, in row-major form. This is B-squared, i.e. it is each column of B multiplied against each other column of B. For a Nx3 B,

b[:,0] = conjugate(B[:,0])*B[:,0]
b[:,1] = conjugate(B[:,0])*B[:,1]
b[:,2] = conjugate(B[:,0])*B[:,2]
b[:,3] = conjugate(B[:,1])*B[:,1]
b[:,4] = conjugate(B[:,1])*B[:,2]
b[:,5] = conjugate(B[:,2])*B[:,2]

BsqColsint

sum(range(NullDim+1)), i.e. number of columns in b

x{float|complex array}

Modified inplace for output. Should be zeros upon entry

Sp,Sjint array

BSR indptr and indices members for matrix, S

Returns:

BtB[i] = B_i.H*B_i in column major format
where B_i is B[colindices,:], colindices = all the nonzero
column indices for block row i in S

Notes

Principle calling routine is energy_prolongation_smoother(…) in smooth.py.

Calculates the following python code:

rows_per_block = Sparsity_Pattern.blocksize[0]
BtB = zeros((Nnodes,NullDim,NullDim), dtype=B.dtype)
S2 = Sparsity_Pattern.tocsr()
for i in range(Nnodes):
    Bi = mat( B[S2.indices[S2.indptr[i*rows_per_block]:S2.indptr[i*rows_per_block + 1]],:] )
    BtB[i,:,:] = Bi.H*Bi

pyamg.amg_core.classical_strength_of_connection_abs()#

Compute a strength of connection matrix using the classical strength: of connection measure by Ruge and Stuben. Both the input and output matrices are stored in CSR format. An off-diagonal nonzero entry A[i,j] is considered strong if:

Otherwise, the connection is weak.

Parameters:

num_rowsint: number of rows in A
thetafloat: strength of connection tolerance
Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array
Sparray: CSR row pointer
Sjarray: CSR index array
Sxarray: CSR data array

Returns:

Nothing, S will be stored in Sp, Sj, Sx

Notes

Storage for S must be preallocated. Since S will consist of a subset of A’s nonzero values, a conservative bound is to allocate the same storage for S as is used by A.

pyamg.amg_core.classical_strength_of_connection_min()#

pyamg.amg_core.cljp_naive_splitting()#

pyamg.amg_core.cluster_center()#

Apply Floyd–Warshall to cluster “a” and use the result to find the cluster center

Parameters:

aint: cluster index to find the center of
num_nodesint: number of nodes
num_clustersint: number of clusters
Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array (edge lengths)
cmarray, num_nodes: cluster index for each node
ICparray, num_clusters+1: CSC column pointer array for I
ICiarray, num_nodes: CSC column indexes for I
Larray, num_nodes: Local index mapping

Returns:

iint: global node index of center of cluster a

References

[1]

Graph Center: https://en.wikipedia.org/wiki/Graph_center

[2]

Floyd-Warshall: https://en.wikipedia.org/wiki/Floyd–Warshall_algorithm

[3]

Graph Distance: https://en.wikipedia.org/wiki/Distance_(graph_theory)

pyamg.amg_core.cluster_node_incidence()#

Compute the incidence matrix for a clustering

Parameters:

num_nodesint: number of nodes
num_clustersint: number of clusters
cmarray, num_nodes: cluster index for each node
ICparrayt, num_clusters+1, inplace: CSC column pointer array for I
ICiarray, num_nodes, inplace: CSC column indexes for I
Larray, num_nodes, inplace: Local index mapping

Notes

I = Incidence matrix between nodes and clusters (num_nodes x num_clusters) I[i,a] = 1 if node i is in cluster a, otherwise 0

Cluster indexes: a,b,c in 1..num_clusters Global node indexes: i,j,k in 1..num_rows Local node indexes: pair (a,m) where a is cluster and m in 1..num_nodes_in_cluster

We store I in both CSC and CSR formats because we want to be able to map global <-> local node indexes. However, I in CSR format is simply the cm array, so we only need to compute CSC format.

IC = (ICp,ICi) = I in CSC format (don’t store ICx because it’s always 1).

IR = (IRa) = (cm) = I in CSR format (don’t store IRp because we have exactly one nonzero entry per row, and don’t store IRx because it’s always 1). This is just the cm array.

Converting local (a,m) -> global i: i = ICi[ICp[a] + m] Converting global i -> local (a,m): a = cm[i], m = L[i]

L is an additional vector (length num_rows) to store local indexes.

pyamg.amg_core.connected_components()#

Compute the connected components of a graph stored in CSR format.

Parameters:

num_rowsint: number of rows in A (number of vertices)
Aparray: CSR row pointer
Ajarray: CSR index array
componentsarray, num_rows: component labels

Notes

Vertices belonging to each component are marked with a unique integer in the range [0,K), where K is the number of components.

pyamg.amg_core.cr_helper()#

Helper function for compatible relaxation to perform steps 3.1d - 3.1f in Falgout / Brannick (2010).

Parameters:

Aparray: Row pointer for sparse matrix in CSR format.
Ajarray: Column indices for sparse matrix in CSR format.
Barray: Target near null space vector for computing candidate set measure.
earray, inplace: Relaxed vector for computing candidate set measure.
indicesarray, inplace: Array of indices, where indices[0] = the number of F indices, nf, followed by F indices in elements 1:nf, and C indices in (nf+1):n.
splittingarray, inplace: Integer array with current C/F splitting of nodes, 0 = C-point, 1 = F-point.
gammaarray, inplace: Preallocated vector to store candidate set measure.
thetacsfloat: Threshold for coarse grid candidates from set measure.

Returns:

Nothing, updated C/F-splitting and corresponding indices modified in place.

pyamg.amg_core.csc_scale_columns()#

Scale the columns of a CSC matrix in place

References

scipy/scipy

pyamg.amg_core.csc_scale_rows()#

Scale the rows of a CSC matrix in place

References

scipy/scipy

pyamg.amg_core.evolution_strength_helper()#

Create strength-of-connection matrix based on constrained min problem of

min( z - B*x ), such that

(B*x)|_i = z|_i, i.e. they are equal at point i: z = (I - (t/k) Dinv A)^k delta_i

Strength is defined as the relative point-wise approximation error between B*x and z. B is the near-nullspace candidates. The constrained min problem is also restricted to consider B*x and z only at the nonzeros of column i of A

Can use either the D_A inner product, or l2 inner-prod in the minimization problem. Using D_A gives scale invariance. This choice is reflected in whether the parameter DB = B or diag(A)*B

This is a quadratic minimization problem with a linear constraint, so we can build a linear system and solve it to find the critical point, i.e. minimum.

Parameters:

Sp{int array}: Row pointer array for CSR matrix S
Sj{int array}: Col index array for CSR matrix S
Sx{float|complex array}: Value array for CSR matrix S. Upon entry to the routine, S = (I - (t/k) Dinv A)^k
nrows{int}: Dimension of S
B{float|complex array}: nrows x NullDim array of near nullspace vectors in col major form, if calling from within Python, take a transpose.
DB{float|complex array}: nrows x NullDim array of possibly scaled near nullspace vectors in col major form. If calling from within Python, take a transpose. For a scale invariant measure, DB = diag(A)*conjugate(B)), corresponding to the D_A inner-product Otherwise, DB = conjugate(B), corresponding to the l2-inner-product
b{float|complex array}: nrows x BDBCols array in row-major form. This array is B-squared, i.e. it is each column of B multiplied against each other column of B. For a Nx3 B, b[:,0] = conjugate(B[:,0])*B[:,0] b[:,1] = conjugate(B[:,0])*B[:,1] b[:,2] = conjugate(B[:,0])*B[:,2] b[:,3] = conjugate(B[:,1])*B[:,1] b[:,4] = conjugate(B[:,1])*B[:,2] b[:,5] = conjugate(B[:,2])*B[:,2]
BDBCols{int}: sum(range(NullDim+1)), i.e. number of columns in b
NullDim{int}: Number of nullspace vectors
tol{float}: Used to determine when values are numerically zero

Returns:

Sx{float|complex array}: Modified inplace and holds strength values for the above minimization problem

Notes

Upon entry to the routine, S = (I - (t/k) Dinv A)^k. However, we only need the values of S at the sparsity pattern of A. Hence, there is no need to completely calculate all of S.

b is used to save on the computation of each local minimization problem

Principle calling routine is evolution_strength_of_connection(…) in strength.py. In that routine, it is used to calculate strength-of-connection for the case of multiple near-nullspace modes.

Examples

See evolution_strength_of_connection(…) in strength.py

pyamg.amg_core.extract_subblocks()#

Extract diagonal blocks from A and insert into a linear array. This is a helper function for overlapping_schwarz_csr.

Parameters:

Aparray: CSR row pointer
Ajarray: CSR index array must be sorted for each row
Axarray: CSR data array, blocks assumed square
Txarray, inplace: Inverse of each diagonal block of A, stored in row major
Tparray: Pointer array into Tx indicating where the diagonal blocks start and stop
Sjarray: Indices of each subdomain must be sorted over each subdomain
Sparray: Pointer array indicating where each subdomain starts and stops
nsdomainsint: Number of subdomains
nrowsint: Number of rows

Returns:

Nothing, Tx will be modified inplace

pyamg.amg_core.filter_matrix_rows()#

Filter matrix rows by diagonal entry, that is set A_ij = 0 if:

|A_ij| < theta * |A_ii|

Parameters:

num_rowsint: number of rows in A
thetafloat: strength of connection tolerance
Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array

Returns:

Nothing, Ax is modified in place

pyamg.amg_core.fit_candidates()#

pyamg.amg_core.gauss_seidel()#

Perform one iteration of Gauss-Seidel relaxation on the linear system Ax = b, where A is stored in CSR format and x and b are column vectors.

Parameters:

Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array
xarray, inplace: approximate solution
barray: right hand side
row_startint: beginning of the sweep
row_stopint: end of the sweep (i.e. one past the last unknown)
row_stepint: stride used during the sweep (may be negative)

Returns:

Nothing, x will be modified inplace

Notes

The unknowns are swept through according to the slice defined by row_start, row_end, and row_step. These options are used to implement standard forward and backward sweeps, or sweeping only a subset of the unknowns. A forward sweep is implemented with gauss_seidel(Ap, Aj, Ax, x, b, 0, N, 1) where N is the number of rows in matrix A. Similarly, a backward sweep is implemented with gauss_seidel(Ap, Aj, Ax, x, b, N, -1, -1).

pyamg.amg_core.gauss_seidel_indexed()#

Perform one iteration of Gauss-Seidel relaxation on the linear system Ax = b, where A is stored in CSR format and x and b are column vectors.

Parameters:

Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array
xarray, inplace: approximate solution
barray: right hand side
Idarray: index array representing the
row_startint: beginning of the sweep (in array Id)
row_stopint: end of the sweep (in array Id)
row_stepint: stride used during the sweep (may be negative)

Returns:

Nothing, x will be modified inplace

Notes

Unlike gauss_seidel, which is restricted to updating a slice of the unknowns (defined by row_start, row_start, and row_step), this method updates unknowns according to the rows listed in an index array. This allows and arbitrary set of the unknowns to be updated in an arbitrary order, as is necessary for the relaxation steps in the Compatible Relaxation method.

In this method the slice arguments are used to define the subset of the index array Id which is to be considered.

pyamg.amg_core.gauss_seidel_ne()#

Perform NE Gauss-Seidel on the linear system A x = b This effectively carries out Gauss-Seidel on A A.H y = b, where x = A.h y.

Parameters:

Aparray: index pointer for CSR matrix A
Ajarray: column indices for CSR matrix A
Axarray: value array for CSR matrix A
xarray: current guess to the linear system
barray: right hand side
Txarray: inverse(diag(A A.H))
omegafloat: relaxation parameter (if not 1.0, then algorithm becomes SOR)
row_start,stop,stepint: controls which rows to iterate over

Returns:

x is modified inplace in an additive, not overwriting fashion

Notes

Primary calling routine is gauss_seidel_ne in relaxation.py

pyamg.amg_core.gauss_seidel_nr()#

Perform NR Gauss-Seidel on the linear system A x = b This effectively carries out Gauss-Seidel on A.H A x = A.H b

Parameters:

Aparray: index pointer for CSC matrix A
Ajarray: row indices for CSC matrix A
Axarray: value array for CSC matrix A
xarray: current guess to the linear system
zarray: initial residual
Txarray: inverse(diag(A.H A))
omegafloat: relaxation parameter (if not 1.0, then algorithm becomes SOR)
col_start,stop,stepint: controls which rows to iterate over

Returns:

x is modified inplace in an additive, not overwriting fashion

Notes

Primary calling routine is gauss_seidel_nr in relaxation.py

pyamg.amg_core.householder_hornerscheme()#

For use after gmres is finished iterating and the least squares solution has been found. This routine maps the solution back to the original space via the Householder reflectors.

Apply Householder reflectors in B to z while also adding in the appropriate value from y, so that we follow the Horner-like scheme to map our least squares solution in y back to the original space

Implements the below python

for j in range(inner,-1,-1):
  z[j] += y[j]
  # Apply j-th reflector, (I - 2.0*w_j*w_j.T)*update
  z = z - 2.0*dot(conjugate(B[j,:]), update)*B[j,:]

Parameters:

zarray: length n vector to be operated on
Barray: n x m matrix of householder reflectors must be in row major form
yarray: solution to the reduced system at the end of GMRES
nint: dimensionality of z
start, stop, stepint: control the choice of vectors in B to use

Returns:

z is modified in place to reflect the application of
the Householder reflectors, B[:,range(start,stop,step)],
and the inclusion of values in y.

Notes

Principle calling routine is gmres(…) and fgmres(…) in krylov.py

References

See pages 164-167 in Saad, “Iterative Methods for Sparse Linear Systems”

pyamg.amg_core.incomplete_mat_mult_bsr()#

Calculate A*B = S, but only at the pre-existing sparsity pattern of S, i.e. do an exact, but incomplete mat-mat mult.

A, B and S must all be in BSR, may be rectangular, but the indices need not be sorted. Also, A.blocksize[0] must equal S.blocksize[0], A.blocksize[1] must equal B.blocksize[0], and B.blocksize[1] must equal S.blocksize[1]

Parameters:

Ap{int array}: BSR row pointer array
Aj{int array}: BSR col index array
Ax{float|complex array}: BSR value array
Bp{int array}: BSR row pointer array
Bj{int array}: BSR col index array
Bx{float|complex array}: BSR value array
Sp{int array}: BSR row pointer array
Sj{int array}: BSR col index array
Sx{float|complex array}: BSR value array
n_brow{int}: Number of block-rows in A
n_bcol{int}: Number of block-cols in S
brow_A{int}: row blocksize for A
bcol_A{int}: column blocksize for A
bcol_B{int}: column blocksize for B

Returns:

Sx is modified in-place to reflect S(i,j) = <A_{i,:}, B_{:,j}>
but only for those entries already present in the sparsity pattern
of S.

Notes

Algorithm is SMMP

Principle calling routine is energy_prolongation_smoother(…) in smooth.py. Here it is used to calculate the descent direction A*P_tent, but only within an accepted sparsity pattern.

Is generally faster than the commented out incomplete_BSRmatmat(…) routine below, except when S has far few nonzeros than A or B.

pyamg.amg_core.incomplete_mat_mult_csr()#

Calculate A*B = S, but only at the pre-existing sparsity pattern of S, i.e. do an exact, but incomplete mat-mat multiply.

A must be in CSR, B must be in CSC and S must be in CSR Indices for A, B and S must be sorted A, B, and S must be square

Parameters:

Ap{int array}: Row pointer array for CSR matrix A
Aj{int array}: Col index array for CSR matrix A
Ax{float|complex array}: Value array for CSR matrix A
Bp{int array}: Row pointer array for CSC matrix B
Bj{int array}: Col index array for CSC matrix B
Bx{float|complex array}: Value array for CSC matrix B
Sp{int array}: Row pointer array for CSR matrix S
Sj{int array}: Col index array for CSR matrix S
Sx{float|complex array}: Value array for CSR matrix S
dimen: {int}: dimensionality of A,B and S

Returns:

Sx{float|complex array}: Modified inplace to reflect S(i,j) = <A_{i,:}, B_{:,j}>

Notes

A must be in CSR, B must be in CSC and S must be in CSR. Indices for A, B and S must all be sorted. A, B and S must be square.

Algorithm is naive, S(i,j) = <A_{i,:}, B_{:,j}> But, the routine is written for the case when S’s sparsity pattern is a subset of A*B, so this algorithm should work well.

Principle calling routine is evolution_strength_of_connection in strength.py. Here it is used to calculate S*S only at the sparsity pattern of the original operator. This allows for BIG cost savings.

Examples

>>> from pyamg.amg_core import incomplete_mat_mult_csr
>>> from scipy import arange, eye, ones
>>> from scipy.sparse import csr_matrix, csc_matrix
>>>
>>> A = csr_matrix(arange(1,10,dtype=float).reshape(3,3))
>>> B = csc_matrix(ones((3,3),dtype=float))
>>> AB = csr_matrix(eye(3,3,dtype=float))
>>> A.sort_indices()
>>> B.sort_indices()
>>> AB.sort_indices()
>>> incomplete_mat_mult_csr(A.indptr, A.indices, A.data, B.indptr, B.indices,
                      B.data, AB.indptr, AB.indices, AB.data, 3)
>>> print "Incomplete Matrix-Matrix Multiplication\n" + str(AB.todense())
>>> print "Complete Matrix-Matrix Multiplication\n" + str((A*B).todense())

pyamg.amg_core.jacobi()#

Perform one iteration of Jacobi relaxation on the linear system Ax = b, where A is stored in CSR format and x and b are column vectors. Damping is controlled by the omega parameter.

Refer to gauss_seidel for additional information regarding row_start, row_stop, and row_step.

Parameters:

Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array
xarray, inplace: approximate solution
barray: right hand side
temp, array: temporary vector the same size as x
row_startint: beginning of the sweep
row_stopint: end of the sweep (i.e. one past the last unknown)
row_stepint: stride used during the sweep (may be negative)
omegafloat: damping parameter

Returns:

Nothing, x will be modified inplace

pyamg.amg_core.jacobi_indexed()#

Perform one iteration of Jacobi relaxation on the linear

system Ax = b for a given set of row indices, where A is stored in CSR format and x and b are column vectors. Damping is controlled by the omega parameter.

Parameters: Ap[] - CSR row pointer Aj[] - CSR index array Ax[] - CSR data array x[] - approximate solution b[] - right hand side temp[] - temporary vector the same size as x indices[] - list of row indices to perform Jacobi on, e.g. F-points omega - damping parameter
Returns:: Nothing, x will be modified in place

pyamg.amg_core.jacobi_ne()#

Perform NE Jacobi on the linear system A x = b This effectively carries out weighted-Jacobi on A^TA x = A^T b (also known as Cimmino’s relaxation)

Parameters:

Aparray: index pointer for CSR matrix A
Ajarray: column indices for CSR matrix A
Axarray: value array for CSR matrix A
xarray, inplace: current guess to the linear system
barray: right hand side
Txarray: scaled residual D_A^{-1} (b - Ax)
temparray: work space
row_startint: controls which rows to start on
row_stopint: controls which rows to stop on
row_stepint: controls which rows to iterate over
omegaarray: size one array that contains the weighted-jacobi parameter. An array must be used to pass in omega to account for the case where omega may be complex

Returns:

x is modified inplace in an additive, not overwriting fashion

Notes

Primary calling routine is jacobi_ne in relaxation.py

pyamg.amg_core.lloyd_cluster()#

Perform one iteration of Lloyd clustering on a distance graph

Parameters:

num_nodesint: number of nodes (number of rows in A)
Aparray: CSR row pointer for adjacency matrix A
Ajarray: CSR index array
Axarray: CSR data array (edge lengths)
num_clustersint: number of clusters (seeds)
darray, num_nodes: distance to nearest seed
cmarray, num_nodes: cluster index for each node
carray, num_clusters: cluster centers

References

[Bell2008]

Nathan Bell, Algebraic Multigrid for Discrete Differential Forms PhD thesis (UIUC), August 2008

pyamg.amg_core.lloyd_cluster_exact()#

Perform one iteration of Lloyd clustering on a distance graph using exact centers

Parameters:

num_nodesint: number of rows in A (number of vertices)
Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array (edge lengths)
num_clustersint: number of clusters = number of seeds
darray, num_nodes: distance to nearest seed
cmarray, num_nodes: cluster index for each node
carray, num_clusters: cluster centers

Notes

This version computes exact cluster centers with Floyd-Warshall and also uses a balanced version of Bellman-Ford to try and find nearly-equal-sized clusters.

pyamg.amg_core.maximal_independent_set_k_parallel()#

Compute MIS-k.

Parameters:

num_rowsint: number of rows in A (number of vertices)
Aparray: CSR row pointer
Ajarray: CSR index array
kint: minimum separation between MIS vertices
xarray, inplace: state of each vertex (1 if in the MIS, 0 otherwise)
yarray: random values used during parallel MIS algorithm
max_itersint: maximum number of iterations to use (default, no limit)

Notes

Compute a distance-k maximal independent set for a graph stored in CSR format using a parallel algorithm. An MIS-k is a set of vertices such that all vertices in the MIS-k are separated by a path of at least K+1 edges and no additional vertex can be added to the set without destroying this property. A standard MIS is therefore a MIS-1.

pyamg.amg_core.maximal_independent_set_parallel()#

Compute a maximal independent set for a graph stored in CSR format: using a variant of Luby’s parallel MIS algorithm

Parameters:

num_rowsint: number of rows in A (number of vertices)
Aparray: CSR row pointer
Ajarray: CSR index array
activefloat: value used for active vertices
Cfloat: value used to mark non-MIS vertices
Ffloat: value used to mark MIS vertices
xarray, output: state of each vertex
yarray: random values for each vertex
max_itersint: maximum number of iterations By default max_iters=-1 and no limit is imposed

Returns:

Nint: The number of nodes in the MIS.

Notes

Only the vertices with values with x[i] == active are considered when determining the MIS. Upon return, all active vertices will be assigned the value C or F depending on whether they are in the MIS or not.

pyamg.amg_core.maximal_independent_set_serial()#

Compute a maximal independent set for a graph stored in CSR format using a greedy serial algorithm

Parameters:

num_rowsint: Number of rows in A (number of vertices)
Aparray: CSR row pointer
Ajarray: CSR index array
activefloat-like: Value used for active vertices
Cfloat-like: Value used to mark non-MIS vertices
Ffloat-like: Value used to mark MIS vertices
xarray, inplace output: State of each vertex

Returns:

Nint: The number of nodes in the MIS.

Notes

pyamg.amg_core.maximum_row_value()#

Compute the maximum in magnitude row value for a CSR matrix

Parameters:

num_rowsint: number of rows in A
xarray, inplace: num_rows array
Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array

Returns:

Nothing, x[i] will hold row i’s maximum magnitude entry

pyamg.amg_core.min_blocks()#

Given a BSR with num_blocks stored, return a linear array of length: num_blocks, which holds each block’s smallest, nonzero, entry

Parameters:

n_blocks{int}: Number of blocks in matrix
blocksize{int}: Size of each block
Sx{float|complex array}: Block data structure of BSR matrix, S Sx is (n_blocks x blocksize) in length
Tx{float|complex array}: modified in place for output

Returns:

Tx{float|complex array}: Modified in place; Tx[i] holds the minimum nonzero value of block i of S

Notes

Principle calling routine is evolution_strength_of_connection(…) in strength.py. In that routine, it is used to assign a strength of connection value between supernodes by setting the strength value to be the minimum nonzero in a block.

Examples

>>> from scipy.sparse import bsr_matrix, csr_matrix
>>> from pyamg.amg_core import min_blocks
>>> from numpy import zeros, array, ravel, round
>>> from scipy import rand
>>> row  = array([0,2,4,6])
>>> col  = array([0,2,2,0,1,2])
>>> data = round(10*rand(6,2,2), decimals=1)
>>> S = bsr_matrix( (data,col,row), shape=(6,6) )
>>> T = zeros(data.shape[0])
>>> print "Matrix before\n" + str(S.todense())
>>> min_blocks(6, 4, ravel(S.data), T)
>>> S2 = csr_matrix((T, S.indices, S.indptr), shape=(3,3))
>>> print "Matrix after\n" + str(S2.todense())

pyamg.amg_core.naive_aggregation()#

Compute aggregates for a matrix A stored in CSR format

Parameters:

n_rowint: number of rows in A
Aparray, n_row + 1: CSR row pointer
Ajarray, nnz: CSR column indices
xarray, n_row, inplace: aggregate numbers for each node
yarray, n_row, inplace: will hold Cpts upon return

Returns:

int: The number of aggregates (== max(x[:]) + 1)

Notes

Differs from standard aggregation. Each dof is considered. If it has been aggregated, skip over. Otherwise, put dof and any unaggregated neighbors in an aggregate. Results in possibly much higher complexities.

pyamg.amg_core.one_point_interpolation()#

Interpolate C-points by value and each F-point by value from its strongest connected C-neighbor.

Parameters:

Rpconst array<int>: Pre-determined row-pointer for P in CSR format
Rjarray<int>: Empty array for column indices for P in CSR format
Cpconst array<int>: Row pointer for SOC matrix, C
Cjconst array<int>: Column indices for SOC matrix, C
Cxconst array<float>: Data array for SOC matrix, C
splittingconst array<int>: Boolean array with 1 denoting C-points and 0 F-points

Returns:

Nothing, Rj[] modified in place.

pyamg.amg_core.overlapping_schwarz_csr()#

Perform one iteration of an overlapping Schwarz relaxation on the linear system Ax = b, where A is stored in CSR format and x and b are column vectors.

Refer to gauss_seidel for additional information regarding row_start, row_stop, and row_step.

Parameters:

Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array, blocks assumed square
xarray, inplace: approximate solution
barray: right hand side
Txarray: Inverse of each diagonal block of A, stored in row major
Tparray: Pointer array into Tx indicating where the diagonal blocks start and stop
Sjarray: Indices of each subdomain must be sorted over each subdomain
Sparray: Pointer array indicating where each subdomain starts and stops
nsdomains: Number of subdomains
nrows: Number of rows
row_startint: The subdomains are processed in this order,
row_stopint: for(i = row_start, i != row_stop, i+=row_step)
row_stepint: {…computation…}

Returns:

Nothing, x will be modified inplace

pyamg.amg_core.pairwise_aggregation()#

Compute aggregates for a matrix S stored in CSR format

Parameters:

n_rowint: number of rows in S
Sparray, n_row+1: CSR row pointer
Sjarray, nnz: CSR column indices
Sxarray, nnz: CSR data array
xarray, n_row, inplace: aggregate numbers for each node
yarray, n_row, inplace: will hold Cpts upon return

Returns:

int: The number of aggregates (== max(x[:]) + 1 )

Notes

S is the strength matrix. Assumes that the strength matrix is for classic strength with min norm.

pyamg.amg_core.pinv_array()#

Replace each block of A with a Moore-Penrose pseudoinverse of that block. Routine is designed to invert many small matrices at once.

Parameters:

AAarray: (m, n, n) array, assumed to be “raveled” and in row major form
m,nint: dimensions of AA
TransAchar: ‘T’ or ‘F’. Decides whether to transpose each nxn block of A before inverting. If using Python array, should be ‘T’.

Returns:

AAarray: AA is modified in place with the pseduoinverse replacing each block of AA. AA is returned in row-major form for Python

Notes

This routine is designed to be called once for a large m. Calling this routine repeatably would not be efficient.

This function offers substantial speedup over native Python code for many small matrices, e.g. 5x5 and 10x10. Tests have indicated that matrices larger than 27x27 are faster if done in native Python.

Examples

>>> from pyamg.amg_core import pinv_array
>>> from scipy import arange, ones, array, dot
>>> A = array([arange(1,5, dtype=float).reshape(2,2), ones((2,2),dtype=float)])
>>> Ac = A.copy()
>>> pinv_array(A, 2, 2, 'T')
>>> print "Multiplication By Inverse\n" + str(dot(A[0], Ac[0]))
>>> print "Multiplication by PseudoInverse\n" + str(dot(Ac[1], dot(A[1], Ac[1])))
>>>
>>> A = Ac.copy()
>>> pinv_array(A,2,2,'F')
>>> print "Changing flag to \'F\' results in different Inverse\n" + str(dot(A[0], Ac[0]))
>>> print "A holds the inverse of the transpose\n" + str(dot(A[0], Ac[0].T))

pyamg.amg_core.remove_strong_FF_connections()#

Remove strong F-to-F connections that do NOT have a common C-point from the set of strong connections. Specifically, set the data value in CSR format to 0. Removing zero entries afterwards will adjust row pointer and column indices.

Parameters:

n_nodesint: Number of rows in A
Sparray: Row pointer for SOC matrix, C
Sjarray: Column indices for SOC matrix, C
Sxarray: Data array for SOC matrix, C
splittingarray: Boolean array with 1 denoting C-points and 0 F-points

Returns:

Nothing, Sx[] is set to zero to eliminate connections.

pyamg.amg_core.rs_cf_splitting()#

Compute a C/F (coarse-fine( splitting using the classical coarse grid selection method of Ruge and Stuben. The strength of connection matrix S, and its transpose T, are stored in CSR format. Upon return, the splitting array will consist of zeros and ones, where C-nodes (coarse nodes) are marked with the value 1 and F-nodes (fine nodes) with the value 0.

Parameters:

n_nodesint: number of rows in A
Sparray: CSR row pointer array for SOC matrix
Sjarray: CSR column index array for SOC matrix
Tparray: CSR row pointer array for transpose of SOC matrix
Tjarray: CSR column index array for transpose of SOC matrix
influencearray: array that influences splitting (values stored here are added to lambda for each point)
splittingarray, inplace: array to store the C/F splitting

Notes

The splitting array must be preallocated

pyamg.amg_core.rs_cf_splitting_pass2()#

pyamg.amg_core.rs_classical_interpolation_pass1()#

First pass of classical AMG interpolation to build row pointer for P based on SOC matrix and CF-splitting.

Parameters:

n_nodesint: Number of rows in A
Sparray: Row pointer for SOC matrix, C
Sjarray: Column indices for SOC matrix, C
splittingarray: Boolean array with 1 denoting C-points and 0 F-points
Pparray: empty array to store row pointer for matrix P

Returns:

Nothing, Pp is modified in place.

pyamg.amg_core.rs_classical_interpolation_pass2()#

Produce a classical AMG interpolation operator for the case in which two strongly connected F -points do NOT have a common C-neighbor. Formula can be found in Sec. 3 Eq. (8) of [1] for modified=False and Eq. (9) for modified=True.

Parameters:

Aparray: Row pointer for matrix A
Ajarray: Column indices for matrix A
Axarray: Data array for matrix A
Sparray: Row pointer for SOC matrix, C
Sjarray: Column indices for SOC matrix, C
Sxarray: Data array for SOC matrix, C – MUST HAVE VALUES OF A
splittingarray: Boolean array with 1 denoting C-points and 0 F-points
Pparray: Row pointer for matrix P
Pjarray: Column indices for matrix P
Pxarray: Data array for matrix P
modifiedbool: Use modified interpolation formula

Returns:

Nothing, Pj[] and Px[] modified in place.

Notes

For modified interpolation, it is assumed that SOC matrix C is passed in WITHOUT any F-to-F connections that do not share a common C-point neighbor. Any SOC matrix C can be set as such by calling remove_strong_FF_connections().

References

[0] V. E. Henson and U. M. Yang, BoomerAMG: a parallel algebraic multigrid

solver and preconditioner, Applied Numerical Mathematics 41 (2002).

[1] “Distance-Two Interpolation for Parallel Algebraic Multigrid,”

De Sterck, R. Falgout, J. Nolting, U. M. Yang, (2008).

pyamg.amg_core.rs_direct_interpolation_pass1()#

Produce the Ruge-Stuben prolongator using “Direct Interpolation”

The first pass uses the strength of connection matrix ‘S’ and C/F splitting to compute the row pointer for the prolongator.

The second pass fills in the nonzero entries of the prolongator

Parameters:

n_nodesint: Number of nodes
Sparray: Strength matrix row pointer array
Sjarray: Strength matrix column index array
splittingarray: C/F splitting
Pparray, inplace: Row pointer array

References

Page 479 of Multigrid

pyamg.amg_core.rs_direct_interpolation_pass2()#

pyamg.amg_core.satisfy_constraints_helper()#

Helper routine for satisfy_constraints routine called by energy_prolongation_smoother(…) in smooth.py

Parameters:

rows_per_blockint: rows per block in the BSR matrix, S
cols_per_blockint: cols per block in the BSR matrix, S
num_block_rowsint: Number of block rows, S.shape[0]/rows_per_block
NullDimint: Null-space dimension, i.e., the number of columns in B
xarray: Conjugate of near-nullspace vectors, B, in row major
yarray: S*B, in row major
zarray: BtBinv, in row major, i.e. z[i] = pinv(B_i.H Bi), where B_i is B restricted to the neighborhood of dof of i.
Sparray: Row pointer array for BSR matrix S
Sjarray: Col index array for BSR matrix S
Sxarray: Value array for BSR matrix S

Returns:

Sx is modified such that S*B = 0. S ends up being the
update to the prolongator in the energy_minimization algorithm.

Notes

Principle calling routine is energy_prolongation_smoother(…) in smooth.py.

This implements the python code:

# U is a BSR matrix, B is num_block_rows x cols_per_block x cols_per_block
# UB is num_block_rows x rows_per_block x cols_per_block,  BtBinv is
     num_block_rows x cols_per_block x cols_per_block
B  = asarray(B).reshape(-1,cols_per_block,B.shape[1])
UB = asarray(UB).reshape(-1,rows_per_block,UB.shape[1])
rows = csr_matrix((U.indices,U.indices,U.indptr), \
        shape=(U.shape[0]/rows_per_block,U.shape[1]/cols_per_block)).tocoo(copy=False).row
for n,j in enumerate(U.indices):
   i = rows[n]
   Bi  = mat(B[j])
   UBi = UB[i]
   U.data[n] -= dot(UBi,dot(BtBinv[i],Bi.H))

pyamg.amg_core.standard_aggregation()#

Compute aggregates for a matrix A stored in CSR format

Parameters:

n_rowint: number of rows in A
Aparray, n_row + 1: CSR row pointer
Ajarray, nnz: CSR column indices
xarray, n_row, inplace: aggregate numbers for each node
yarray, n_row, inplace: will hold Cpts upon return

Returns:

int: The number of aggregates (== max(x[:]) + 1)

Notes

It is assumed that A is symmetric.
A may contain diagonal entries (self loops)
Unaggregated nodes are marked with a -1

pyamg.amg_core.symmetric_strength_of_connection()#

Compute a strength of connection matrix using the standard symmetric Smoothed Aggregation heuristic. Both the input and output matrices are stored in CSR format. A nonzero connection A[i,j] is considered strong if:

The strength of connection matrix S is simply the set of nonzero entries of A that qualify as strong connections.

Parameters:

num_rowsint: number of rows in A
thetafloat: strength of connection tolerance
Aparray: CSR row pointer
Ajarray: CSR index array
Axarray: CSR data array
Sparray, inplace: CSR row pointer
Sjarray, inplace: CSR index array
Sxarray, inplace: CSR data array

Notes

Storage for S must be preallocated. Since S will consist of a subset of A’s nonzero values, a conservative bound is to allocate the same storage for S as is used by A.

pyamg.amg_core.truncate_rows_csr()#

Truncate the entries in A, such that only the largest (in magnitude)

k entries per row are left. Smaller entries are zeroed out.

Parameters: n_row - number of rows in A k - number of entries per row to keep Sp[] - CSR row pointer Sj[] - CSR index array Sx[] - CSR data array
Returns:: Nothing, A will be stored in Sp, Sj, Sx with some entries zeroed out

pyamg.amg_core.vertex_coloring_LDF()#

Compute a vertex coloring of a graph using the parallel Largest-Degree-First (LDF) algorithm

Parameters:

num_rowsint: number of rows in A (number of vertices)
Aparray: CSR row pointer
Ajarray: CSR index array
xarray: color of each vertex
yarray: initial random values for each vertex

References

[LDF]

J. R. Allwright and R. Bordawekar and P. D. Coddington and K. Dincer and C. L. Martin A Comparison of Parallel Graph Coloring Algorithms DRAFT SCCS-666 http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.4650

pyamg.amg_core.vertex_coloring_jones_plassmann()#

Compute a vertex coloring of a graph using the Jones-Plassmann algorithm

Parameters:

num_rowsint: number of rows in A (number of vertices)
Aparray: CSR row pointer
Ajarray: CSR index array
xarray, inplace: color of each vertex
yarray: initial random values for each vertex

Notes

Arrays x and y will be overwritten

References

[Jones92]

Mark T. Jones and Paul E. Plassmann A Parallel Graph Coloring Heuristic SIAM Journal on Scientific Computing 14:3 (1993) 654–669 http://citeseer.ist.psu.edu/jones92parallel.html

pyamg.amg_core.vertex_coloring_mis()#: Compute a vertex coloring for a graph stored in CSR format.

The coloring is computed by removing maximal independent sets of vertices from the graph. Specifically, at iteration i an independent set of the remaining subgraph is constructed and assigned color i.

Returns the K, the number of colors used in the coloring. On return x[i] in [0,1, …, K - 1] will contain the color of the i-th vertex.