New blog post on banded matrices and ordinary differential equations using BandedMatrices.jl

This is called a BandedBlockBandedMatrix in BlockBandedMatrices.jl, which is the topic of the next blog post.

1 Like

Dump question: by using a change of variable, can I bring the far band closer to the diagonal?

If it comes from finite differences or similar, you can’t. But there’s Cuthill-McKee algorithm that can calculate the minimi bandwidth permutation for general sparse matrices, which is yet to be implemented: https://github.com/JuliaMatrices/BandedMatrices.jl/issues/76

I tried to use your structure BandedBlockBandedMatrix for the 2d problem (see 2 answers above) I stated above but I could not do it. Is there a simplified call where you give the bands (which can be far appart) and it automatically asign it to the good blocks in BandedBlockBandedMatrix?

I thought you are just doing the FD matrix? That’s pretty easy:

using BlockBandedMatrices, BandedMatrices, LazyArrays, FillArrays

function finitedifference_2d(n)
    h = 1/n
    DΒ² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))/h^2
    D_xx = BandedBlockBandedMatrix(Kron(DΒ², Eye(n)))
    D_yy = BandedBlockBandedMatrix(Kron(Eye(n), DΒ²))
    D_xx + D_yy
end

No, it is not FD and it is not a kron. The band are 0,1,-1,Nx,-Nx where the matrix line (column) size is Nx * Nx. I know the vecteur for each band, say v0,v1,vm1,vnx,vmnx.

Oh that won’t work in BandedBlockBandedMatrix unless, for example, A[Nx+1,Nx] == 0. But since it looks like its an operator acting on a tensor grid, it would be very strange if A[Nx+1,Nx] β‰  0.

the condition A[Nx+1,Nx] == 0. is satisfied. I guess my question is how do I pass my bands to BandedBlockBandedMatrix?

I suspect it’s actually easier to generate block bands instead of bands (you are probably combining multiple block-bands into single bands at some point). So for example the following creates FD but you can change your data as appropriate:

julia> Nx = 3; A = BandedBlockBandedMatrix(Zeros(Nx^2,Nx^2), (Fill(Nx,Nx), Fill(Nx,Nx)), (1,1), (1,1));

julia> for K = 1:Nx
          view(A,Block(K,K))[band(1)] .= 1
          view(A,Block(K,K))[band(0)] .= -4
          view(A,Block(K,K))[band(-1)] .= 1
       end

julia> for K = 1:Nx-1
          view(A,Block(K,K+1))[band(0)] .= 1
          view(A,Block(K+1,K))[band(0)] .= 1   
       end

julia> A
9Γ—9 BandedBlockBandedMatrix{Float64,BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2}}}:
 -4.0   1.0    β‹…   β”‚   1.0   0.0    β‹…   β”‚    β‹…     β‹…     β‹… 
  1.0  -4.0   1.0  β”‚   0.0   1.0   0.0  β”‚    β‹…     β‹…     β‹… 
   β‹…    1.0  -4.0  β”‚    β‹…    0.0   1.0  β”‚    β‹…     β‹…     β‹… 
 ──────────────────┼────────────────────┼──────────────────
  1.0   0.0    β‹…   β”‚  -4.0   1.0    β‹…   β”‚   1.0   0.0    β‹… 
  0.0   1.0   0.0  β”‚   1.0  -4.0   1.0  β”‚   0.0   1.0   0.0
   β‹…    0.0   1.0  β”‚    β‹…    1.0  -4.0  β”‚    β‹…    0.0   1.0
 ──────────────────┼────────────────────┼──────────────────
   β‹…     β‹…     β‹…   β”‚   1.0   0.0    β‹…   β”‚  -4.0   1.0    β‹… 
   β‹…     β‹…     β‹…   β”‚   0.0   1.0   0.0  β”‚   1.0  -4.0   1.0
   β‹…     β‹…     β‹…   β”‚    β‹…    0.0   1.0  β”‚    β‹…    1.0  -4.0

Thank you a lot! Unfortunately, I had not seen error…

Huh?

I need to solve A x = rhs… I guess I can do as in your example and try Gauss Seidel