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 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).

@@            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.

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.

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.

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.