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

[Simon Urbanek]

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