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

#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?

0 Likes

#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

0 Likes

#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`?

0 Likes

#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
``````
0 Likes

#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`.

0 Likes

#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`.

0 Likes

#28

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

0 Likes

#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
``````
0 Likes

#30

Thank you a lot! Unfortunately, I had not seen errorβ¦

0 Likes

#31

Huh?

0 Likes

#32

I need to solve `A x = rhs`β¦ I guess I can do as in your example and try Gauss Seidel

0 Likes