matmult_block_diagonal function

Multiply Block Diagonal Matrices in Parallel