# Initialise 2x2 block-diagonal matrix
**URL:** https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369
**Category:** New to Julia
**Tags:** arrays
**Created:** [June 4, 2021, 12:16am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369 "2021-06-04T00:16:39Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@baptnz](https://discourse.julialang.org/u/baptnz)
#### Post date: [June 4, 2021, 12:16am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/1 "2021-06-04T00:16:40Z")
I’m writing code using some linear algebra, where I need to set up a 6Nx6N matrix (called F), with a 2x2 block layout, each 3Nx3N block being itself composed of N^2 3x3 sub-blocks.
The initialisation of F should be such the diagonal of each 3Nx3N block is filled with ones.
Here’s the code I’m using:
```julia
T = Float64
N = 2
using LinearAlgebra
F = Matrix{T}(I, 6N, 6N) # T is inferred from the type of another variable
# main diagonal is filled with ones, good
# how do I do the same with the 2 off-diagonal blocks, without pre-allocating dummy identity blocks?
using StaticArrays
dummyI = SMatrix{3,3}(I)
for i=1:N
F[3N .+ (3i-2:3i), 3i-2:3i] = dummyI
F[3i-2:3i, 3N .+ (3i-2:3i)] = dummyI
end
# ... later on
# I do something with this matrix, filling non-trivial elements etc.
for i in 1:N
for j in i+1:N
# F[3i-2:3i, 3j-2:3j] = ... fill specific off-diagonal 3x3 subblocks
end
end
```
I wonder what would be the best way to initialise F without the dummy 3x3 identity block and the for loop to assign those to the two off-diagonal blocks.
BlockArrays seems like a logical option but I’m not too sure how to avoid temporary intermediates either (mortar((A,B),(C,D)) would be even worse, if I need to expand the blocks A,B,C,D).
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 4, 2021, 2:29am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/2 "2021-06-04T02:29:56Z")
> [@baptnz](#):
>
> I’m writing code using some linear algebra, where I need to set up a 6Nx6N matrix (called F), with a 2x2 block layout, each 3Nx3N block being itself composed of N^2 3x3 sub-blocks.
> The initialisation of F should be such the diagonal of each 3Nx3N block is filled with ones.
Sounds like you just want:
```julia
I₃ = Matrix{T}(I, 3N, 3N)
F = [I₃ I₃; I₃ I₃]
```
or equivalently
```julia
F = diagm(0 => ones(T,6N), 3N => ones(T,3N), -3N => ones(T,3N))
```
or equivalently
```julia
F = kron(ones(T,2,2), I(3N))
```
All of these match the `F` constructed by your code above.
---
### Author:  [@baptnz](https://discourse.julialang.org/u/baptnz)
#### Post date: [June 4, 2021, 2:46am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/3 "2021-06-04T02:46:45Z")
Thanks for this; I was hoping for a way to avoid expanding the 3Nx3N matrix as in the first solution (N could be in the thousands) – I am guessing the second and third solutions do just that? Both seem pretty neat.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 5, 2021, 12:12am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/4 "2021-06-05T00:12:53Z")
The `diagm` solution will be the most efficient if that is a concern.
Of course, if your matrix is that big, you should think about using some kind of sparse-matrix structure.
---
### Author:  [@baptnz](https://discourse.julialang.org/u/baptnz)
#### Post date: [June 5, 2021, 12:19am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/5 "2021-06-05T00:19:09Z")
Thanks again; I am curious why `diagm` will be more efficient than `kron` – is the latter creating temporary objects?
I don’t think I would gain from using sparse matrices since the block-identity structure is only for the initialisation; later in the code it becomes a fully-dense matrix ([electromagnetic coupled dipole interaction matrix](http://nano-optics.ac.nz/CoupledDipole.jl/dev/theory/#Coupled-dipole-equations) – I imagine that’ll speak to you ; )
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 5, 2021, 12:37am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/6 "2021-06-05T00:37:53Z")
> [@baptnz](#):
>
> Thanks again; I am curious why `diagm` will be more efficient than `kron` – is the latter creating temporary objects?
Because `diagm` is just populating the matrix, whereas `kron` has to do computations.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 5, 2021, 12:40am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/7 "2021-06-05T00:40:36Z")
> [@baptnz](#):
>
> I don’t think I would gain from using sparse matrices since the block-identity structure is only for the initialisation; later in the code it becomes a fully-dense matrix ([electromagnetic coupled dipole interaction matrix](http://nano-optics.ac.nz/CoupledDipole.jl/dev/theory/#Coupled-dipole-equations) – I imagine that’ll speak to you ; )
In that case (DDA, a form of VIE) it’s more a matter of using matrix-free algorithms like fast-multipole… There has been some recent nice work on fast circulant preconditioners for this ([[1903.07884] Circulant preconditioning in the volume integral equation method for silicon photonics](https://arxiv.org/abs/1903.07884))
---
### Author:  [@baptnz](https://discourse.julialang.org/u/baptnz)
#### Post date: [June 5, 2021, 12:49am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/8 "2021-06-05T00:49:54Z")
DDA is different in that it discretises bodies in a regular cubic lattice of point dipoles; in CDA (same equations, naturally – [Markel has good recent reviews on it](http://whale.med.upenn.edu/vmarkel/EPUBS/JQSRT-2019-236-106611_AccMan.pdf)) the dipoles can be placed (and oriented) arbitrarily, and typically one dipole represents a particle (in the Rayleigh-Gans regime).
The block-Toeplitz form of the DDA interaction matrix arises from the periodicity in dipole positions, but it’s not the case in CDA.
---
### Author:  [@baptnz](https://discourse.julialang.org/u/baptnz)
#### Post date: [June 5, 2021, 12:52am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/9 "2021-06-05T00:52:06Z")
Makes sense thanks, I guess I should look at the source code to see if those calculations are likely time-consuming. This was more for general knowledge though; in my use-case I only create the matrix once and it’s definitely not a bottleneck – just trying to build good habits in julia, and learn more about the language at the same time.
---
### Author:  [@baptnz](https://discourse.julialang.org/u/baptnz)
#### Post date: [June 5, 2021, 1:02am UTC](https://discourse.julialang.org/t/initialise-2x2-block-diagonal-matrix/62369/10 "2021-06-05T01:02:59Z")
I’ve been looking for a good iterative solver for these coupled dipole equations (no Toeplitz structure, so circulant won’t work), especially in the case where the right-hand side has multiple columns (can be in the hundreds for good orientation-averaging of cross-sections by spherical cubature. Haven’t been successful so far ([this seems promising though](https://juliahub.com/docs/KrylovMethods/qy696/0.6.0/autodocs/#KrylovMethods.blockBiCGSTB-Union%7BTuple%7BT%7D,%20Tuple%7BFunction,Array%7BT,N%7D%20where%20N%7D%7D%20where%20T)), so I use brute-force solve.