a1 = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
a2 = np.array([[2, 2, 2], [2, 2, 2], [2, 2, 2]])
a3 = np.array([[3, 3, 3], [3, 3, 3], [3, 3, 3]])
import scipy.linalg
scipy.linalg.block_diag(a1, a2, a3)
r = np.array([[1, 1, 1, 0, 0, 0, 0, 0, 0], [1, 1, 1, 0, 0, 0, 0, 0, 0], [1,
    1, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 2, 2, 2, 0, 0, 0], [0, 0, 0, 2, 2, 2,
    0, 0, 0], [0, 0, 0, 2, 2, 2, 0, 0, 0], [0, 0, 0, 0, 0, 0, 3, 3, 3], [0,
    0, 0, 0, 0, 0, 3, 3, 3], [0, 0, 0, 0, 0, 0, 3, 3, 3]])
(scipy.linalg.block_diag(a1, a2, a3) == r).all()
How can I transform blocks into a blockdiagonal matrix (NumPy)