Hello Rcpp

This past weekend I discovered the wonders of c++ thanks to this datacamp course. Although c++ syntax is different, knowing Fortran made this much easier.

Filling a matrix with c++

The following code creates a function that can be called from R to fill a matrix. Something that is different than in Fortran is that to make loops more efficient you have to do right (j) to left (i) instead of left to right. The other big difference, c++ starts counting at 0 instead of 1! This one got me more than once…

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix my_matrix(int I, int J) {
  NumericMatrix A(I,J);
  int iter = 0;
  for(int j = 0; j < J; j++) {
    for(int i = 0; i < I; i++){
      iter++;
      A(i,j) = iter ;
    }
  }
  return A;
}

my_matrix(6,5)

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    7   13   19   25
## [2,]    2    8   14   20   26
## [3,]    3    9   15   21   27
## [4,]    4   10   16   22   28
## [5,]    5   11   17   23   29
## [6,]    6   12   18   24   30

Similarly to Fortran, you can use OpenMP do speed things up very easily:

#include <Rcpp.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>

// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix my_matrix(int I, int J, int nthreads) {
  NumericMatrix A(I,J);
  // create a thread safe accessor for A
  RcppParallel::RMatrix<double> a(A);
  int tid;
  omp_set_num_threads(nthreads);
#pragma omp parallel for private(tid)
  for(int j = 0; j < J; j++) {
    for(int i = 0; i < I; i++) {
      tid = omp_get_thread_num();
      a(i, j) = tid ;
    }
  }

  return A;
}

my_matrix(6,5,4)

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    1    2    3
## [2,]    0    0    1    2    3
## [3,]    0    0    1    2    3
## [4,]    0    0    1    2    3
## [5,]    0    0    1    2    3
## [6,]    0    0    1    2    3

In this case a is a safe accessor for A. For more details check RcppParallel.

My conclusion: although Fortran is faster than c++, writing c++ code using Rcpp is way easier than writing Fortran code, and my time is way more valuable than the small difference that I may gain by using Fortran.