diff --git a/tensorPredictors/R/tensor_times_matrix.R b/tensorPredictors/R/tensor_times_matrix.R
deleted file mode 100644
index 6b7af26..0000000
--- a/tensorPredictors/R/tensor_times_matrix.R
+++ /dev/null
@@ -1,72 +0,0 @@
-#' Tensor Times Matrix (n-mode tensor matrix product)
-#'
-#' @param T array of order at least \code{mode}
-#' @param M matrix, the right hand side of the mode product such that
-#'  \code{ncol(M)} equals \code{dim(T)[mode]}
-#' @param mode the mode of the product in the range \code{1:length(dim(T))}
-#'
-#' @returns multi-dimensional array of the same order as \code{T} with the
-#'  \code{mode} dimension equal to \code{nrow(M)}
-#'
-#' @export
-ttm <- function(T, M, mode = length(dim(T))) {
-    mode <- as.integer(mode)
-    dims <- dim(T)
-
-    if (length(dims) < mode) {
-        stop(sprintf("Mode (%d) must be smaller equal the tensor order %d",
-                     mode, length(dims)))
-    }
-    if (dims[mode] != ncol(M)) {
-        stop(sprintf("Dim. missmatch, mode %d has dim %d but ncol is %d.",
-                     mode, dims[mode], ncol(M)))
-    }
-
-    # Special case of mode being equal to tensor order, then an alternative
-    # (but more efficient) version is Z M' where Z is only the reshaped but
-    # no permutation of elements is required (as in the case of mode 1)
-    if (mode == length(dims)) {
-        # Convert tensor to matrix (similar to matricization)
-        dim(T) <- c(prod(dims[-mode]), dims[mode])
-
-        # Equiv matrix product
-        C <- tcrossprod(T, M)
-
-        # Shape back to tensor
-        dim(C) <- c(dims[-mode], nrow(M))
-
-        C
-    } else {
-        # Matricize tensor T
-        if (mode != 1L) {
-            perm <- c(mode, seq_along(dims)[-mode])
-            T <- aperm(T, perm)
-        }
-        dim(T) <- c(dims[mode], prod(dims[-mode]))
-
-        # Perform equivalent matrix multiplication
-        C <- M %*% T
-
-        # Reshape and rearrange matricized result back to a tensor
-        dim(C) <- c(nrow(M), dims[-mode])
-        if (mode == 1L) {
-            C
-        } else {
-            aperm(C, order(perm))
-        }
-    }
-
-}
-
-#' @rdname ttm
-#' @export
-`%x_1%` <- function(T, M) ttm(T, M, 1L)
-#' @rdname ttm
-#' @export
-`%x_2%` <- function(T, M) ttm(T, M, 2L)
-#' @rdname ttm
-#' @export
-`%x_3%` <- function(T, M) ttm(T, M, 3L)
-#' @rdname ttm
-#' @export
-`%x_4%` <- function(T, M) ttm(T, M, 4L)
diff --git a/tensorPredictors/R/ttm.R b/tensorPredictors/R/ttm.R
new file mode 100644
index 0000000..5b3c8a4
--- /dev/null
+++ b/tensorPredictors/R/ttm.R
@@ -0,0 +1,28 @@
+#' Tensor Times Matrix (n-mode tensor matrix product)
+#'
+#' @param T array of order at least \code{mode}
+#' @param M matrix, the right hand side of the mode product such that
+#'  \code{ncol(M)} equals \code{dim(T)[mode]}
+#' @param mode the mode of the product in the range \code{1:length(dim(T))}
+#'
+#' @returns multi-dimensional array of the same order as \code{T} with
+#'  \code{mode} dimension equal to \code{nrow(M)}
+#'
+#' @export
+ttm <- function(T, M, mode = length(dim(T))) {
+    storage.mode(T) <- storage.mode(M) <- "double"
+    .Call("C_ttm", T, M, as.integer(mode))
+}
+
+#' @rdname ttm
+#' @export
+`%x_1%` <- function(T, M) ttm(T, M, 1L)
+#' @rdname ttm
+#' @export
+`%x_2%` <- function(T, M) ttm(T, M, 2L)
+#' @rdname ttm
+#' @export
+`%x_3%` <- function(T, M) ttm(T, M, 3L)
+#' @rdname ttm
+#' @export
+`%x_4%` <- function(T, M) ttm(T, M, 4L)
diff --git a/tensorPredictors/src/init.c b/tensorPredictors/src/init.c
new file mode 100644
index 0000000..5640700
--- /dev/null
+++ b/tensorPredictors/src/init.c
@@ -0,0 +1,23 @@
+#include <R.h>
+#include <Rinternals.h>
+#include <R_ext/Rdynload.h>
+
+// extern SEXP FastPOI_C_sub(SEXP in_B, SEXP in_Delta,
+//     SEXP in_lambda, SEXP in_maxit, SEXP in_tol
+// );
+
+/* Tensor Times Matrix a.k.a. Mode Product */
+extern SEXP ttm(SEXP A, SEXP X, SEXP mode);
+
+/* List of registered routines (e.g. C entry points) */
+static const R_CallMethodDef CallEntries[] = {
+    // {"FastPOI_C_sub", (DL_FUNC) &FastPOI_C_sub, 5},  // NOT USED
+    {"C_ttm", (DL_FUNC) &ttm, 3},
+    {NULL, NULL, 0}
+};
+
+/* Restrict C entry points to registered routines. */
+void R_init_tensorPredictors(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}
diff --git a/tensorPredictors/src/ttm.c b/tensorPredictors/src/ttm.c
new file mode 100644
index 0000000..760dc43
--- /dev/null
+++ b/tensorPredictors/src/ttm.c
@@ -0,0 +1,115 @@
+// The need for `USE_FC_LEN_T` and `FCONE` is due to a Fortran character string
+// to C incompatibility. See: Writing R Extentions: 6.6.1 Fortran character strings
+#define USE_FC_LEN_T
+#include <R.h>
+#include <Rinternals.h>
+#include <R_ext/BLAS.h>
+#ifndef FCONE
+#define FCONE
+#endif
+
+/**
+ * Tensor Times Matrix a.k.a. Mode Product
+ *
+ * @param A multi-dimensionl array
+ * @param B matrix
+ * @param m mode index (1-indexed)
+ */
+extern SEXP ttm(SEXP A, SEXP B, SEXP m) {
+
+    // get zero indexed mode
+    int mode = asInteger(m) - 1;
+
+    // get dimension attribute of A
+    SEXP dim = getAttrib(A, R_DimSymbol);
+
+    // validate mode (mode must be smaller than the nr of dimensions)
+    if (mode < 0 || length(dim) <= mode) {
+        error("Illegal mode");
+    }
+
+    // and check if B is a matrix of non degenetate size
+    if (!isMatrix(B)) {
+        error("Expected a matrix as second argument");
+    }
+    if (!nrows(B) || !ncols(B)) {
+        error("Zero dimension detected");
+    }
+
+    // check matching of dimensions
+    if (INTEGER(dim)[mode] != ncols(B)) {
+        error("Dimension missmatch (mode dim not equal to ncol)");
+    }
+
+    // calc nr of response elements `prod(dim(X)[-mode]) * ncol(X)`,
+    int leny = 1;
+    // and the strides
+    //  `stride[0] <- prod(dim(X)[seq_len(mode - 1)])`
+    //  `stride[1] <- dim(X)[mode]`
+    //  `stride[2] <- prod(dim(X)[-seq_len(mode)])`
+    int stride[3] = {1, INTEGER(dim)[mode], 1};
+    for (int i = 0; i < length(dim); ++i) {
+        int size = INTEGER(dim)[i];
+        // check for non-degenetate dimensions
+        if (!size) {
+            error("Zero dimension detected");
+        }
+        leny *= (i == mode) ? nrows(B) : size;
+        stride[0] *= (i < mode) ? size : 1;
+        stride[2] *= (i > mode) ? size : 1;
+    }
+
+    // as well as the matrix rows (cols only needed once)
+    int nrow = nrows(B);
+
+    // create response object C
+    SEXP C = PROTECT(allocVector(REALSXP, leny));
+
+    // raw data access pointers
+    double* a = REAL(A);
+    double* b = REAL(B);
+    double* c = REAL(C);
+
+    // Tensor Times Matrix / Mode Product
+    const double zero = 0.0;
+    const double one  = 1.0;
+    if (mode == 0) {
+        // mode 1: (A x_1 B)_(1) = B A_(1) as a single Matrix-Matrix prod
+        F77_CALL(dgemm)("N", "N", &nrow, &stride[2], &stride[1], &one,
+            b, &nrow, a, &stride[1],
+            &zero, c, &nrow FCONE FCONE);
+    } else {
+        // Other modes can be written as blocks of matrix multiplications
+        for (int i2 = 0; i2 < stride[2]; ++i2) {
+            F77_CALL(dgemm)("N", "T", &stride[0], &nrow, &stride[1], &one,
+                &a[i2 * stride[0] * stride[1]], &stride[0], b, &nrow,
+                &zero, &c[i2 * stride[0] * nrow], &stride[0] FCONE FCONE);
+        }
+    }
+    /*
+    // Tensor Times Matrix / Mode Product (reference implementation)
+    memset(c, 0, leny * sizeof(double));
+    for (int i2 = 0; i2 < stride[2]; ++i2) {
+        for (int i1 = 0; i1 < stride[1]; ++i1) { // stride[1] == ncols(B)
+            for (int j = 0; j < nrow; ++j) {
+                for (int i0 = 0; i0 < stride[0]; ++i0) {
+                    c[i0 + (j + i2 * nrow) * stride[0]] +=
+                        a[i0 + (i1 + i2 * stride[1]) * stride[0]] * b[j + i1 * nrow];
+                }
+            }
+        }
+    }
+    */
+
+    // finally, set result dimensions
+    SEXP newdim = PROTECT(allocVector(INTSXP, length(dim)));
+    for (int i = 0; i < length(dim); ++i) {
+        INTEGER(newdim)[i] = (i == mode) ? nrows(B) : INTEGER(dim)[i];
+    }
+    setAttrib(C, R_DimSymbol, newdim);
+
+    // release C to the hands of the garbage collector
+    UNPROTECT(2);
+
+    return C;
+}