[Rcpp-devel] Rcpp Parallel and Rcpp Armadillo

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/
<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)