Merge pull request #20 from kevinushey/feature/boost-simd

add in boost.SIMD

Author: jjallaire

.Rbuildignore

@@ -8,3 +8,4 @@
 ^README.md
 ^doc$
 ^src/.*\.o$
+^examples/

.gitignore

@@ -4,3 +4,5 @@
 .DS_Store
 src-i386/
 src-x64/
+.Rprofile
+inst/doc

DESCRIPTION

@@ -18,13 +18,17 @@ Description: High level functions for doing parallel programming with 'Rcpp'.
   For example, the parallelFor() function can be used to convert the work of 
   a standard serial "for" loop into a parallel one and the parallelReduce() 
   function can be used for accumulating aggregate or other values.
-Suggests: Rcpp, RUnit
-LinkingTo: BH
+Suggests:
+    Rcpp,
+    RUnit,
+    knitr,
+    rmarkdown
+LinkingTo: BH (>= 1.60.0-1)
 SystemRequirements: GNU make, Windows: cmd.exe and cscript.exe, Solaris: g++ is required
 License: GPL-2
 URL: http://rcppcore.github.io/RcppParallel, https://github.com/RcppCore/RcppParallel
 Collate: 
    'build.R'
    'hooks.R'
    'options.R'
-
+VignetteBuilder: knitr

R/build.R

@@ -17,16 +17,25 @@ RcppParallelLibs <- function() {
 inlineCxxPlugin <- function() {
    list(
       env = list(
-         PKG_CXXFLAGS = tbbCxxFlags(),
+         PKG_CXXFLAGS = paste(tbbCxxFlags(), mtuneFlags()),
          PKG_LIBS = tbbLdFlags()
       ),
       includes = "#include <RcppParallel.h>",
-      LinkingTo = "RcppParallel",
-      body = function( x ) x,
+      LinkingTo = c("RcppParallel", "BH"),
+      body = function(x) x,
       Depends = "RcppParallel"
    )
 }
 
+mtuneFlags <- function() {
+   switch(Sys.info()[["sysname"]],
+          "Linux"   = "-mtune=native",
+          "Darwin"  = "-mtune=core2",
+          "Windows" = "-mtune=core2",
+          ""
+   )
+}
+
 tbbCxxFlags <- function() {
    flags <- "$(CXX1XSTD)"
    if (Sys.info()['sysname'] == "Windows")

TODO

@@ -0,0 +1,10 @@
+Boost.SIMD
+==========
+
+- Make it easy for client packages to use Boost.SIMD (plugin? Also need LinkingTo: BH?)
+- Enable 'parallelFor', 'parallelReduce'-style calls to SIMD
+- See 'Wizard' (http://www.evanmiller.org/) for inspiration
+- Need to think about NA (idiom that handled NA)
+- Setting of compiler flags for CRAN
+- Don't overlap with Armadillo, Eigen
+

examples/boost-simd-abssum.cpp

@@ -0,0 +1,49 @@
+// [[Rcpp::depends(RcppParallel)]]
+#include <RcppParallelSIMD.h>
+#include <Rcpp.h>
+using namespace Rcpp;
+
+#include <boost/simd/sdk/simd/pack.hpp>
+#include <boost/simd/memory/allocator.hpp>
+#include <boost/simd/include/functions/sum.hpp>
+#include <boost/simd/include/functions/load.hpp>
+#include <boost/simd/include/functions/plus.hpp>
+#include <boost/simd/include/functions/multiplies.hpp>
+#include <boost/simd/include/functions/abs.hpp>
+
+// [[Rcpp::export]]
+double simd_abssum(NumericVector x)
+{
+   using boost::simd::pack;
+   using boost::simd::load;
+   using boost::simd::aligned_load;
+   
+   typedef std::vector< double, boost::simd::allocator<double> > vector_type;
+   typedef pack<double> packed_type;
+   
+   vector_type data(x.begin(), x.end());
+   
+   packed_type packed;
+   double* it = &data[0];
+   double* end = &data[0] + data.size();
+   
+   while (it != end)
+   {
+      packed_type loaded = load<packed_type>(it);
+      packed = packed + boost::simd::abs(loaded);
+      it += packed_type::static_size;
+   }
+   
+   return sum(packed);
+}
+
+/*** R
+n <- 1024 * 1000
+x <- rnorm(n)
+
+library(microbenchmark)
+microbenchmark(
+   R = sum(abs(x)),
+   simd = simd_abssum(x)
+)
+*/

examples/boost-simd-accumulate.cpp

@@ -0,0 +1,38 @@
+// [[Rcpp::depends(RcppParallel)]]
+#include <RcppParallelSIMD.h>
+#include <Rcpp.h>
+using namespace Rcpp;
+
+#include <boost/simd/sdk/simd/algorithm.hpp>
+
+struct plus
+{
+   template <class T>
+   T operator()(const T& lhs, const T& rhs) const
+   {
+      return lhs + rhs;
+   }
+};
+
+// [[Rcpp::export]]
+double simd_sum(NumericVector x)
+{
+   return boost::simd::accumulate(x.begin(), x.end(), 0.0, plus());
+}
+
+// [[Rcpp::export]]
+double cpp_sum(NumericVector x)
+{
+   return std::accumulate(x.begin(), x.end(), 0.0, plus());
+}
+
+/***R
+n <- 1024 * 1000
+data <- rnorm(n)
+simd_sum(data)
+microbenchmark::microbenchmark(
+   R    = sum(data),
+   simd = simd_sum(data),
+   cpp  = cpp_sum(data)
+)
+*/

examples/boost-simd-capabilities.cpp

@@ -0,0 +1,37 @@
+// [[Rcpp::depends(RcppParallel)]]
+#include <RcppParallelSIMD.h>
+#include <Rcpp.h>
+using namespace Rcpp;
+
+#include <boost/simd/sdk/simd/extensions/meta/tags.hpp>
+#include <boost/simd/sdk/config/is_supported.hpp>
+#include <boost/config.hpp>
+
+// Try playing around with the '-march' flag to see what support is detected.
+// Note that most modern processors will implement instructions up to SSE3,
+// but certain AVX and FMA instructions are only available on newer CPUs.
+
+// [[Rcpp::export]]
+void simd_capabilities() {
+   
+   using boost::simd::is_supported;
+   using namespace boost::simd::tag;
+   
+   std::cout << "SIMD Capabilities\n"
+             << "=================\n\n";
+   
+   std::cout << "AVX2:   " << is_supported<avx2_>()   << "\n"
+             << "AVX:    " << is_supported<avx_>()    << "\n"
+             << "FMA4:   " << is_supported<fma4_>()   << "\n"
+             << "FMA3:   " << is_supported<fma3_>()   << "\n"
+             << "SSE4a:  " << is_supported<sse4a_>()  << "\n"
+             << "SSE4_2: " << is_supported<sse4_2_>() << "\n"
+             << "SSE4_1: " << is_supported<sse4_1_>() << "\n"
+             << "SSE3:   " << is_supported<sse3_>()   << "\n"
+             << "SSE2:   " << is_supported<sse2_>()   << "\n"
+             << "SSE:    " << is_supported<sse_>()    << "\n";
+}
+
+/***R
+simd_capabilities()
+*/

examples/boost-simd-dot.cpp

@@ -0,0 +1,69 @@
+// [[Rcpp::depends(RcppParallel)]]
+#include <RcppParallelSIMD.h>
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// http://nt2.metascale.fr/doc/html/tutorials/processing_data_the_simd_way.html
+#include <boost/simd/memory/allocator.hpp>
+#include <boost/simd/sdk/simd/pack.hpp>
+#include <boost/simd/include/functions/sum.hpp>
+#include <boost/simd/include/functions/load.hpp>
+#include <boost/simd/include/functions/plus.hpp>
+#include <boost/simd/include/functions/multiplies.hpp>
+
+template <typename Value>
+Value simd_dot(Value* first1, Value* last1, Value* first2)
+{
+   using boost::simd::sum;
+   using boost::simd::pack;
+   
+   typedef pack<Value> type;
+   type tmp;
+   
+   // Let's consider that (last1-first1) is divisible by the size of the pack.
+   while (first1 != last1)
+   {
+      // Load current values from the datasets
+      pack<Value> x1 = boost::simd::aligned_load<Value>(first1);
+      pack<Value> x2 = boost::simd::aligned_load<Value>(first2);
+      
+      // Computation
+      tmp = tmp + x1 * x2;
+      
+      // Advance to the next SIMD vector
+      first1 += type::static_size;
+      first2 += type::static_size;
+   }
+   
+   return sum(tmp);
+}
+
+// [[Rcpp::export]]
+double dot(NumericVector lhs, NumericVector rhs)
+{
+   // construct simd vectors
+   typedef std::vector< double, boost::simd::allocator<double> > vector_t;
+   vector_t a(lhs.begin(), lhs.end());
+   vector_t b(rhs.begin(), rhs.end());
+   
+   // call dot function
+   double result = simd_dot(&a[0], &a[0] + a.size(), &b[0]);
+   
+   return result;
+}
+
+/*** R
+n <- 1024
+lhs <- rnorm(n)
+rhs <- rnorm(n)
+result <- dot(lhs, rhs)
+all.equal(result, sum(lhs * rhs))
+
+library(microbenchmark)
+lhs <- rnorm(n * 1000)
+rhs <- rnorm(n * 1000)
+microbenchmark(
+   simd = dot(lhs, rhs),
+   R = sum(lhs * rhs)
+)
+*/

examples/boost-simd-hello-world.cpp

@@ -0,0 +1,30 @@
+// [[Rcpp::depends(RcppParallel)]]
+#include <RcppParallelSIMD.h>
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// http://nt2.metascale.fr/doc/html/tutorials/simd_hello_world.html
+#include <boost/simd/sdk/simd/pack.hpp>
+#include <boost/simd/sdk/simd/io.hpp>
+#include <boost/simd/include/functions/splat.hpp>
+#include <boost/simd/include/functions/plus.hpp>
+#include <boost/simd/include/functions/multiplies.hpp>
+#include <iostream>
+
+// [[Rcpp::export]]
+void HelloWorld()
+{
+   typedef boost::simd::pack<float> p_t;
+   
+   p_t res;
+   p_t u(10);
+   p_t r = boost::simd::splat<p_t>(11);
+   
+   res = (u + r) * 2.f;
+   
+   Rcout << res << std::endl;
+}
+
+/*** R
+HelloWorld()
+*/