Construct new BlockSkylineMatrix from BlockSkylineMatrix and UnitRanges

[johnomotani]

Take an existing BlockSkylineMatrix (or BlockBandedMatrix), and select a sub-matrix by passing UnitRanges for the rows and columns to select, returning a BlockSkylineMatrix (or BlockBandedMatrix).

[johnomotani]

I think I'm going to want a feature like this (actually possibly, eventually, even a more complicated version where I'd allow Vector{<:Int} indexing instead of the UnitRange...) to implement a matrix-solve algorithm based on a domain decomposition. Would be happy to just keep a function like this in a separate repo somewhere, but maybe it's generally useful?

Needs tests - I can probably work on some if there's interest in merging this PR.

Question: I've attempted to make the data copying efficient, but it's pretty ugly, so possibly fragile. Maybe an experienced BlockBandedMatrices.jl developer can see a better way?

[codecov]

Codecov Report

❌ Patch coverage is 0% with 57 lines in your changes missing coverage. Please review.
✅ Project coverage is 85.23%. Comparing base (dbb1e50) to head (df22d67).

Files with missing lines	Patch %	Lines
src/BlockSkylineMatrix.jl	0.00%	57 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #229      +/-   ##
==========================================
- Coverage   88.31%   85.23%   -3.09%     
==========================================
  Files          11       11              
  Lines        1113     1172      +59     
==========================================
+ Hits          983      999      +16     
- Misses        130      173      +43     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[dlfivefifty]

These should be Int, not Int64

[dlfivefifty]

This seems very questionable. I think you might want to be using blockaxes or blocklengths.

I'm not going to review the rest since it all seems overly complicated.

[dlfivefifty]

You probably want AbstractUnitRange

[dlfivefifty]

I'm not convinced this is the right name. The comment isn't exactly clear what this doing.

[dlfivefifty]

Adding tests would make it clear what exactly this is doing.

I think this seems overly complicated. I'm pretty sure you can accomplish the same thing in a couple of lines using blockaxes. Once I understand what you are doing I can make some more suggestions.

[johnomotani]

Thanks for the pointers @dlfivefifty. I've been thinking harder about what exactly my use case is - I wrote a script to mock up the essential features, and now I think I understand better exactly what I want to do. I think you're right that this PR as it stands is not well-conceived.

The original idea was that for some BlockSkylineMatrix bm, any slice bm[rowinds,colinds] also has the structure of a BlockSkylineMatrix, but as far as I can see there is no function to explicitly convert that slice into a BlockSkylineMatrix. The intention was to define an efficient function where BlockSkylineMatrix(bm, rowinds,colinds) returns a BlockSkylineMatrix constructed so that BlockSkylineMatrix(bm, rowinds,colinds) == bm[rowinds,colinds]. However, I'm no longer sure that that's the most useful thing to do.

My particular motivation is that I will have a BlockSkylineMatrix where the blocks are pretty large (maybe 1000x1000) but sparse (think of the overall matrix as something like a three-dimensional derivative stencil). It seems like too much work (and not necessarily efficient, as the indexing would get much more complicated) to have a specialised sparse matrix type with the exact structure of my matrix. A BlockSkylineMatrix where I use the known structure in one 'outermost' dimension to identify large blocks that I know are zero seems like a good compromise with relatively simple indexing but also giving a large reduction in memory usage compared to a dense Matrix. I then want to divide the matrix into sub-blocks according to some domain decomposition - I think in the end the details of how I do that are not relevant, the end state is that I want to work with some sub-matrix (that isn't just one or several blocks of the original matrix, because the blocks also get divided up). OK, so I can just get a view of the original matrix and work with that.

My problem then is that I will want to do something like convert that sub-matrix view to a SparseMatrixCSC (like I was talking about in issue #228). Assuming the sub-matrix is large enough that falling back on the AbstractArray version of sparse() is not acceptable, I want some optimised version. I suspect the best way of doing this is to iterate through the non-zero parts of the block-columns of the submatrix convert each of those (since they're dense, but possibly containing many zeros in my application) to a SparseMatrixCSC, and then append the SparseMatrixCSC for the block-column to the SparseMatrixCSC that's being built up for the full sub-matrix (this is what I was trying to do in #228, but working from a BlockSkylineMatrix, not a view/SubArray). So I think what I need to know, is how to:

1. get the block structure of a view, v, of a BlockSkylineMatrix. As far as I can see at the moment, blockaxes(v), blocksizes(v), blocklengths(v) all treat v as though it was a single block? Maybe related to Return a BandedMatrix from a view of a BandedBlockBandedMatrix #223?
2. get either a copy or a view of the data-vector for (the non-zero entries of) one block-column of the view v.
3. get the row/column indices of the first point in a block-column, and either the number of rows/columns or the indices of the last point.

I guess I should probably be able to figure out points 2 and 3 given a solution to 1.

Edited to add: the reason for wanting to work with block-columns instead of just iterating over non-zero entries is that I'm assuming the blocks are large, so creating a SparseMatrixCSC one entry at a time would be slow compared to working out the offsets, etc. for converting a block.

[dlfivefifty]

It sounds like you want a BlockedArray wrapping a SparseMatrixCSC ?

If you need to know the block bandwidths we can add functionality for computing it from the I and J

are you sure you don’t want to use block indexing like A[Block.(2:3) , Block.(3:4)] ? This automatically preserves the block structure.

[johnomotani]

Thanks for the time you've already spent @dlfivefifty, so apologies if the following ends up being too much detail... (Maybe see TL;DR at the bottom.)

The problem I'm working on is to construct a Jacobian matrix (to be used as a preconditioner), and then decomposing/factorizing it in a way that can be parallelised efficiently, using knowledge of the very special structure of this particular matrix. The vector x in M*x=b represents a variable on a three-dimensional grid (possibly 4d or 5d at some point in the future). I don't want to start with a SparseMatrixCSC (or a blocked collection of them) because inserting entries into the matrix is much easier to do with simple row/column indexing, which would be inefficient with a SparseMatrixCSC. The strategy for parallelising is to decompose the matrix by splitting recursively, until we reach a point where there is one matrix chunk per process, which will then be converted to a SparseMatrixCSC and LU-factorized.

The splitting is where things start to get complicated. I've attached a sketch where I try to show an example where we split the grid in half, by dividing each of the three dimensions (say x, y and z) in half in turn. [The things in the sketch are not all to scale, as I've made some parts bigger/smaller to make it readable]

• Splitting the first dimension, x, is relatively straightforward, and nicely block structured (going from top-left to top-right in my sketch). There is a set of points (all positions in y and z for the single x value half way along the grid) that divides the grid in two - we remove those points (which will be handled specially to capture the coupling, using a Schur complement method that we don't need to go into here) and the rest of the matrix splits into two decoupled blocks $A_1^1$ and $A_2^1$.
• Next we want to split the y-dimension (top-right to bottom-right in my sketch). Logically this is the same operation, but is not as compatible with the vector ordering. For each point in x, there is a block corresponding to the y-z grid. Each of these blocks gets split in half (in a similar way as x was split), but now, for example, the top-left block of the split is coupled to other top-left blocks at different values of x (the adjacent values of x). So conceptually, we reassemble the coupled blocks into partial matrices $A_1^2$, $A_2^2$, $B^2$, $C^2$, $D^2$ that are still block-structured but each block is (roughly) half the size on each side as the blocks in the pre-split matrix.
• At the third level, we want to split the z-dimension (indicated as the blow-up to the left of the $A_1^2$ matrix). Now we have a larger number of even smaller blocks, one for each position in x and y. Each of these small blocks we split in the same way, and again need to reassemble all the small top-left blocks into a single matrix (and similar for other blocks).

TL;DR

In the end, the result of all this decomposition is that for each of many small submatrices (that I then want to LU-factorize) I have set of row and column indices that select that sub-matrix out of the original matrix. I also know the structure of the non-zero elements in the sub-matrix. So maybe the best solution (and it seems like it should be simple to implement) is just to create a view representing the sub-matrix, and then loop over the non-zero elements (which in principle I know the positions of) to fill a SparseMatrixCSC to be LU-factorized. With that approach no code needs adding to this library, and in principle it should be as efficient as possible (although I have to possibly do a bit more work to keep track of the indices of the decomposed sub-matrices).

[johnomotani]

I think we can close this now, unless you have suggestions/questions?

Also, thanks for the discussion! This was very helpful for me.