How to safely use OpenMP pragma inside a .C() function?

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
13 messages Options
Reply | Threaded
Open this post in threaded view
|

How to safely use OpenMP pragma inside a .C() function?

Alireza Mahani
This post was updated on .
I am trying to parallelize part of a C function that is called from R (via .C) using OpenMP's "parallel for" pragma. I get mixed results: some runs finish with no problem, but some lead to R crashing after issuing a long error message involving memory violations.

I found this post, which describes how a .Call() function can be made to avoid crashing R by raising the stack limit:

http://www.r-bloggers.com/using-openmp-ized-c-code-with-r/

However, trying this in my .C function doesn't help things. Any suggestions/tips on whether I can safely use OpenMP inside a .C function, and if yes how?

Thank you,
Alireza Mahani
Reply | Threaded
Open this post in threaded view
|

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

Simon Urbanek

On Aug 29, 2011, at 7:48 PM, Alireza Mahani wrote:

> I am trying to parallelize part of a C function that is called from R (via
> .C) using OpenMP's "parallel for" pragma. I get mixed results: some runs
> finish with no problem, but some lead to R crashing after issuing a long
> error message involving memory violations.
>

You'll have to provide the code. In general it works (even R uses it itself), but there are strict requirements (no R API calls) that you must adhere to.


> I found this post, which describes how a .Call() function can be made to
> avoid crashing R by raising the stack limit:
>
> http://www.r-bloggers.com/using-openmp-ized-c-code-with-r/
>

I skimmed through the post and all of the examples are broken - they will only work (incidentally) as R internals, not officially (and they are unnecessary inefficient).


> However, trying this in my .C function doesn't help things. Any
> suggestions/tips on whether I can safely use OpenMP inside a .C function,
> and if yes how?
>

There are issues with OpenMP on some platforms in general (in fact pretty much every platform had some issue at some point in time), but apart from those you only have to make sure that you declare shared/private variables properly and don't use *any* R API calls in the parallel part (this includes things like LENGTH, REAL, ...).

Cheers,
Simon



> Thank you,
> Alireza Mahani
>
> --
> View this message in context: http://r.789695.n4.nabble.com/How-to-safely-using-OpenMP-pragma-inside-a-C-function-tp3777036p3777036.html
> Sent from the R devel mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

Alireza Mahani
I have no R API calls inside the parallelized block. I will work on creating a self-contained example and post it for your review. Thanks! -Alireza
Reply | Threaded
Open this post in threaded view
|

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

pawelm
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?
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;
        }
Reply | Threaded
Open this post in threaded view
|

Re: 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.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

Simon Urbanek
In reply to this post by pawelm
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?
> What are my other options if I want to work with SEXP structures in my
> parallel code?
>

I forgot to answer this one in my previous response, sorry.

The short answer is none. If you can't access what you need before the omp loop you may be in trouble. For example you can't allocate anything in the parallel code, so you can't create character elements or assign objects. There are some things (very few) you can do like using ISNA() or ISNAN() but the trouble is that without looking at the sources you don't know what you can do (note that the pragma below declares R_NaReal although it actually uses NA_REAL -- in theory that is only possible in internal R code because it knows that NA_REAL is defined as R_NaReal variable and not a fixed constant - which is entirely a matter of implementation and thus not under your control).

Cheers,
Simon


> 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.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

pawelm
Simon,

This is very useful example and explanation. Thank you very, very much. The icing on the cake would be some guidelines how to set up the number of threads. R source code uses global variable R_num_math_threads. Can we use that? Or each openmp-enabled R package would have it's own mechanism?

Related question: I was trying to find out if R_num_math_threads could be set up by a user and I found this line in the ./src/main/names.c file:

{"setNumMathThreads", do_setnumthreads, 0, 11, 1, {PP_FUNCALL, PREC_FN, 0}}

It suppose to be a function name but all I get in R command line is:

Error: could not find function "setNumMathThreads"

Any thoughts?

Thank you
Regards
Reply | Threaded
Open this post in threaded view
|

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

pawelm
I just found this (performance improvement of the "dist" function when using openmp):

.Internal(setMaxNumMathThreads(1)); .Internal(setNumMathThreads(1)); m <- matrix(rnorm(810000),900,900); system.time(d <- dist(m))

  user  system elapsed
  3.510   0.013   3.524

.Internal(setMaxNumMathThreads(5)); .Internal(setNumMathThreads(5)); m <- matrix(rnorm(810000),900,900); system.time(d <- dist(m));

   user  system elapsed
  3.536   0.007   1.321

Works great! Just the question stays if it's a good practice to use "R_num_math_threads" in external packages?

Thanks
Reply | Threaded
Open this post in threaded view
|

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

Simon Urbanek
Pawel,

On Aug 31, 2011, at 4:46 PM, pawelm wrote:

> I just found this (performance improvement of the "dist" function when using
> openmp):
>
> .Internal(setMaxNumMathThreads(1)); .Internal(setNumMathThreads(1)); m <-
> matrix(rnorm(810000),900,900); system.time(d <- dist(m))
>
>  user  system elapsed
>  3.510   0.013   3.524
>
> .Internal(setMaxNumMathThreads(5)); .Internal(setNumMathThreads(5)); m <-
> matrix(rnorm(810000),900,900); system.time(d <- dist(m));
>
>   user  system elapsed
>  3.536   0.007   1.321
>
> Works great! Just the question stays if it's a good practice to use
> "R_num_math_threads" in external packages?
>

Normally you don't need to mess with all this and I would recommend not to do so. The R internals use a different strategy since they need to cope with the fall-back case, but packages should not worry about that. The default number of threads is defined by the OMP_NUM_THREADS environment variable and that is the documented way in OpenMP, so my recommendation would be to not mess with num_threads() which is precisely why I did not use it in the example I gave you.

That said, R-devel has new facilities for parallelization so things may change in the future.

Cheers,
Simon

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

Prof Brian Ripley
Note that currently R internals do not actually use multiple threads
in OpenMP, and there is no documented way to make them do so.

The main issue is that there is insufficient knowlege of where they
are worthwhile (which is both OS and platform-dependent: we don't even
have reliable cross-platform ways to decide a reasonable number of
threads, and the number of virtual cores on a multi-user platform
definitely is not reasonable).  Luke Tierney reported that the
crossover point for a speed-up on Mac OS X was much larger matrices
than on Linux, for example, and there is currently no OpenMP support
in the Windows toolchain.

The current implementation is a trial: there are more places planned
to use OpenMP as and when the uncertainties are resolved.

This will change at some point: given the current instability in
thread support in the MinGW-w64 project this may or may not be before
R 2.14.0.

On Wed, 31 Aug 2011, Simon Urbanek wrote:

> Pawel,
>
> On Aug 31, 2011, at 4:46 PM, pawelm wrote:
>
>> I just found this (performance improvement of the "dist" function when using
>> openmp):

You failed to describe the platform!  See the posting guide (which
asked you to do so 'at a minimum').

>> .Internal(setMaxNumMathThreads(1)); .Internal(setNumMathThreads(1)); m <-
>> matrix(rnorm(810000),900,900); system.time(d <- dist(m))
>>
>>  user  system elapsed
>>  3.510   0.013   3.524
>>
>> .Internal(setMaxNumMathThreads(5)); .Internal(setNumMathThreads(5)); m <-
>> matrix(rnorm(810000),900,900); system.time(d <- dist(m));
>>
>>   user  system elapsed
>>  3.536   0.007   1.321
>>
>> Works great! Just the question stays if it's a good practice to use
>> "R_num_math_threads" in external packages?

Most definitely not: it is never good practice to use undocumented
non-API variables. See 'Writing R Extensions'.

> Normally you don't need to mess with all this and I would recommend
> not to do so. The R internals use a different strategy since they
> need to cope with the fall-back case, but packages should not worry
> about that. The default number of threads is defined by the
> OMP_NUM_THREADS environment variable and that is the documented way
> in OpenMP, so my recommendation would be to not mess with
> num_threads() which is precisely why I did not use it in the example
> I gave you.

I'd be cautious there.  OMP_NUM_THREADS affects all the OpenMP code in
the R session, and possibly others which use it (some parallel BLAS do
too).

>
> That said, R-devel has new facilities for parallelization so things
> may change in the future.
>
> Cheers,
> Simon

--
Brian D. Ripley,                  [hidden email]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

pawelm
>>> .Internal(setMaxNumMathThreads(1)); .Internal(setNumMathThreads(1)); m <-
>>> matrix(rnorm(810000),900,900); system.time(d <- dist(m))
>>>
>>>  user  system elapsed
>>>  3.510   0.013   3.524
>>>
>>> .Internal(setMaxNumMathThreads(5)); .Internal(setNumMathThreads(5)); m <-
>>> matrix(rnorm(810000),900,900); system.time(d <- dist(m));
>>>
>>>  user  system elapsed
>>>  3.536   0.007   1.321

These results are obtained on:

  Model Name: Mac Pro
  Model Identifier: MacPro4,1
  Processor Name: Quad-Core Intel Xeon
  Processor Speed: 2.93 GHz
  Number Of Processors: 2
  Total Number Of Cores: 8
  L2 Cache (per core): 256 KB
  L3 Cache (per processor): 8 MB
  Memory: 12 GB
  Processor Interconnect Speed: 6.4 GT/s
  Boot ROM Version: MP41.0081.B07
  SMC Version (system): 1.39f5
  SMC Version (processor tray): 1.39f5

Thank you for your comments
Regards
--
Pawel Matykiewicz
http://www.neuron.m4u.pl
http://www.linkedin.com/in/pawelmatykiewicz

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

ghostwheel
In reply to this post by Prof Brian Ripley
This is probably obvious, but I just wanted to say that it should be possible to turn off multithreading even when on a machine with multiple cores.
Reasons could be because you run in a cluster, and are given just one core for yourself.
Or, if you have a setup with trivial parallelization (i.e. you run almost the same task 100 times), where getting a speedup of 100 fold on 100 cores is easy, and multithreading would slow you down.

Michael


On 1 Sep 2011, at 8:34AM, Prof Brian Ripley wrote:

> Note that currently R internals do not actually use multiple threads in OpenMP, and there is no documented way to make them do so.
>
> The main issue is that there is insufficient knowlege of where they are worthwhile (which is both OS and platform-dependent: we don't even have reliable cross-platform ways to decide a reasonable number of threads, and the number of virtual cores on a multi-user platform definitely is not reasonable).  Luke Tierney reported that the crossover point for a speed-up on Mac OS X was much larger matrices than on Linux, for example, and there is currently no OpenMP support in the Windows toolchain.
>
> The current implementation is a trial: there are more places planned to use OpenMP as and when the uncertainties are resolved.
>
> This will change at some point: given the current instability in thread support in the MinGW-w64 project this may or may not be before R 2.14.0.
>
> On Wed, 31 Aug 2011, Simon Urbanek wrote:
>
>> Pawel,
>>
>> On Aug 31, 2011, at 4:46 PM, pawelm wrote:
>>
>>> I just found this (performance improvement of the "dist" function when using
>>> openmp):
>
> You failed to describe the platform!  See the posting guide (which asked you to do so 'at a minimum').
>
>>> .Internal(setMaxNumMathThreads(1)); .Internal(setNumMathThreads(1)); m <-
>>> matrix(rnorm(810000),900,900); system.time(d <- dist(m))
>>>
>>> user  system elapsed
>>> 3.510   0.013   3.524
>>>
>>> .Internal(setMaxNumMathThreads(5)); .Internal(setNumMathThreads(5)); m <-
>>> matrix(rnorm(810000),900,900); system.time(d <- dist(m));
>>>
>>>  user  system elapsed
>>> 3.536   0.007   1.321
>>>
>>> Works great! Just the question stays if it's a good practice to use
>>> "R_num_math_threads" in external packages?
>
> Most definitely not: it is never good practice to use undocumented non-API variables. See 'Writing R Extensions'.
>
>> Normally you don't need to mess with all this and I would recommend not to do so. The R internals use a different strategy since they need to cope with the fall-back case, but packages should not worry about that. The default number of threads is defined by the OMP_NUM_THREADS environment variable and that is the documented way in OpenMP, so my recommendation would be to not mess with num_threads() which is precisely why I did not use it in the example I gave you.
>
> I'd be cautious there.  OMP_NUM_THREADS affects all the OpenMP code in the R session, and possibly others which use it (some parallel BLAS do too).
>
>>
>> That said, R-devel has new facilities for parallelization so things may change in the future.
>>
>> Cheers,
>> Simon
>
> --
> Brian D. Ripley,                  [hidden email]
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford,             Tel:  +44 1865 272861 (self)
> 1 South Parks Road,                     +44 1865 272866 (PA)
> Oxford OX1 3TG, UK                Fax:  +44 1865 272595
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

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

Martin Morgan
In reply to this post by Simon Urbanek
On 08/29/2011 05:36 PM, Simon Urbanek wrote:

>
> On Aug 29, 2011, at 7:48 PM, Alireza Mahani wrote:
>
>> I am trying to parallelize part of a C function that is called from
>> R (via .C) using OpenMP's "parallel for" pragma. I get mixed
>> results: some runs finish with no problem, but some lead to R
>> crashing after issuing a long error message involving memory
>> violations.
>>
>
> You'll have to provide the code. In general it works (even R uses it
> itself), but there are strict requirements (no R API calls) that you
> must adhere to.

Hi Simon et al.,

I'm trying to resolve this dictum (no R API calls) with my understanding
of OpenMP's shared memory model. For instance I would have thought R API
calls inside critical sections would be safe, as this (though obviously
gratuitous) appears to be on my own machine (which I know is no way to
arrive at a general understanding!)

#include <omp.h>
#include <Rdefines.h>

SEXP omp()
{
     const int n = 40;
     SEXP res = PROTECT(allocVector(INTSXP, n));
     const int it = 1000;
     int j;
     for (j = 0; j < Rf_length(res); ++j)
         INTEGER(res)[j] = 0;

     omp_set_num_threads(n);
     j = 0;
#pragma omp parallel for
     for (int i = 0; i < it; ++i)
#pragma omp critical
     {
         j  = omp_get_thread_num();
         SEXP elt = PROTECT(Rf_ScalarInteger(1));
         INTEGER(res)[j] = INTEGER(res)[j] + INTEGER(elt)[0];
         UNPROTECT(1);
     }

     UNPROTECT(1);
     return res;
}

and

 > dyn.load("omp.so"); table(x <- .Call("omp"))

25
40


?

Martin

 > sessionInfo()
R Under development (unstable) (2011-09-12 r56997)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
  [7] LC_PAPER=C                 LC_NAME=C
  [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  utils     datasets  grDevices methods   base

loaded via a namespace (and not attached):
[1] tools_2.14.0


>
>
>> I found this post, which describes how a .Call() function can be
>> made to avoid crashing R by raising the stack limit:
>>
>> http://www.r-bloggers.com/using-openmp-ized-c-code-with-r/
>>
>
> I skimmed through the post and all of the examples are broken - they
> will only work (incidentally) as R internals, not officially (and
> they are unnecessary inefficient).
>
>
>> However, trying this in my .C function doesn't help things. Any
>> suggestions/tips on whether I can safely use OpenMP inside a .C
>> function, and if yes how?
>>
>
> There are issues with OpenMP on some platforms in general (in fact
> pretty much every platform had some issue at some point in time), but
> apart from those you only have to make sure that you declare
> shared/private variables properly and don't use *any* R API calls in
> the parallel part (this includes things like LENGTH, REAL, ...).
>
> Cheers, Simon
>
>
>
>> Thank you, Alireza Mahani
>>
>> -- View this message in context:
>> http://r.789695.n4.nabble.com/How-to-safely-using-OpenMP-pragma-inside-a-C-function-tp3777036p3777036.html
>>
>>
Sent from the R devel mailing list archive at Nabble.com.
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>>
>
> ______________________________________________ [hidden email]
> mailing list https://stat.ethz.ch/mailman/listinfo/r-devel


--
Computational Biology
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109

Location: M1-B861
Telephone: 206 667-2793

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel