# 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:
**x**arrayn-vector to be operated on

**B**arrayEach 4 entries represent a Givens rotation length nrot*4

**n**intdimensionality of x

**nrot**intnumber 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:
**z**arraylength n vector to be operated on

**B**arrayn x m matrix of householder reflectors must be in row major form

**n**intdimensionality of z

**start, stop, step**intcontrol 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:
**Rp**array<int>Empty row-pointer for R

**Cp**const array<int>Row pointer for SOC matrix, C

**Cj**const array<int>Column indices for SOC matrix, C

**Cpts**array<int>List of global C-point indices

**splitting**const array<int>Boolean array with 1 denoting C-points and 0 F-points

**distance**int, default 2Distance 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:
**Rp**const array<int>Pre-determined row-pointer for R in CSR format

**Rj**array<int>Empty array for column indices for R in CSR format

**Rx**array<float>Empty array for data for R in CSR format

**Ap**const array<int>Row pointer for matrix A

**Aj**const array<int>Column indices for matrix A

**Ax**const array<float>Data array for matrix A

**Cp**const array<int>Row pointer for SOC matrix, C

**Cj**const array<int>Column indices for SOC matrix, C

**Cx**const array<float>Data array for SOC matrix, C

**Cpts**array<int>List of global C-point indices

**splitting**const array<int>Boolean array with 1 denoting C-points and 0 F-points

**distance**int, default 2Distance of F-point neighborhood to consider, options are 1 and 2.

**use_gmres**bool, default 0Use GMRES for local dense solve

**maxiter**int, default 10Maximum GMRES iterations

**precondition**bool, default TrueDiagonally 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_nodes**intnumber of nodes (number of rows in A)

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR data array (edge lengths)

**d**array, inplacedistance to nearest center

**cm**array, inplacecluster 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 CSR format.

- Parameters:
**Rp**const array<int>Pre-determined row-pointer for R in CSR format

**Rj**array<int>Empty array for column indices for R in CSR format

**Rx**array<float>Empty array for data for R in CSR format

**Ap**const array<int>Row pointer for matrix A

**Aj**const array<int>Column indices for matrix A

**Ax**const array<float>Data array for matrix A

**Cp**const array<int>Row pointer for SOC matrix, C

**Cj**const array<int>Column indices for SOC matrix, C

**Cx**const array<float>Data array for SOC matrix, C

**Cpts**array<int>List of global C-point indices

**splitting**const array<int>Boolean array with 1 denoting C-points and 0 F-points

**blocksize**intBlocksize of matrix (assume square blocks)

**distance**int, default 2Distance of F-point neighborhood to consider, options are 1 and 2.

**use_gmres**bool, default 0Use GMRES for local dense solve

**maxiter**int, default 10Maximum GMRES iterations

**precondition**bool, default TrueDiagonally 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:
**Ap**arrayBSR row pointer

**Aj**arrayBSR index array

**Ax**arrayBSR data array, blocks assumed square

**x**array, inplaceapproximate solution

**b**arrayright hand side

**Tx**arrayInverse of each diagonal block of A stored as a (n/blocksize, blocksize, blocksize) array

**row_start**intbeginning of the sweep

**row_stop**intend of the sweep (i.e. one past the last unknown)

**row_step**intstride used during the sweep (may be negative)

**blocksize**intdimension 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:
**Ap**arrayBSR row pointer

**Aj**arrayBSR index array

**Ax**arrayBSR data array, blocks assumed square

**x**array, inplaceapproximate solution

**b**arrayright hand side

**Tx**arrayInverse of each diagonal block of A stored as a (n/blocksize, blocksize, blocksize) array

**temp**arraytemporary vector the same size as x

**row_start**intbeginning of the sweep

**row_stop**intend of the sweep (i.e. one past the last unknown)

**row_step**intstride used during the sweep (may be negative)

**omega**floatdamping parameter

**blocksize int**dimension of sqare 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:
**Ap**arrayBSR row pointer

**Aj**arrayBSR index array

**Ax**arrayBSR data array, blocks assumed square

**x**arrayapproximate solution

**b**arrayright hand side

**Tx**arrayInverse of each diagonal block of A stored as a (n/blocksize, blocksize, blocksize) array

**indices**arrayIndices

**omega**floatdamping parameter

**blocksize**intdimension 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_rows**intnumber of rows in A (number of vertices)

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**order**array, num_rows, inplacerecords the order in which vertices were searched

**level**array, num_rows, inplacerecords 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:
**Ap**arrayBSR row pointer

**Aj**arrayBSR index array

**Ax**arrayBSR data array

**x**array, inplaceapproximate solution

**b**arrayright hand side

**row_start**intbeginning of the sweep (block row index)

**row_stop**intend of the sweep (i.e. one past the last unknown)

**row_step**intstride used during the sweep (may be negative)

**blocksize**intBSR 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:
**Ap**arrayBSR row pointer

**Aj**arrayBSR index array

**Ax**arrayBSR data array

**x**array, inplaceapproximate solution

**b**arrayright hand side

**temp**array, inplacetemporary vector the same size as x

**row_start**intbeginning of the sweep (block row index)

**row_stop**intend of the sweep (i.e. one past the last unknown)

**row_step**intstride used during the sweep (may be negative)

**blocksize**intBSR blocksize (blocks must be square)

**omega**floatdamping 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:
**Ap**arrayBSR row pointer

**Aj**arrayBSR index array

**Ax**arrayBSR data array

**x**arrayapproximate solution

**b**arrayright hand side

**indices**arraylist of row indices to perform Jacobi on, e.g., F-points. Note, it is assumed that indices correspond to blocks in A.

**blocksize**intBSR blocksize (blocks must be square)

**omega**floatdamping parameter

- Returns:
- Nothing, x will be modified in place

- pyamg.amg_core.calc_BtB()#
Helper routine for energy_prolongation_smoother

- Parameters:
**NullDim**intNumber of near nullspace vectors

**Nnodes**intNumber of nodes, i.e. number of block rows in BSR matrix, S

**cols_per_block**intColumns per block in S

**b**arrayNnodes 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]

**BsqCols**intsum(range(NullDim+1)), i.e. number of columns in b

**x**{float|complex array}Modified inplace for output. Should be zeros upon entry

**Sp,Sj**int arrayBSR 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_rows**intnumber of rows in A

**theta**floatstength of connection tolerance

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR data array

**Sp**arrayCSR row pointer

**Sj**arrayCSR index array

**Sx**arrayCSR 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:
**a**intcluster index to find the center of

**num_nodes**intnumber of nodes

**num_clusters**intnumber of clusters

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR data array (edge lengths)

**cm**array, num_nodescluster index for each node

**ICp**array, num_clusters+1CSC column pointer array for I

**ICi**array, num_nodesCSC column indexes for I

**L**array, num_nodesLocal index mapping

- Returns:
**i**intglobal 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_nodes**intnumber of nodes

**num_clusters**intnumber of clusters

**cm**array, num_nodescluster index for each node

**ICp**arrayt, num_clusters+1, inplaceCSC column pointer array for I

**ICi**array, num_nodes, inplaceCSC column indexes for I

**L**array, num_nodes, inplaceLocal 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_rows**intnumber of rows in A (number of vertices)

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**components**array, num_rowscomponent 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:
**Ap**arrayRow pointer for sparse matrix in CSR format.

**Aj**arrayColumn indices for sparse matrix in CSR format.

**B**arrayTarget near null space vector for computing candidate set measure.

**e**array, inplaceRelaxed vector for computing candidate set measure.

**indices**array, inplaceArray 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.

**splitting**array, inplaceInteger array with current C/F splitting of nodes, 0 = C-point, 1 = F-point.

**gamma**array, inplacePreallocated vector to store candidate set measure.

**thetacs**floatThreshold 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

- pyamg.amg_core.csc_scale_rows()#
Scale the rows of a CSC matrix

*in place*References

- 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:
**Ap**arrayCSR row pointer

**Aj**arrayCSR index array must be sorted for each row

**Ax**arrayCSR data array, blocks assumed square

**Tx**array, inplaceInverse of each diagonal block of A, stored in row major

**Tp**arrayPointer array into Tx indicating where the diagonal blocks start and stop

**Sj**arrayIndices of each subdomain must be sorted over each subdomain

**Sp**arrayPointer array indicating where each subdomain starts and stops

**nsdomains**intNumber of subdomains

**nrows**intNumber 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_rows**intnumber of rows in A

**theta**floatstength of connection tolerance

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR 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:
**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR data array

**x**array, inplaceapproximate solution

**b**arrayright hand side

**row_start**intbeginning of the sweep

**row_stop**intend of the sweep (i.e. one past the last unknown)

**row_step**intstride 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:
**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR data array

**x**array, inplaceapproximate solution

**b**arrayright hand side

**Id**arrayindex array representing the

**row_start**intbeginning of the sweep (in array Id)

**row_stop**intend of the sweep (in array Id)

**row_step**intstride 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:
**Ap**arrayindex pointer for CSR matrix A

**Aj**arraycolumn indices for CSR matrix A

**Ax**arrayvalue array for CSR matrix A

**x**arraycurrent guess to the linear system

**b**arrayright hand side

**Tx**arrayinverse(diag(A A.H))

**omega**floatrelaxation parameter (if not 1.0, then algorithm becomes SOR)

**row_start,stop,step**intcontrols which rows to iterate over

- Returns:
- x is modified inplace in an additive, not overwriting fashion

Notes

Primary calling routine is gass_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:
**Ap**arrayindex pointer for CSC matrix A

**Aj**arrayrow indices for CSC matrix A

**Ax**arrayvalue array for CSC matrix A

**x**arraycurrent guess to the linear system

**z**arrayinitial residual

**Tx**arrayinverse(diag(A.H A))

**omega**floatrelaxation parameter (if not 1.0, then algorithm becomes SOR)

**col_start,stop,step**intcontrols 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:
**z**arraylength n vector to be operated on

**B**arrayn x m matrix of householder reflectors must be in row major form

**y**arraysolution to the reduced system at the end of GMRES

**n**intdimensionality of z

**start, stop, step**intcontrol 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:
**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR data array

**x**array, inplaceapproximate solution

**b**arrayright hand side

**temp, array**temporary vector the same size as x

**row_start**intbeginning of the sweep

**row_stop**intend of the sweep (i.e. one past the last unknown)

**row_step**intstride used during the sweep (may be negative)

**omega**floatdamping 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:
**Ap**arrayindex pointer for CSR matrix A

**Aj**arraycolumn indices for CSR matrix A

**Ax**arrayvalue array for CSR matrix A

**x**array, inplacecurrent guess to the linear system

**b**arrayright hand side

**Tx**arrayscaled residual D_A^{-1} (b - Ax)

**temp**arraywork space

**row_start**intcontrols which rows to start on

**row_stop**intcontrols which rows to stop on

**row_step**intcontrols which rows to iterate over

**omega**arraysize 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_nodes**intnumber of nodes (number of rows in A)

**Ap**arrayCSR row pointer for adjacency matrix A

**Aj**arrayCSR index array

**Ax**arrayCSR data array (edge lengths)

**num_clusters**intnumber of clusters (seeds)

**d**array, num_nodesdistance to nearest seed

**cm**array, num_nodescluster index for each node

**c**array, num_clusterscluster 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_nodes**intnumber of rows in A (number of vertices)

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR data array (edge lengths)

**num_clusters**intnumber of clusters = number of seeds

**d**array, num_nodesdistance to nearest seed

**cm**array, num_nodescluster index for each node

**c**array, num_clusterscluster 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_rows**intnumber of rows in A (number of vertices)

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**k**intminimum separation between MIS vertices

**x**array, inplacestate of each vertex (1 if in the MIS, 0 otherwise)

**y**arrayrandom values used during parallel MIS algorithm

**max_iters**intmaximum 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_rows**intnumber of rows in A (number of vertices)

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**active**floatvalue used for active vertices

**C**floatvalue used to mark non-MIS vertices

**F**floatvalue used to mark MIS vertices

**x**array, outputstate of each vertex

**y**arrayrandom values for each vertex

**max_iters**intmaximum number of iterations By default max_iters=-1 and no limit is imposed

- Returns:
**N**intThe 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_rows**intNumber of rows in A (number of vertices)

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**active**float-likeValue used for active vertices

**C**float-likeValue used to mark non-MIS vertices

**F**float-likeValue used to mark MIS vertices

**x**array, inplace outputState of each vertex

- Returns:
**N**intThe 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.maximum_row_value()#
Compute the maximum in magnitude row value for a CSR matrix

- Parameters:
**num_rows**intnumber of rows in A

**x**array, inplacenum_rows array

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR 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_row**intnumber of rows in A

**Ap**array, n_row + 1CSR row pointer

**Aj**array, nnzCSR column indices

**x**array, n_row, inplaceaggregate numbers for each node

**y**array, n_row, inplacewill 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:
**Rp**const array<int>Pre-determined row-pointer for P in CSR format

**Rj**array<int>Empty array for column indices for P in CSR format

**Cp**const array<int>Row pointer for SOC matrix, C

**Cj**const array<int>Column indices for SOC matrix, C

**Cx**const array<float>Data array for SOC matrix, C

**splitting**const 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:
**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR data array, blocks assumed square

**x**array, inplaceapproximate solution

**b**arrayright hand side

**Tx**arrayInverse of each diagonal block of A, stored in row major

**Tp**arrayPointer array into Tx indicating where the diagonal blocks start and stop

**Sj**arrayIndices of each subdomain must be sorted over each subdomain

**Sp**arrayPointer array indicating where each subdomain starts and stops

**nsdomains**Number of subdomains

**nrows**Number of rows

**row_start**intThe subdomains are processed in this order,

**row_stop**intfor(i = row_start, i != row_stop, i+=row_step)

**row_step**int{…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_row**intnumber of rows in S

**Sp**array, n_row+1CSR row pointer

**Sj**array, nnzCSR column indices

**Sx**array, nnzCSR data array

**x**array, n_row, inplaceaggregate numbers for each node

**y**array, n_row, inplacewill 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:
**AA**array(m, n, n) array, assumed to be “raveled” and in row major form

**m,n**intdimensions of AA

**TransA**char‘T’ or ‘F’. Decides whether to transpose each nxn block of A before inverting. If using Python array, should be ‘T’.

- Returns:
**AA**arrayAA 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_nodes**intNumber of rows in A

**Sp**arrayRow pointer for SOC matrix, C

**Sj**arrayColumn indices for SOC matrix, C

**Sx**arrayData array for SOC matrix, C

**splitting**arrayBoolean 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_nodes**intnumber of rows in A

**Sp**arrayCSR row pointer array for SOC matrix

**Sj**arrayCSR column index array for SOC matrix

**Tp**arrayCSR row pointer array for transpose of SOC matrix

**Tj**arrayCSR column index array for transpose of SOC matrix

**influence**arrayarray that influences splitting (values stored here are added to lambda for each point)

**splitting**array, inplacearray 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_nodes**intNumber of rows in A

**Sp**arrayRow pointer for SOC matrix, C

**Sj**arrayColumn indices for SOC matrix, C

**splitting**arrayBoolean array with 1 denoting C-points and 0 F-points

**Pp**arrayempty 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:
**Ap**arrayRow pointer for matrix A

**Aj**arrayColumn indices for matrix A

**Ax**arrayData array for matrix A

**Sp**arrayRow pointer for SOC matrix, C

**Sj**arrayColumn indices for SOC matrix, C

**Sx**arrayData array for SOC matrix, C – MUST HAVE VALUES OF A

**splitting**arrayBoolean array with 1 denoting C-points and 0 F-points

**Pp**arrayRow pointer for matrix P

**Pj**arrayColumn indices for matrix P

**Px**arrayData array for matrix P

**modified**boolUse 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_nodes**intNumber of nodes

**Sp**arrayStrength matrix row pointer array

**Sj**arrayStrength matrix column index array

**splitting**arrayC/F splitting

**Pp**array, inplaceRow 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_block**introws per block in the BSR matrix, S

**cols_per_block**intcols per block in the BSR matrix, S

**num_block_rows**intNumber of block rows, S.shape[0]/rows_per_block

**NullDim**intNull-space dimension, i.e., the number of columns in B

**x**arrayConjugate of near-nullspace vectors, B, in row major

**y**arrayS*B, in row major

**z**arrayBtBinv, 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.

**Sp**arrayRow pointer array for BSR matrix S

**Sj**arrayCol index array for BSR matrix S

**Sx**arrayValue 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_row**intnumber of rows in A

**Ap**array, n_row + 1CSR row pointer

**Aj**array, nnzCSR column indices

**x**array, n_row, inplaceaggregate numbers for each node

**y**array, n_row, inplacewill 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_rows**intnumber of rows in A

**theta**floatstength of connection tolerance

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**Ax**arrayCSR data array

**Sp**array, inplaceCSR row pointer

**Sj**array, inplaceCSR index array

**Sx**array, inplaceCSR 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_rows**intnumber of rows in A (number of vertices)

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**x**arraycolor of each vertex

**y**arrayinitial 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_rows**intnumber of rows in A (number of vertices)

**Ap**arrayCSR row pointer

**Aj**arrayCSR index array

**x**array, inplacecolor of each vertex

**y**arrayinitial 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.