#' Multi Kronecker Multiplication
#'
#' \deqn{C = A (B_1\otimes ...\otimes B_r)}{%
#'       C = A (B_1 %x% ... %x% B_r)}
#'
#' @examples
#' n <- 17
#' p <- c(2, 7, 11)
#' q <- c(3, 5, 13)
#'
#' A <- matrix(rnorm(n * prod(p)), n)
#' Bs <- Map(matrix, Map(rnorm, p * q), p)
#'
#' stopifnot(all.equal(
#'     A %*% Reduce(`%x%`, Bs),
#'     mkm(A, Bs)
#' ))
#'
#' @export
mkm <- function(A, Bs) {
    # reshape
    dim(A) <- c(nrow(A), rev(mapply(nrow, Bs)))

    # perform equiv Multi Linear Multiplication
    C <- mlm(A, rev(Bs), seq_along(Bs) + 1, transposed = TRUE)

    # reshape back
    dim(C) <- c(nrow(C), prod(dim(C)[-1]))

    C
}