How do I efficiently construct a block matrix of the following form in JULIA?

4
votes

So I am working on a code that requires me to construct a large Matrix M using smaller 'square' matrices J AND M of sizes n x n each, repeated such that:

This is what the matrix should look like!

i.e with the dimensions of M such that M is repeated 'L' times along the diagonal, J' is repeated along the upper second diagonal and J on the lower second diagonal.

Note that I am working on Julia v 1.0.0 and as far as I understand there are no direct ways of assigning Block Matrices in Julia, unlike Mathematica.

I tried to use Kronecker products to solve my problem:

 ????=Diagonal(ones(L)) #IDENTITY matrix of L x L size
 ????=kron(????,M)

Doing so, I can make a Block diagonal matrix M with small matrix M repeated along its diagonal. But now how do I place matrices J and J' along its second diagonals as required?

2
matrixjulia

2 Answers

7
votes

BlockArrays.jl (and maybe BlockBandedMatrices.jl) should be what you are looking for as it makes handling block matrix structures very convenient.

An example (with Stringss):

julia> using BlockArrays

julia> M = fill("M", 2,2);

julia> J = fill("J", 2,2);

julia> B = BlockArray{String}(undef_blocks, [2,2], [2,2])
4×4 BlockArray{String,2,Array{String,2}}:
 #undef  #undef  │  #undef  #undef
 #undef  #undef  │  #undef  #undef
 ────────────────┼────────────────
 #undef  #undef  │  #undef  #undef
 #undef  #undef  │  #undef  #undef

julia> setblock!(B, M, 1,1);

julia> setblock!(B, M, 2,2);

julia> setblock!(B, J, 1,2);

julia> setblock!(B, J, 2,1);

julia> B
4×4 BlockArray{String,2,Array{String,2}}:
 "M"  "M"  │  "J"  "J"
 "M"  "M"  │  "J"  "J"
 ──────────┼──────────
 "J"  "J"  │  "M"  "M"
 "J"  "J"  │  "M"  "M"

For more information, check out the package's documentation.

0
votes

You can also do this without any external packages like this:

# assign stand-in values
L = 4;
n = 2;
M = reshape(1:n^2, (n, n));
J = -M;

# put together the larger matrix
μ = [zeros(Int64, n, n) for i in 1:L, j in 1:L];
μ[diagind(μ)] .= [M];
μ[diagind(μ, -1)] .= [J];
μ[diagind(μ, 1)] .= [J'];
bM = reduce(vcat, [reduce(hcat, μ[i, :]) for i in 1:L]);

The final matrix bM is then:

8×8 Array{Int64,2}:
  1   3  -1  -2   0   0   0   0
  2   4  -3  -4   0   0   0   0
 -1  -3   1   3  -1  -2   0   0
 -2  -4   2   4  -3  -4   0   0
  0   0  -1  -3   1   3  -1  -2
  0   0  -2  -4   2   4  -3  -4
  0   0   0   0  -1  -3   1   3
  0   0   0   0  -2  -4   2   4

Not sure how the effiency compares to other answers, but it may at least be nice to know it can be done without too much complication.