Add implementation of the submatrix method

The submatrix method allows to approximately solve matrix functions for
large, sparse matrices by solving the same function on a set of smaller,
dense matrices. While the original publication mainly targets inverse
p-th roots, the method can similarly applied to other matrix functions.

This implementation of the method takes a dbcsr matrix as input and
assigns a set of submatrices to each of the ranks. Non-zero blocks of
the original matrix are then communicated to ranks that need those
blocks to generate their submatrices.

Operations on the submatrices themselves can be performed locally
without any additional communication. They can be trivially
parallelized. The mapping of submatrices to ranks happens in a way that
balances load between the ranks, assuming that operations performed on
the submatrices take O(n^3) time.

Results from all ranks are finally communicated back to generate a dbcsr
result matrix.

Original publication of the method:

Michael Lass, Stephan Mohr, Hendrik Wiebeler, Thomas D. Kühne, and
Christian Plessl. 2018. A Massively Parallel Algorithm for the
Approximate Calculation of Inverse p-th Roots of Large Sparse Matrices.
In Proceedings of the Platform for Advanced Scientific Computing
Conference (PASC ’18). Association for Computing Machinery, New York,
NY, USA, Article 7, 1–11. DOI:https://doi.org/10.1145/3218176.3218231