Skip to content

Distributed Regridding

Jump to bottom
Julia Sloan edited this page Jul 13, 2023 · 14 revisions

Sparse Matrix Representation in Julia (SparseMatrixCSC)

A simple example

Consider the following dense matrix:

$$dense\_matrix = \begin{bmatrix} 6 & 0 & 7 & 0\\0 & 0 & 1 & 0\\-8 & 3 & 0 & 0\\0 & 0 & 5 & -2 \end{bmatrix}$$

We can use Julia's sparse matrix representation to store the dense matrix as a collection of integers and vectors instead:

using SparseArrays
dense_matrix = [6 0 7 0; 0 0 1 0; -8 3 0 0; 0 0 5 -2]
sparse_matrix = SparseMatrixCSC(dense_matrix)

sparse_matrix.m: 4

sparse_matrix.n: 4

sparse_matrix.nzval: [-6 7 1 -8 3 5 -2], with length 7

sparse_matrix.rowval: [1 3 3 1 2 4 4] with length 7 (= |sparse_matrix.nzval|)

sparse_matrix.colptr: [1 3 4 7 8] with length 5 (= sparse_matrix.n + 1)

Here, m and n are the number of rows and columns in sparse_matrix, respectively.

sparse_matrix.nzval is a vector containing the nonzero values of dense_matrix.

sparse_matrix.rowval is a vector containing the row indices of the nonzero values of dense_matrix. That is, sparse_matrix.rowval[i] gives the row index of the value sparse_matrix.nzval[i] in the original dense_matrix.

sparse_matrix.colptr is a vector containing the indices into sparse_matrix.nzval of the first nonzero value in each column of the dense_matrix, as well as a final bound which is 1 larger than the largest index of sparse_matrix.nzval (i.e. |sparse_matrix.nzval| + 1). Values from column j in the original dense_matrix will have indices in sparse_matrix.nzval in the range [sparse_matrix.colptr[j], sparse_matrix.colptr[j+1]). Because this vector contains one index for each column of the original matrix and one extra bound value, its length is always n + 1. This is perhaps a less intuitive approach than storing all the column indices into the original dense_matrix would be (as is done for rowval), but it encodes sufficient information for the user to be able to recover the dense matrix from the sparse representation, and requires storing fewer numbers than storing all column indices would.

These fields provide all the necessary information to index into a sparse matrix, or to convert a sparse matrix to the corresponding dense matrix.

A more complicated example

Consider the case where we want to remap from a source space with 1 element in the vertical, 2 elements in the horizontal, and 3 nodes per element (nex = 1, ney = 2, nq = 3), to a target space with 1 element in the vertical, 3 elements in the horizontal, and 3 nodes per element (nex = 1, ney = 3, nq = 3). The source space has 2 * 3^2 = 18 nodes total, and the target space has 3 * 3^2 = 27 nodes total.

The weight matrix used to map between these spaces will accordingly have 18 columns (corresponding to the nodes of the source space) and 27 rows (corresponding to the nodes of the target space).

The weight matrix for the mapping can be represented as a sparse matrix, and will have the following properties:

weights.m = 27

weights.n = 18

length(weights.colptr) = 19 = weights.n + 1

length(weights.rowval) = length(weights.nzval) = 96

See also

Julia's SparseMatrixCSC documentation

Indexing Conventions

Element Indexing

Global Indices

  • Integer providing element ordering globally across all processes

Local Indices

  • Integer providing element ordering locally within each process
  • For a serial space, this is the same as the global index
  • Note that when working on a distributed space, each element is fully contained within one process. There is no overlap of elements across multiple processes.
  • The number of local indices across all processes is equal to the number of global indices for a distributed space.

Cartesian Indices

  • A tuple (elem_x, elem_y) that provides a unique indexing for each element

Nodal Indexing

Redundant Indexing

  • Constructed by traversing all nodes, keeping a counter of node number. For each node, assign it the node number and increment by 1 after each assignment
  • This results in multiple number assignments for nodes which exist at the boundary of more than one element, since these nodes are effectively duplicated for each element they appear within.

Unique Indexing

  • Constructed similarly to redundant indices, except that nodes which already have numbers are skipped
  • This results in a single unique index for each node, even those which appear in multiple elements

TODO add diagram showing indices

Conventions used by ClimaCore, TempestRemap

ClimaCore

  • Local element indices (usually denoted by lidx)
  • Redundant node indices
  • Typically applies space filling curve for element ordering

TempestRemap

  • Unique global indices - one integer per node
  • Does not provide information about element indexing
  • Uses ordering of input topology
Clone this wiki locally