[Rcpp-devel] Rcpp Parallel and Rcpp Armadillo

[Maxime To]

Hi,

I started to work with Rcpp and I have some trouble to adapt some code. I copy below the following code:

- I want to write a function named contrib1, that is simplified version here, but the idea is that that function makes a large use of the pnorm (and qnorm) on marices and I would like to make it faster by applying the function using Rcpp parallel.

- Then I started with an example that I tried to adapt (http://gallery.rcpp.org/articles/parallel-matrix-transform/) to build the structure Pnorm and the function pPnorm. I also created the function f (I wasn't able to use directly the pnorm function as would Sugar allow...)

The example doesn't work. And I don't really understand why. I guess that it is a problem of memory allocation... Can anyone help me on that point?? Moreover, I was able to build a previous version that was working, but once I started to iterate on the contrib1 function (I execute an optimization algorithm), the RAM memory saturated I wasn't able to figure out why the program was not freeing memory after each execution of the function.

I attach the code below. I am a Rcpp novice, so don't hesitate to comment on strange things....
Thanks for your answers!

Maxime

--------------------------------------

#include <RcppArmadillo.h>
#include <cmath>
#include <algorithm>
#include <RcppParallel.h>

using namespace Rcpp;
using namespace arma;
using namespace std;
using namespace RcppParallel;

// [[Rcpp::depends(RcppArmadillo, RcppParallel)]]

inline double f(double x) { return ::Rf_pnorm5(x, 0.0, 1.0, 1, 0); }

struct Pnorm : public Worker
{
   // source matrix
   mat input;
   
   // destination matrix
   mat output;
   
   // initialize with source and destination
   Pnorm(mat input, mat output) 
      : input(input), output(output) {}
   
   // take the square root of the range of elements requested
   void operator()(std::size_t begin, std::size_t end) {
      std::transform(input.begin() + begin, 
                     input.begin() + end, 
                     output.begin() + begin, 
                     f);
   }
};

mat pPnorm(mat x) {
  
    // allocate the output matrix
    mat output(x.n_rows, x.n_cols);
  
    // Pnorm functor (pass input and output matrixes)
    Pnorm ppnorm(x, output);
  
    // call parallelFor to do the work
    parallelFor(0, x.n_elem, ppnorm);
  
    // return the output matrix
    mat outmat(output.begin(), output.n_rows, output.n_cols);
      return outmat;
}

// [[Rcpp::export]]
mat contrib1(mat my_matrix) {

        mat pbound2    = pPnorm(my_matrix);
    return pbound2;
}


------------------------------------------
### Main R code

rm(list=ls(all=TRUE))

library(Rcpp)
library(RcppGSL)
library(RcppParallel)
library(RcppArmadillo)
setwd('U:/testParallel')

sourceCpp("test.cpp")

asd = matrix(runif(1000), ncol = 100)
contrib1(asd)


 		 	   		  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20141209/8fdd907f/attachment.html>