def diag_block_mat_slicing(L):
    shp = L[0].shape
N = len(L)
r = list(range(N))
out = np.zeros((N, shp[0], N, shp[1]), dtype=int)
out[(r), :, (r), :] = L
return out.reshape(np.asarray(shp) * N)
python numpy: how to construct a big diagonal array(matrix) from two small array