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_rowint: Dimension of matrix, S.
epsilonfloat: Drop tolerance.
Sparray: Row pointer array for CSR matrix S.
Sjint array: Col index array for CSR matrix S.
Sxarray: Value array for CSR matrix S.

Returns:

array: Modified in place such that the above dropping strategy has been applied. There will be explicit zero entries for each weak connection.

See also

distance_strength_of_connection: Principle calling routines are strength of connection routines.

Examples

>>> from scipy.sparse import csr_array
>>> 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_array( (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_rowint: Dimension of matrix, S.
epsilonfloat: Drop tolerance.
Sparray: Row pointer array for CSR matrix S.
Sjarray: Col index array for CSR matrix S.
Sxarray: Value array for CSR matrix S.

Returns:

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_array
>>> 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_array( (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: Vector of length n 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:

None: 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 routines are gmres and fgmres.

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: Vector of length n to be operated on.
Barray: Matrix of size (n,n) 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:

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

Notes

Principle calling routines are gmres and fgmres.

pyamg.amg_core.approx_ideal_restriction_pass1()#

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

Parameters:

Rparray: Empty row-pointer for R.
Cparray: Row pointer for SOC matrix, C.
Cjarray: Column indices for SOC matrix, C.
Cptsarray: List of global C-point indices.
splittingarray: 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:

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

Rparray: Pre-determined row-pointer for R in CSR format.
Rjarray: Empty array for column indices for R in CSR format.
Rxarray: Empty array for data for R in CSR format.
Aparray: Row pointer for matrix A.
Ajarray: Column indices for matrix A.
Axarray: Data array for matrix A.
Cparray: Row pointer for SOC matrix, C.
Cjarray: Column indices for SOC matrix, C.
Cxarray: Data array for SOC matrix, C.
Cptsarray: List of global C-point indices.
splittingarray: 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:

None: 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).
carray: Cluster center.
darray, inplace: Distance to nearest center.
marray, inplace: Cluster index for each node.
parray, inplace: Predecssor on the shortest path to center.

See also

pyamg.graph.bellman_ford

Notes

There are no checks within this kernel.
Ax is assumed to be positive

Initializations:

- d[i] = 0 if i is a center, else inf
- m[i] = 0 .. num_clusters if in a cluster, else -1
- p[i] = -1

References

[1]

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

pyamg.amg_core.bellman_ford_balanced()#

Bellman-Ford with a heuristic to balance cluster sizes

This version is modified to break distance ties by assigning nodes to the cluster with the fewest points, while preserving cluster connectivity. This results in more balanced cluster sizes.

Parameters:

num_nodes(IN) number of nodes (number of rows in A)
Ap[](IN) CSR row pointer for A (num_nodes x 1)
Aj[](IN) CSR column index for A (num_edges x 1)
Ax[](IN) CSR data array (edge weights) (num_edges x 1): c : (INOUT) cluster center (num_clusters x 1) d : (INOUT) distance to cluster center (num_nodes x 1) m : (INOUT) cluster index (num_nodes x 1) p : (INOUT) predecessor on shortest path to center (num_nodes x 1) pc : (INOUT) number of predecessors (num_nodes x 1) s : (INOUT) cluster size (num_clusters x 1)

See also

pyamg.graph.bellman_ford

Notes

There are no checks within this kernel.
Ax > 0 is assumed

pyamg.amg_core.block_approx_ideal_restriction_pass2()#

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

Parameters:

Rparray: Pre-determined row-pointer for R in CSR format.
Rjarray: Empty array for column indices for R in CSR format.
Rxarray: Empty array for data for R in CSR format.
Aparray: Row pointer for matrix A.
Ajarray: Column indices for matrix A.
Axarray: Data array for matrix A.
Cparray: Row pointer for SOC matrix, C.
Cjarray: Column indices for SOC matrix, C.
Cxarray: Data array for SOC matrix, C.
Cptsarray: List of global C-point indices.
splittingarray: 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:

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

Notes

Rx[] must be passed in initialized to zero.

pyamg.amg_core.block_gauss_seidel()#

Block Gauss-Seidel iteration.

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

Returns:

None: Result in place.

pyamg.amg_core.block_jacobi()#

Block Jacobi iteration.

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: 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.
blocksizeint: Dimension of square blocks in BSR matrix A.

Returns:

None: Result in place.

pyamg.amg_core.block_jacobi_indexed()#

Indexed Block Jacobi iteration.

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:

None: Array x will be modified in place.

pyamg.amg_core.breadth_first_search()#

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

Returns:

None: In place.

Notes

The values of the level must be initialized to -1.

pyamg.amg_core.bsr_gauss_seidel()#

Gauss-Seidel iteration with BSR arrays.

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

None: Array x will be modified inplace.

pyamg.amg_core.bsr_jacobi()#

Weighted Jacobi iteration on BSR arrays.

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: Approximate solution.
barray: Right hand side.
temparray: 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:

None: Array x will be modified inplace.

pyamg.amg_core.bsr_jacobi_indexed()#

Indexed weighted Jacobi on BSR arrays.

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:

None: Array 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.center_nodes()#

Update center nodes for a cluster

Parameters:

num_nodes(IN) number of nodes (number of rows in A): Ap[] : (IN) CSR row pointer for A (num_nodes x 1) Aj[] : (IN) CSR column index for A (num_edges x 1) Ax[] : (IN) CSR data array (edge weights) (num_edges x 1)
Cptr[](INOUT) ptr to start of indices in C for each cluster (num_clusters x 1): D[] : (INOUT) FW distance array (max_size x max_size) P[] : (INOUT) FW predecessor array (max_size x max_size) C[] : (INOUT) FW global index for current cluster (num_nodes x 1) L[] : (INOUT) FW local index for current cluster (num_nodes x 1) q : (INOUT) FW work array for D**2 (max_size x max_size) c : (INOUT) cluster center (num_clusters x 1) d : (INOUT) distance to cluster center (num_nodes x 1) m : (INOUT) cluster index (num_nodes x 1) p : (INOUT) predecessor on shortest path to center (num_nodes x 1) pc : (INOUT) predecessor count (num_nodes x 1) s : (INOUT) cluster size (num_clusters x 1)

Returns:

changedflag to indicate a change in arrays D or P

Notes

sort into clusters first O(n)
s: [4 2 4 ….

Cptr: [0 4 6 11 …

| | |

v v v v

C: [87 99 4 6 82 13 15 9 12 55 66 77 …]
^ ^ ^ | |________ |_____ |_____ | |

| |

L: [ 2 3 1
^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ | | | | | | | | | | | … 0 1 2 3 4 5 6 7 8 9 10 …

pass pointer to start of each C[start,…., start+N]
N is the cluster size

pyamg.amg_core.classical_strength_of_connection_abs()#

Classical strength of connection.

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:

|A[i,j]| >= theta * max( |A[i,k]| )   where k != i

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:

None: Array S is 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()#

Classical strength of connection.

A[i,j] >= theta * max( -A[i,k] )   where k != i

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:

None: Array S is 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.cljp_naive_splitting()#

Compute a CLJP splitting.

Parameters:

nint: Number of rows in A (number of vertices).
Sparray: CSR row pointer (strength matrix).
Sjarray: CSR index array.
Tparray: CSR row pointer (transpose of the strength matrix).
Tjarray: CSR index array.
splittingarray, inplace: Array to store the C/F splitting.
colorflagint: Flag to indicate coloring.

Returns:

None: In place.

Notes

The splitting array must be preallocated. CLJP naive since it requires the transpose.

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.

Returns:

None: In place.

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.

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:

None: 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.

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:

Sparray

Row pointer array for CSR matrix S.

Sjarray

Col index array for CSR matrix S.

Sxarray

Value array for CSR matrix S. Upon entry to the routine, S = (I - (t/k) Dinv A)^k.

nrowsint

Dimension of S.

Barray

Array of size (nrows, NullDim) of near nullspace vectors in col major form, if calling from within Python, take a transpose.

DBarray

Array of size (nrows, NullDim) 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.

barray

Array of size (nrows, BDBCols) 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]

BDBColsint

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

NullDimint

Number of nullspace vectors.

tolfloat

Used to determine when values are numerically zero.

Returns:

array: Modified inplace and holds strength values for the above minimization problem.

See also

evolution_strength_of_connection

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.

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

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

None: Array 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:

None: Nothing, Ax is modified in place.

pyamg.amg_core.fit_candidates()#

pyamg.amg_core.floyd_warshall()#

Floyd-Warshall on a subgraph or cluster of nodes in A.

Parameters:

num_nodesint: Number of nodes (number of rows in A).
Aparray: CSR row pointer for A, (num_nodes, 1).
Ajarray: CSR column index for A, (num_edges, 1).
Axarray: CSR data array (edge weights), (num_edges, 1).
Darray: FW distance array, (max_size, max_size).
Parray: FW predecessor array, (max_size, max_size).
Carray: FW global index for current cluster, (N, 1).
Larray: FW local index for current cluster, (num_nodes, 1).
marray: Cluster index, (num_nodes, 1).
aarray: Center of current cluster.
Nint: Size of current cluster.

Notes

There are no checks within this kernel
There is no initialization within this kernel
Ax > 0 is assumed
Only a slice of C is passed to Floyd–Warshall. See center_nodes.
C[i] is the global index of i for i=0, …, N in the current cluster
N = |C|
L = local indices, nx1 (-1 if not in the cluster)
assumes a fully connected (directed) graph

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.gauss_seidel()#

Gauss-Seidel iteration.

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

None: Array 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()#

Indexed Gauss-Seidel iteration.

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

None: Array 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()#

Gauss-Seidel NE iteration.

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:

None: Array x is modified inplace in an additive, not overwriting fashion.

pyamg.amg_core.gauss_seidel_nr()#

Gauss-Seidel NR iteration.

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:

None: Array x is modified inplace in an additive, not overwriting fashion.

pyamg.amg_core.householder_hornerscheme()#

Householder Horner Scheme.

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: Vector of length n to be operated on.
Barray: Matrix of size (n,m) 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:

None: 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 are gmres and fgmres.

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

pyamg.amg_core.incomplete_mat_mult_bsr()#

Mat-mul over a sparsity pattern.

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:

Aparray: BSR row pointer array.
Ajarray: BSR col index array.
Axarray: BSR value array.
Bparray: BSR row pointer array.
Bjarray: BSR col index array.
Bxarray: BSR value array.
Sparray: BSR row pointer array.
Sjarray: BSR col index array.
Sxarray: BSR value array.
n_browint: Number of block-rows in A.
n_bcolint: Number of block-cols in S.
brow_Aint: Row blocksize for A.
bcol_Aint: Column blocksize for A.
bcol_Bint: 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 a pre-existing sparsity.

Use the 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:

Aparray: Row pointer array for CSR matrix A.
Ajarray: Col index array for CSR matrix A.
Axarray: Value array for CSR matrix A.
Bparray: Row pointer array for CSC matrix B.
Bjarray: Col index array for CSC matrix B.
Bxarray: Value array for CSC matrix B.
Sparray: Row pointer array for CSR matrix S.
Sjarray: Col index array for CSR matrix S.
Sxarray: Value array for CSR matrix S.
dimenint: Dimensionality of A,B and S.

Returns:

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
>>> import numpy as np
>>> from scipy.sparse import csr_array, csc_array
>>> A = csr_array(np.arange(1,10,dtype=float).reshape(3,3))
>>> B = csc_array(np.ones((3,3),dtype=float))
>>> AB = csr_array(np.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()#

Weighted Jacobi iteration.

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: Approximate solution.
barray: Right hand side.
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.

Returns:

None: Array x will be modified inplace.

pyamg.amg_core.jacobi_indexed()#

Indexed weighted Jacobi iteration.

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:

Aparray: CSR row pointer.
Ajarray: CSR index array.
Axarray: CSR data array.
xarray: Approximate solution.
barray: Right hand side.
temparray: Temporary vector the same size as x.
indicesarray: List of row indices to perform Jacobi on, e.g. F-points.
omegaarray: Damping parameter.

Returns:

None: Array x will be modified in place.

pyamg.amg_core.jacobi_ne()#

Jacobi NE iteration.

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

None: x is modified inplace in an additive, not overwriting fashion.

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

Returns:

None: In place.

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()#

Parallel maximal independent set.

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:

int: 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()#

Serial maximal independent set.

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: Value used for active vertices.
Cfloat: Value used to mark non-MIS vertices.
Ffloat: Value used to mark MIS vertices.
xarray: State of each vertex.

Returns:

int: 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: Num_rows array.
Aparray: CSR row pointer.
Ajarray: CSR index array.
Axarray: CSR data array.

Returns:

None: Array x[i] will hold row i’s maximum magnitude entry and is modified in place.

pyamg.amg_core.min_blocks()#

Find the size of the smallest entry in each block.

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_blocksint: Number of blocks in matrix.
blocksizeint: Size of each block.
Sxarray: Block data structure of BSR matrix, S. Sx is (n_blocks x blocksize) in length.
Txarray: Modified in place for output.

Returns:

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_array, csr_array
>>> from pyamg.amg_core import min_blocks
>>> from numpy import zeros, array, ravel, round
>>> from numpy 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_array( (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_array((T, S.indices, S.indptr), shape=(3,3))
>>> print("Matrix after\n" + str(S2.todense()))

pyamg.amg_core.most_interior_nodes()#

Find the most interior nodes.

Parameters:

num_nodesint: Number of nodes (number of rows in A).
Aparray: CSR row pointer for adjacency matrix A, (num_nodes, 1).
Ajarray: CSR index array, (num_edges, 1).
Axarray: CSR data array (edge lengths), (num_edges, 1).
carray, num_cluster: Cluster centers, (num_clusters, 1).
darray, num_nodes: Distance to nearest seed, (num_nodes, 1).
marray, num_nodes: Cluster index for each node, (num_nodes, 1).
parray: Predecessor on shortest path to center, (num_nodes, 1).

Notes

There are no checks within this kernel.
Ax is assumed to be positive

References

[Bell2008]

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

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 and each F-point from its strongest connected C-neighbor.

Parameters:

Rparray: Pre-determined row-pointer for P in CSR format.
Rjarray: Empty array for column indices for P in CSR format.
Cparray: Row pointer for SOC matrix, C.
Cjarray: Column indices for SOC matrix, C.
Cxarray: Data array for SOC matrix, C.
splittingarray: Boolean array with 1 denoting C-points and 0 F-points.

Returns:

None: Nothing, Rj[] modified in place.

pyamg.amg_core.overlapping_schwarz_csr()#

Overlapping Schwarz iteration.

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: 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.
nsdomainsint: Number of subdomains.
nrowsint: Number of rows.
row_startint: Subdomain processing start index.
row_stopint: Subdomain processing stop index.
row_stepint: Subdomain processing step index.

Returns:

None: Array 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: Array of size (m, n, n), 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:

None: 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
>>> import numpy as np
>>> A = np.array([np.arange(1,5, dtype=float).reshape(2,2), np.ones((2,2),dtype=float)])
>>> Ac = A.copy()
>>> pinv_array(A, 2, 2, 'T')
>>> print "Multiplication By Inverse\n" + str(np.dot(A[0], Ac[0]))
>>> print "Multiplication by PseudoInverse\n" + str(np.dot(Ac[1], np.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(np.dot(A[0], Ac[0]))
>>> print "A holds the inverse of the transpose\n" + str(np.dot(A[0], Ac[0].T))

pyamg.amg_core.remove_strong_FF_connections()#

Remove strong F-to-F 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:

None: Array Sx is set to zero to eliminate connections.

pyamg.amg_core.rs_cf_splitting()#

Ruge-Stuben 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()#

Ruge-Stuben splitting pass 2.

Parameters:

n_nodesint: Number of rows in A.
Sparray: CSR row pointer array for SOC matrix.
Sjarray: CSR column index array for SOC matrix.
splittingarray, inplace: Array to store the C/F splitting.

Notes

The splitting array must be preallocated.

pyamg.amg_core.rs_classical_interpolation_pass1()#

First pass of classical AMG interpolation.

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:

None: Pp is modified in place.

pyamg.amg_core.rs_classical_interpolation_pass2()#

RS classical interpolation pass 2.

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:

None: Arrays 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.

Returns:

None: In place.

Notes

See page 479 of Multigrid.

pyamg.amg_core.rs_direct_interpolation_pass2()#

RS direct interpolation pass 2.

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.

Returns:

None: Arrays Pj and Px modified in place.

pyamg.amg_core.satisfy_constraints_helper()#

Helper routine for satisfy_constraints routine.

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:

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

See also

energy_prolongation_smoother

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_array((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.sor_gauss_seidel()#

SOR iteration.

Perform one iteration of SOR 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: 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).
omegafloat: Relaxation parameter for SOR.

Returns:

None: Array x will be modified inplace.

Notes

Nearly identical to gauss_seidel with a relaxation parameter omega.

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

abs(A[i,j]) >= theta * sqrt( abs(A[i,i]) * abs(A[j,j]) )

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.

Only the largest (in magnitude) k entries per row are left. Smaller entries are zeroed out.

Parameters:

n_rowint: Number of rows in A.
kint: Number of entries per row to keep.
Sparray: CSR row pointer.
Sjarray: CSR index array.
Sxarray: CSR data array.

Returns:

None: 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 parallel Largest-Degree-First (LDF).

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.
zarray: Initial random values for each vertex.

Notes

Arrays x and z 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.

pyamg.amg_core#

This Page