-
Notifications
You must be signed in to change notification settings - Fork 16
Distributed Regridding
Consider the following dense matrix:
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.
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