BlockDiagonalFactors

(coming soon!)

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