This package allows you to solve linear systems of the type M * x = b
where M
is block diagonal (sparse or not).
It is particularly efficient if some of the blocks of M
are repeated, because it will only compute the factorizations of these repeated objects once.
Consider the blockdiagonal matrix
M = [A ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ A ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ B ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ A ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ C ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ A ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ C ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ B ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ A]
Instead of creating that big matrix, factorizing it whole, or factorizing each block, you can create a BlockFactors
or SparseBlockFactors
object (depending if A
, B
, and C
are sparse) via the following syntax
# Form an array of the matrices
Ms = [A, B, C]
# and an array of "repetition" indices
indices = [1, 1, 2, 1, 3, 1, 3, 2, 1]
# to create the Block Diagonal Factors (BDF) object
BDF = factorize(Ms, indices)
This way A
, B
, and C
are factorized only once.
Then, you can solve for linear system M * x = b

via backslash (same as
M \ b
)x = BDF \ b

via the inplace (same as
ldiv!(M, b)
)ldiv!(BDF, b)

or via the inplace (same as
ldiv!(x, M, b)
)ldiv!(x, BDF, b)
The package simply creates two new types, BlockFactors
or SparseBlockFactors
, which look like
struct (Sparse)BlockFactors{Tv}
factors::Vector
indices::Vector{<:Int}
m::Int
n::Int
end
and overloads factorize
, lu
, and other factorization functions to create those objects from an array of matrices and the repeating indices.
In order to solve linear systems it also overloads \
and ldiv!
.
That's it!
If you use this package directly, please cite it! Use the CITATION.bib, which contains a bibtex entry for the software (coming soon).