What's the best way to build a block sparse matrix whose entries are diagonal matrices?

I want to construct a square block matrix out of square matrices, whose entries all lie on their respective diagonals. For example,

clear
N = 3;
M = 5;
MFull = zeros(N * M, N * M);
    
for n = 1 : M
    for m = 1 : M
        Mblock = rand(N,1) * m * n;
        MFull((n - 1) * N + 1:n * N, (m - 1) * N + 1:m * N) = diag(Mblock);
    end
end

When this matrix becomes big, this is clearly inefficient. I would like to define MFull as a sparse matrix to avoid memory issues and speed things up. Any suggestions are appreciated! To clarify: I would like to avoid creating MFull as a full matrix, and then converting it to a sparse matrix.

Answer 1

[David Cyncynates]

Here's a solution that seems to work for me:

clear
N = 3;
M = 5;
BlockContainer = zeros(N * M, 2 * M - 1);
    
for n = 1 : M
    for m = 1: M
        row = m;
        column = M + m - n;
        
        Mblock = rand(N,1) * m * n;
        BlockContainer((row - 1) * N + 1:(row) * N,column) = Mblock;
    end
end
MFull = spdiags(BlockContainer,-(M - 1) * N:N:(M - 1) * N, N * M, N * M)

[David Cyncynates]

I included a factor of m * n to help make the structure of the matrix easier to see.