AdjacencyList offset array can overflow

[chrisrichardson]

How to reproduce the bug

Running a problem of sufficient size on one process (e.g. Q2 problem on mesh with 60 million dofs) can overwhelm the range of AdjacencyList.

• When running e.g. on a GPU system, we may have many degrees of freedom on one process.
• If the number of non-zeros in the local sparsity pattern exceeds 2^31, the AdjacencyList offsets used when building the sparsity pattern will overflow.

Possible solution: use 64-bit indexing for offsets in AdjacencyList or template it.

Minimal Example (Python)

No response

Output (Python)

No response

Version

main branch

DOLFINx git commit

No response

Installation

No response

Additional information

No response

[jorgensd]

Oh dear.
I would say change it to int64, as I think we have introduce too many templates already (and I do not see a large benefit of this being templated).

[chrisrichardson]

I'm going to think about other solutions: we only really need it in a few places (sparsity pattern and MatrixCSR) so I don't think we want to make it int64 everywhere...

[jorgensd]

I'm going to think about other solutions: we only really need it in a few places (sparsity pattern and MatrixCSR) so I don't think we want to make it int64 everywhere...

Then I guess we want to template it, but maybe make some wrappers already at the C++ level to avoid everything being moved into header files.

[chrisrichardson]

I think I'm going to close this. On reflection, it is only useful when we have >50M dofs, and the associated matrix will occupy around 64Gb of memory, which is close to the limit, at the moment. We should however, have some kind of test to make sure we are not overflowing when creating sparsity patterns or matrices.

[garth-wells]

This should be fixed at some point. The simplest approach might be (1) avoid AdjacencyList when we don't really need it (see #2632), and then (2) store the offsets in a std::vector<std::size_t>.