#' Higher Order Singular Value Decomposition
#'
#' @param X multi-dimensional array (at least a matrix)
#' @param nu Number of Singula Vector per mode. Defaults to a complete HO-SVD.
#' @param eps tolerance for detecting linear dependence in columns of a matrix.
#'  Used for rank deduction and passed to \code{\link{qr}}.
#'
#' @export
HOSVD <- function(X, nu = NULL, eps = 1e-07) {
    if (!is.null(nu)) {
        stopifnot(all(nu <= dim(X)))
    }

    # Compute per mode singular vectors
    Us <- Map(function(i) {
        xx <- mcrossprod(X, , i)
        La.svd(xx, if (is.null(nu)) qr(xx, tol = eps)$rank else nu[i], 0)$u
    }, seq_along(dim(X)))
    # Compute Core tensor
    C <- mlm(X, Map(t, Us))

    list(C = C, Us = Us)
}

SVD <- function(A) .Call("C_svd", A)