Effective 3d array transposition with R/C++ (Rcpp) Task is reshape 3d array from [row, col, slice] to [slice,row,col]. I tried implement base::aperm analog with Rcpp.

// [[Rcpp::export]]
Rcpp::NumericVector array_perm(const Rcpp::NumericVector&amp; x) {
    if (Rf_isNull(x.attr("dim"))) {
        throw std::runtime_error("'x' does not have 'dim' attibute.");
    }
    Rcpp::Dimension d = x.attr("dim");
    if (d.size() != 3) {
        throw std::runtime_error("'x' must have 3 dimensions.");
    }
    std::size_t n = d[2];
    std::size_t n_vec = d[0] * d[1];
    std::size_t n_out = n_vec * n;
    Rcpp::NumericVector res = Rcpp::no_init(n_out);
    std::size_t ind_from = 0;
    for (std::size_t i = 0; i &lt; n; ++i) {
        std::size_t ind_to = i;
        for (std::size_t j = 0; j &lt; n_vec; ++j) {
            res[ind_to] = x[ind_from];
            ind_to += n;
            ind_from += 1;
        }
    }
    res.attr("dim") = Rcpp::Dimension(d[2], d[0], d[1]);
    return res;
}


How can I improve it?

For testing code:

# Observations
n &lt;- 1000
# Dimension
d &lt;- 100
# Array
a &lt;- replicate(n, matrix(seq_len(d * d), nrow = d, ncol = d))
# Desired result
res &lt;- aperm(a, c(3, 1, 2))
res


Small benchmark current variant of the code:

b &lt;- bench::mark(aperm(a, c(3, 1, 2)), array_perm(a))
b[, c("expression", "min", "median", "max", "itr/sec")]
#&gt;   expression                min   median      max `itr/sec`
#&gt;   &lt;chr&gt;                &lt;bch:tm&gt; &lt;bch:tm&gt; &lt;bch:tm&gt;     &lt;dbl&gt;
#&gt; 1 aperm(a, c(3, 1, 2))   84.9ms   85.1ms   85.5ms     11.7 
#&gt; 2 array_perm(a)         124.8ms  125.2ms  127.2ms      7.96


RcppArmadillo solution

// [[Rcpp::export]]
arma::Cube&lt;double&gt; array_perm2(const arma::Cube&lt;double&gt;&amp; x) {
    std::size_t rows = x.n_rows;
    std::size_t cols = x.n_cols;
    std::size_t slices = x.n_slices;
    arma::Cube&lt;double&gt; res(slices, rows, cols);
    for(std::size_t i = 0; i &lt; rows; ++i) {
        for(std::size_t j = 0; j &lt; cols; ++j) {
            for(std::size_t k = 0; k &lt; slices; ++k) {
                res(k, i, j) = x(i, j, k);
            }
        }
    }
    return res;
}


Benchmark:

  expression                min   median      max `itr/sec`
  &lt;chr&gt;                &lt;bch:tm&gt; &lt;bch:tm&gt; &lt;bch:tm&gt;     &lt;dbl&gt;
1 aperm(a, c(3, 1, 2))   85.8ms   86.4ms   87.7ms     11.6 
2 array_perm(a)         116.1ms  127.3ms  129.6ms      8.08
3 array_perm2(a)        217.4ms  219.7ms  222.1ms      4.55
If you can change the signature, to directly get cols, rows and slices count, you will not have to check for their validity.

Rcpp::NumericVector array_perm(const Rcpp::NumericVector&amp; input, 
                               std::size_t rows, std::size_t cols, std::size_t slices);


Otherwise, if you can change the error management, I think you'll get a speed boost. Exceptions handling come with a cost. Maybe return an empty vector? I don't know R , so I don't know possibilities.

You can also try to flattening your loops, here you have multiples options:

With the computations into output indexing

Rcpp::NumericVector array_perm(const Rcpp::NumericVector&amp; input, std::size_t rows, std::size_t cols, std::size_t slices) {
    // Think about the error management here...
    auto output = Rcpp::NumericVector(Rcpp::no_init(10));

    auto rc = rows * cols;
    auto size = rc * slices;
    for (std::size_t i = 0; i &lt; size; ++i) {
        output[i / rc + i % rc * slices] = input[i];
    }
    return output;
}


With the computations into intput looking

//...
    for (std::size_t i = 0; i &lt; size; ++i) output[i] = input[i/slices + i % slices * rc];
//...


Or in reverse order

//...
    while (size--) output[size] = input[size/slices + size % slices * rc];
//...


Or a mix

//...
    while (size--)  output[size / rc + size % rc * slices] = input[size];
//...


Or even a range-based for loop

//...
    std::size_t i = 0;
    for (auto e : input) {
        output[i / rc + i % rc * slices] = e;
        ++i;
    }
//...


PS: Did you tried with another contiguous storage type? (std::vector, plain old C array, ...)

PPS: I don't have R environment, so I only tested transposition algorithms with c++

Edit: 

Two other way:

//...
    for(std::size_t i = 0, j = 0; i &lt; size; ++i, j+=rc) {
        if (size &lt;= j) j -= size - 1;
        output[i] = input[j];
    }
//...
// or
//...
    for(std::size_t i = size, j = i-1; i--; j -= rc) {
        if (size &lt; j)  j += size-1;
        output[i] = input[j];
    }
//...