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

Specific Domains Numerics
21

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

1 Like
22

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

23

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

24

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?

25

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
26

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.

27

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.

28

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

29

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
30

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

31

Huh?

32

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