This is called a `BandedBlockBandedMatrix`

in BlockBandedMatrices.jl, which is the topic of the next blog post.

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

**dlfivefifty**#21

**rveltz**#22

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

**dlfivefifty**#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

**rveltz**#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`

?

**dlfivefifty**#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
```

**rveltz**#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`

.

**dlfivefifty**#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`

.

**rveltz**#28

the condition `A[Nx+1,Nx] == 0.`

is satisfied. I guess my question is how do I pass my bands to `BandedBlockBandedMatrix`

?

**dlfivefifty**#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
```

**rveltz**#32

I need to solve `A x = rhs`

β¦ I guess I can do as in your example and try Gauss Seidel