Add function to remove (explicit) zero entries from a CSR matrix

[IgorBaratta]

Describe new/missing feature

Explicit zero values can be added to the entries of CSR matrix.
For some types of elements, the element tensor is not dense, but the sparsity pattern doesn't have this information.

Example using scipy:

V = dolfinx.fem.FunctionSpace(mesh, element)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.inner(u, v) * ufl.dx
a = dolfinx.fem.form(a, jit_options=jit_options)
A = dolfinx.fem.assemble_matrix(a)


# Convert to scipy sparse matrix
from scipy.sparse import csr_matrix
D = csr_matrix((A.data, A.indices, A.indptr))


print(D.shape, D.nnz, D.nnz/(D.shape[0]*D.shape[1]) * 100, "%")

D.data[np.abs(D.data) < 1e-15] = 0
D.eliminate_zeros()
D.prune()

print(D.shape, D.nnz, D.nnz/(D.shape[0]*D.shape[1]) * 100, "%")

Output example:

# Before
(64, 64) 512 12.5 %
# After
(64, 64) 256 6.25 %

Maybe a similar function should be implemented for MatrixCSR.h.

Suggestion user interface

class MatrixCSR
{
...
...

void eliminate_zeros(Scalar tol)
}

[chrisrichardson]

This will require MatrixCSR to recompute its sparsity and some of its internal data (insert position for finalize for example).
We also should check (in debug mode) that insertion to indices which are not in the sparsity cause appropriate exceptions.

[IgorBaratta]

This will require MatrixCSR to recompute its sparsity and some of its internal data (insert position for finalize for example).

Yes, I think that's the main idea of the issue.
Maybe instead of changing the matrix in place, we could create a new matrix with a new sparsity pattern.

We also should check (in debug mode) that insertion to indices which are not in the sparsity cause appropriate exceptions.

This something that we want to do anyway, right? Even without removing zeros, since we call std::lower_bound.

[garth-wells]

PETSc now has a function that does this for its matrix types, https://petsc.org/release/manualpages/Mat/MatFilter/.