[LinAlg] Static sparsity layout for sparse matrices

[maxfirmbach]

Description

Currently the sparsity pattern used for our sparse matrix implementation with Epetra_CrsMatrix can be changed dynamically by setting the option StaticProfile = false. The effects are the following: we specify an approximate upper bound as number of non-zeros per row. This number is not a "hard" bound it can be considered as a hint. Afterwards we are still able to add more values in a row if necessary. Usually this will still throw a warning, if overused (e.g. one can specify zero nnz per row and we would still be able to insert values).

There are some drawbacks to this however:

• I would like to have a hard assert for warnings and errors from Epetra, which currently fails for insert_into_global_values() as we on a regular basis insert values without having the allocation ready, resulting in a lot of re-allocation during value insertion.
• The mechanism of un_complete called on a matrix does not help. We pre-allocate the memory and complete, but afterwards modify everything again and usually insert more values than allocated.

=> Violating allocation limits and dynamically changing memory slows down operations.

It makes sense to be cautious about when and where to allocate memory for a sparse matrix. Now to the main topic though. In Tpetra dynamic memory allocation is deprecated for matrices (mainly due to GPU computing). Therefore we need to switch to a static memory profile StaticProfile = true for our matrix implementation to make a later transition to Tpetra as smooth as possible.

I assume this will be a lot of work and code will break in several places.

Possible Solution and Definition of Done

• Turn the StaticProfile option in the Epetra_CrsMatrix constructor to true in all places.
• Run a complete ctest routine and document the failing tests.
• Get all tests running again by either reworking how many entries are pre-allocated for a specific matrix and/or search for a upper bound.

[maxfirmbach]

After a first initial test it seems like ALL tests do fail. We might need to rethink how we want to tackle this issue properly.

[mayrmt]

@maxfirmbach if all tests are failing, that could also be just one problem in the matrix class -- at least I hope so.

Do you know, if the error can be traced back to a single spot in our matrix implementation or if we really need to work on every construction of a matrix throughout the entire code base?

[maxfirmbach]

I'm 90% convinced that we have to touch each sparse matrix constructor throughout the code base. Our matrix class (from a constructor perspective) only hands over values (in most cases).

[maxfirmbach]

I think the initial problem is that we take for example "magic numbers" a little too easy, while handling the upper bound for the number of non-zeros per row. Plain numbers like "81" might work and this is of course highly dependent on the physics, but a more general concept would be better.

For example:
upper bound for nnz per row = (stencil size) * (DOFs per node)

But this of course gets tricky for coupling matrices etc.

[mayrmt]

The magic number 81 is based on exactly such an estimate, in particular for 3D solid mechanics with a pure HEX8 mesh. Similarly, fluid uses 108. But yes, these are just estimates based on very specific assumptions, that often do not hold in practice.

[maxfirmbach]

Yes, there is also nothing wrong about this in the first place (besides copy and pasta of these numbers to different inappropriate places). I think that we should still be a bit more delicate with it and also more explicit, plus, we really need a concrete upper bound that holds true for a whole simulation run and not just an estimate.