def block_diag_multiply(result, value, blocks, block_cumul_dim_0,
    block_cumul_dim_1):
    """result = result + block_diag(blocks)@value"""
    for i in range(len(blocks)):
        start_0 = block_cumul_dim_0[i]
        end_0 = block_cumul_dim_0[i + 1]
        start_1 = block_cumul_dim_1[i]
        end_1 = block_cumul_dim_1[i + 1]
        result[start_0:end_0] += blocks[i] @ value[start_1:end_1]
