Re: [Rcpp-devel] Parallel R::bessel_k function in Rcpp
Hi, Check also the Bessel package from Martin M., https://cran.r-project.org/web/packages/Bessel/index.html It wraps some C code; maybe a C++ interface wouldn't be too hard to write. There's also the GSL package apparently. Cheers, baptiste On 24 March 2016 at 13:20, Hoang Nguyen wrote: > Thank Dirk, > I check Boost and they only have the Modified Bessel Functions of the > First and Second Kinds. I have to read more to find how to derive the third > kind. > Best regards. > > On Thu, Mar 24, 2016 at 1:17 AM, Dirk Eddelbuettel wrote: > >> >> PS Looks like Boost has Bessel functions: >> >> >> http://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/math_toolkit/bessel/bessel_over.html >> >> So maybe try that. We should have access to that in the BH package. >> >> Dirk >> >> -- >> http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org >> > > > ___ > Rcpp-devel mailing list > Rcpp-devel@lists.r-forge.r-project.org > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel > ___ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Re: [Rcpp-devel] Parallel R::bessel_k function in Rcpp
Hi, I found in wiki that Modified Bessel Functions of the Second Kinds has different names and one of it is the third kind. Sorry for this trouble. Thank you. Best. On Thu, Mar 24, 2016 at 1:20 AM, Hoang Nguyen wrote: > Thank Dirk, > I check Boost and they only have the Modified Bessel Functions of the > First and Second Kinds. I have to read more to find how to derive the third > kind. > Best regards. > > On Thu, Mar 24, 2016 at 1:17 AM, Dirk Eddelbuettel wrote: > >> >> PS Looks like Boost has Bessel functions: >> >> >> http://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/math_toolkit/bessel/bessel_over.html >> >> So maybe try that. We should have access to that in the BH package. >> >> Dirk >> >> -- >> http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org >> > > ___ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Re: [Rcpp-devel] Parallel R::bessel_k function in Rcpp
Thank Dirk, I check Boost and they only have the Modified Bessel Functions of the First and Second Kinds. I have to read more to find how to derive the third kind. Best regards. On Thu, Mar 24, 2016 at 1:17 AM, Dirk Eddelbuettel wrote: > > PS Looks like Boost has Bessel functions: > > > http://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/math_toolkit/bessel/bessel_over.html > > So maybe try that. We should have access to that in the BH package. > > Dirk > > -- > http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org > ___ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Re: [Rcpp-devel] Parallel R::bessel_k function in Rcpp
PS Looks like Boost has Bessel functions: http://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/math_toolkit/bessel/bessel_over.html So maybe try that. We should have access to that in the BH package. Dirk -- http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org ___ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Re: [Rcpp-devel] Parallel R::bessel_k function in Rcpp
On 24 March 2016 at 01:02, Hoang Nguyen wrote:
| Hi,
| I am working with Rcpp and need to call function bessel_k from R::bessel_k in
| parallel. I got usually crash in Rstudio, and still not find out which is the
| mistake
You _cannot call back to R_ from parallel code.
| Here is my code in R:
|
| Sys.setenv("PKG_CXXFLAGS"="-fopenmp")
| Sys.setenv("PKG_LIBS"="-fopenmp")
| library(Rcpp)
| sourceCpp("Bessel.cpp")
| xseq = rt(1,df = 8)
| param <- c(0.000, 2.8284271, -0.4954229, 8.000)
| BasselCm <- BasselC(xseq,param)
|
| here is my code in Bessel.cpp:
| if I set num threads(2) it was crash but run smooth if I set back to 1.
Sorry but I fear you either
-- must write a new C++ implementation of the Bessel_k, or
-- if you want to use R's, you need to use process-parallelism (ie
something like OpenMPI to separate R processes).
Dirk
|
| Thank you so much for your consideration.
|
| Best regards.
|
|
|
|
| > sessionInfo()
| R version 3.2.4 Revised (2016-03-16 r70336)
| Platform: x86_64-pc-linux-gnu (64-bit)
| Running under: Ubuntu 14.04.4 LTS
|
| locale:
| [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=
| es_ES.UTF-8
| [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=
| en_US.UTF-8
| [7] LC_PAPER=es_ES.UTF-8 LC_NAME=C LC_ADDRESS=
| C
| [10] LC_TELEPHONE=C LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=
| C
|
| attached base packages:
| [1] stats graphics grDevices utils datasets methods base
|
| other attached packages:
| [1] Rcpp_0.12.3
|
| loaded via a namespace (and not attached):
| [1] tools_3.2.4 RcppArmadillo_0.6.600.4.0
|
|
| --
|
|
| [DELETED ATTACHMENT Bessel.cpp, text/x-c++src]
| ___
| Rcpp-devel mailing list
| Rcpp-devel@lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
--
http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org
___
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
[Rcpp-devel] Parallel R::bessel_k function in Rcpp
Hi,
I am working with Rcpp and need to call function bessel_k from *R::bessel_k*
in parallel. I got usually crash in Rstudio, and still not find out which
is the mistake
Here is my code in R:
Sys.setenv("PKG_CXXFLAGS"="-fopenmp")
Sys.setenv("PKG_LIBS"="-fopenmp")
library(Rcpp)
sourceCpp("Bessel.cpp")
xseq = rt(1,df = 8)
param <- c(0.000, 2.8284271, -0.4954229, 8.000)
BasselCm <- BasselC(xseq,param)
here is my code in Bessel.cpp:
if I set num threads(2) it was crash but run smooth if I set back to 1.
Thank you so much for your consideration.
Best regards.
> sessionInfo()
R version 3.2.4 Revised (2016-03-16 r70336)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
LC_TIME=es_ES.UTF-8
[4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=es_ES.UTF-8
LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=es_ES.UTF-8 LC_NAME=C
LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=es_ES.UTF-8
LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] Rcpp_0.12.3
loaded via a namespace (and not attached):
[1] tools_3.2.4 RcppArmadillo_0.6.600.4.0
--
#include
#include
#include
// [[Rcpp::plugins(openmp)]]
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma ;
using namespace Rcpp ;
// [[Rcpp::export]]
List BasselC(SEXP dataxseq_, SEXP param_) {
BEGIN_RCPP
int seed = 126;
std::srand(seed);
// Init hyperparams
Rcpp::NumericVector param(param_);
double mu = param[0];
double delta = param[1];
double beta = param[2];
double nu = param[3];
// Rcout << "mu " << mu << " delta " << delta << " beta " << beta << " nu " << nu << std::endl;
arma::vec xseq = as(dataxseq_);
arma::vec Besselxseq = xseq;
int maxrows = xseq.n_rows;
#ifdef _OPENMP
omp_set_num_threads(2);
#endif
#pragma omp parallel for default(none) shared(xseq,mu,delta,beta,nu,maxrows,Besselxseq)
for (int i =0; i < maxrows; i++){
for (int j =0; j < 5000; j++){ //Only for making more computation
Besselxseq(i) = log(R::bessel_k(sqrt(pow(beta,2)*(pow(delta,2) + pow(xseq(i) - mu,2))),
(nu + 1)/2, 2)) ;
}
}
List z= List::create(Rcpp::Named("Besselxseq") = Besselxseq);
return z;
PutRNGstate();
END_RCPP
}___
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
