Skip to content

Conversation

johnomotani
Copy link

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

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
Copy link
Author

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?

Copy link

codecov bot commented Oct 13, 2025

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.

bi_start = findblockindex.(axes(A), (rows.start, cols.start))
bi_stop = findblockindex.(axes(A), (rows.stop, cols.stop))

first_block = Int64.(BlockArrays.block.(bi_start))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These should be Int, not Int64

last_blockindices = Int64.(BlockArrays.blockindex.(bi_stop))

orig_blocksizes = blocksizes(A)
orig_row_sizes = [bs[1] for bs view(orig_blocksizes, :, 1)]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

not guaranteed that the `data` for the result is a contiguous subset of `A.data`.
"""
function BlockSkylineMatrix(A::BlockSkylineMatrix{T,DATA,BS}, rows::UnitRange,
cols::UnitRange) where {T,DATA,BS}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You probably want AbstractUnitRange

`rows` and `columns`. This function allocates new memory and copies data, because it is
not guaranteed that the `data` for the result is a contiguous subset of `A.data`.
"""
function BlockSkylineMatrix(A::BlockSkylineMatrix{T,DATA,BS}, rows::UnitRange,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@dlfivefifty
Copy link
Member

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
Copy link
Author

johnomotani commented Oct 14, 2025

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.

@johnomotani johnomotani marked this pull request as draft October 14, 2025 14:50
@dlfivefifty
Copy link
Member

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
Copy link
Author

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

matrix-splitting

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
Copy link
Author

johnomotani commented Oct 15, 2025

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

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants