Re: [Rd] How to safely using OpenMP pragma inside a .C() function?

From: Simon Urbanek <simon.urbanek_at_r-project.org>
Date: Wed, 31 Aug 2011 08:19:13 -0400

On Aug 30, 2011, at 12:57 PM, pawelm wrote:

> Simon,
>
> I found that files R-2.13.1/src/library/stats/src/distance.c and
> R-2.13.1/src/main/array.c have openmp code (example below). I have couple
> questions regarding best practices when using R internals and openmp.
>
> Can we use R-2.13.1/src/library/stats/src/distance.c and
> R-2.13.1/src/main/array.c as an example how to interact with R code and R
> internals?

Technically, close, but not quite, because R internals actually use different code than R packages: in R internals calls to REAL, LENGTH etc. are macros and thus operate directly on the object (and thus optimizable), whereas in packages those are function calls (thus they do change the stack and are not optimizable!). This implies that in theory you can't use any R API calls (and things like REAL are function calls) since you have absolutely no idea what they may touch (e.g., they could - in theory - move things around since everything is serial in R which would bomb any OMP part). In practice they are more harmless, but my point is that you may want to be on the safe side, especially in the simple cases where you can avoid them.

For a trivial example you can use something like:

#include <math.h>
#include <Rinternals.h>

SEXP omp_dnorm(SEXP s_x, SEXP s_mu, SEXP s_sigma) {

    double *x = REAL(s_x);
    R_len_t n = LENGTH(s_x);
    const double mu = asReal(s_mu), sigma = asReal(s_sigma), div = sigma * sqrt(2 * M_PI);     SEXP res = allocVector(REALSXP, n);
    double *y = REAL(res);
    R_len_t i;

#pragma omp parallel for firstprivate(x, y, n) default(none)

    for (i = 0; i < n; i++)

        y[i] = exp(-0.5 * ((x[i] - mu)/sigma) * ((x[i] - mu)/sigma)) / div;

    return res;
}

> dyn.load("ompx.so")
> x = 1:1e7/1e6
> system.time(.Call("omp_dnorm", x, 0, 1))

   user system elapsed
  1.680 0.020 0.095
> system.time(dnorm(x, 0, 1))

   user system elapsed
  2.410 0.020 0.658

Cheers,
Simon

> What are my other options if I want to work with SEXP structures in my
> parallel code?
>
> Thank you
> Regards
>
> =============
>
> #ifdef HAVE_OPENMP
> /* This gives a spurious -Wunused-but-set-variable error */
> if (R_num_math_threads > 0)
> nthreads = R_num_math_threads;
> else
> nthreads = 1; /* for now */
> #pragma omp parallel for num_threads(nthreads) default(none) \
> private(j, i, ix, rx) \
> firstprivate(x, ans, n, p, type, cnt, sum, \
> NaRm, keepNA, R_NaReal, R_NaInt, OP)
> #endif
> for (j = 0; j < p; j++) {
> switch (type) {
> case REALSXP:
> rx = REAL(x) + n*j;
> if (keepNA)
> for (sum = 0., i = 0; i < n; i++) sum += *rx++;
> else {
> for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
> if (!ISNAN(*rx)) {cnt++; sum += *rx;}
> else if (keepNA) {sum = NA_REAL; break;}
> }
> break;
> case INTSXP:
> ix = INTEGER(x) + n*j;
> for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
> if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
> else if (keepNA) {sum = NA_REAL; break;}
> break;
> case LGLSXP:
> ix = LOGICAL(x) + n*j;
> for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
> if (*ix != NA_LOGICAL) {cnt++; sum += *ix;}
> else if (keepNA) {sum = NA_REAL; break;}
> break;
> default:
> /* we checked the type above, but be sure */
> UNIMPLEMENTED_TYPEt("do_colsum", type);
> }
> if (OP == 1) {
> if (cnt > 0) sum /= cnt; else sum = NA_REAL;
> }
> REAL(ans)[j] = sum;
> }
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/How-to-safely-use-OpenMP-pragma-inside-a-C-function-tp3777036p3779214.html
> Sent from the R devel mailing list archive at Nabble.com.
>
> ______________________________________________
> R-devel_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed 31 Aug 2011 - 12:22:48 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

All messages

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 31 Aug 2011 - 14:20:25 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive