Create submatrix for sparse matrix formats

[pratikvn]

This PR adds support for creating submatrices from spans for sparse matrix formats. Dense already has this functionality.

An additional interface is added to ensure the same function signatures for the different matrix formats.

Only CSR is supported for now, but other formats could be easily added.

TODO

• Merge Add reduction for arrays. #831

[codecov]

Codecov Report

Merging #885 (7494456) into develop (b4fc0c7) will increase coverage by 0.02%.
The diff coverage is 98.50%.

@@             Coverage Diff             @@
##           develop     #885      +/-   ##
===========================================
+ Coverage    94.80%   94.82%   +0.02%     
===========================================
  Files          443      443              
  Lines        36459    36592     +133     
===========================================
+ Hits         34565    34699     +134     
+ Misses        1894     1893       -1     

Impacted Files	Coverage Δ
core/device_hooks/common_kernels.inc.cpp	0.00% <0.00%> (ø)
include/ginkgo/core/matrix/csr.hpp	45.36% <ø> (ø)
core/matrix/csr.cpp	98.88% <100.00%> (+0.07%)	⬆️
omp/matrix/csr_kernels.cpp	94.81% <100.00%> (+0.28%)	⬆️
omp/test/matrix/csr_kernels.cpp	100.00% <100.00%> (ø)
reference/matrix/csr_kernels.cpp	99.57% <100.00%> (+0.02%)	⬆️
reference/test/matrix/csr_kernels.cpp	99.79% <100.00%> (+0.01%)	⬆️
omp/reorder/rcm_kernels.cpp	97.53% <0.00%> (-0.61%)	⬇️
omp/components/format_conversion.hpp	100.00% <0.00%> (ø)
core/base/extended_float.hpp	92.23% <0.00%> (+0.97%)	⬆️
... and 2 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b4fc0c7...7494456. Read the comment docs.

[upsj]

Looks pretty useful! I'm guessing this is related to RAS overlaps? Does this degree of genericity suffice, i.e. are ranges enough to describe overlapping regions, or can there be gaps?

The code itself also looks good, only minor comments and questions (I will review the array reduce part in the other PR). Nonzero filtering is also part of the ParILUT implementation, not sure though whether it makes sense to try and merge those two.

[pratikvn]

Yes, ideally, I think for Schwarz type algorithms you would need a superset of this functionality, which i think could probably take a set of ranges, for example an IndexSet and build the submatrix from that. But I think that is slightly more complicated to implement. I will look into that once we have the IndexSet merged.

Nevertheless, I think this simple one with just a continuous range is also useful.

[Slaedr]

Nice work. This seems like it might be a nice feature. Could you describe the use case you have in mind for this?

[upsj]

LGTM in both implementation and tests, really nice job there. The warp-parallel implementation would be nice to have, but not necessary in this PR from my side.
I only don't think creating a separate class requiring a CRTP parameter gives any advantage over directly adding the function to Csr though. Also, I would like to avoid the addition of another constructor to Csr, since we already have a way to do this with little more core.

[Slaedr]

I would still like to know what the use case is for this functionality. I think that context is important in deciding whether this is a useful feature or not. For example, I think this might be useful if the create_submatrix could work on more general IndexSets rather than spans.

[pratikvn]

Extraction of submatrices could be used in many algorithms:

1. A quick test for Schur complement type algorithms.
2. In Domain decomposition algorithms where you would need to extract blocks, possibly overlapping of a global matrix.
3. Maybe in other blocking algorithms (atleast for the generation stage), as this allows you to extract any sub-block of a matrix, given a row and a column span.

While IndexSet would be more general, the kernels added here would be used for creating sub matrices which take an IndexSet as well. Additionally, IndexSet sub-matrix creation might not be the one that has a lot of use cases, because that much generality might be unnecessary, but a simple span based sub-matrix extraction will be more useful.

[upsj]

On the domain decomposition algorithms, is it normally possible to represent blocks and overlaps as contiguous ranges? In a straightforward 5pt stencil implementation, you can only get slices that way, and if it was stored in a blocked fashion, the overlaps probably need to consist of many separate ranges except for a single contiguous one?

[Slaedr]

@pratikvn I guess you have fine-grain domain decomposition in mind, where you would have many subdomains in a single process? In that case, I suppose you plan to reorder the DOFs such that those within a subdomain are always contiguous?

[pratikvn]

In general, overlaps are just extensions of the diagonal block by a few elements along the diagonal in either direction and then seeing the entire matrix with the rows to the left/right and columns to the top/bottom of the overlap elements. You would like to try to get the entire matrix (including the overlap) as a single block as that improves performance due to just a single apply call.

@Slaedr, yes, both fine-grained and coarse grained should be possible. For fine grained implmentations, you probably would not use the core side implementations, but call the kernels themselves in your algorithms. But to compose the algorithms for testing and prototyping purposes, this would be nevertheless useful.

[Slaedr]

I missed @upsj 's point about overlaps. If you need to extract overlapping sub-blocks, it would not always be possible to have the subdomains contiguous in the global ordering. I guess that significantly decreases the current interface's utility for domain decomposition?

You might also consider having IndexSet in the interface, but implement it only for a single contiguous range for now, and expand it later if need be.

[pratikvn]

@Slaedr , yes, I plan to add some functions which take in a IndexSet, but will probably do it another PR, as that will take some more time.

[Slaedr]

Except for the question of IndexSet vs span, I think this looks good. I agree with you that this is relevant for all matrix formats.

In general, I think IndexSet makes sense for this. You could add a small convenience function to get a vector of spans from an index set, and possibly also to also return the first span of an index set.

Stale

[greole]

Some very minor things from my side.

[upsj]

LGTM! I have only a few naming/readability/type nits for all implementations.

[sonarqubecloud]

Kudos, SonarCloud Quality Gate passed!   

0 Bugs
0 Vulnerabilities
0 Security Hotspots
10 Code Smells

70.6% Coverage
9.3% Duplication

[Slaedr]

LGTM. I still don't quite know good uses for this, but I guess that's okay. It might at least come in handy for playing around (initial experiments) with stuff, and maybe some testing, as you said.

What would be interesting is creating sparse submatrix views. That could actually be used in production code for things like multi-physics problems, but that's a whole other story. We could have 4-array (-view) version of CSR format for that, for example.