<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/r-gpu-programming-for-all-with-gpur/"></div></div><div><a data-url="http://www.r-bloggers.com/r-gpu-programming-for-all-with-gpur/" data-counturl="http://www.r-bloggers.com/r-gpu-programming-for-all-with-gpur/" data-count="horizontal" data-text="R – GPU Programming for All with ‘gpuR’" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.parallelr.com/r-gpu-programming-for-all-with-gpur/"> Blog – ParallelR</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<h3>Introduction</h3>
<p>GPUs (Graphic Processing Units) have become much more popular in recent years for computationally intensive calculations.  Despite these gains, the use of this hardware has been very limited in the R programming language.  Although possible, the prospect of programming in either OpenCL or CUDA is difficult for many programmers unaccustomed to working with such a low-level interface.  Creating bindings for R’s high-level programming that abstracts away the complex GPU code would make using GPUs far more accessible to R users.  This is the core idea behind the <a target="_blank" href="https://cran.r-project.org/web/packages/gpuR/index.html">gpuR</a> package.  There are three novel aspects of <code>gpuR</code>:</p>
<ol>
<li>Applicable on ‘ALL’ GPUs</li>
<li>Abstracts away CUDA/OpenCL code to easily incorporate in to existing R algorithms</li>
<li>Separates copy/compute functions to allow objects to persist on GPU</li>
</ol>
<h3><strong>Broad application</strong>:</h3>
<p>The ‘gpuR’ package was created to bring the power of GPU computing to any R user with a GPU device.  Although there are a handful of packages that provide some GPU capability (e.g. <a target="_blank" href="https://cran.r-project.org/web/packages/gputools/index.html">gputools</a>, <a target="_blank" href="http://cran.fhcrc.org/web/packages/cudaBayesreg/index.html">cudaBayesreg</a>, <a target="_blank" href="https://cran.r-project.org/web/packages/HiPLARM/index.html">HiPLARM</a>, <a target="_blank" href="https://cran.r-project.org/web/packages/HiPLARb/index.html">HiPLARb</a>, and <a target="_blank" href="https://cran.r-project.org/web/packages/gmatrix/index.html">gmatrix</a>) all are strictly limited to NVIDIA GPUs.  As such, a backend that is based upon OpenCL would allow all users to benefit from GPU hardware.  The ‘gpuR’ package therefore utilizes the <a target="_blank" href="http://viennacl.sourceforge.net/">ViennaCL</a> linear algebra library which contains auto-tuned OpenCL kernels (among others) that can be leveraged for GPUs.  The headers have been conveniently repackaged in the <a target="_blank" href="https://cran.r-project.org/web/packages/RViennaCL/index.html">RViennaCL </a>package.  It also allows for a CUDA backend for those with NVIDIA GPUs that may see further improved performance (contained within the companion <a target="_blank" href="https://github.com/cdeterman/gpuRcuda">gpuRcuda</a> package not yet formally released).</p>
<h3><strong>Abstract away GPU code:</strong></h3>
<p>The <code>gpuR</code> package uses the S4 object oriented system to have explicit classes and methods that all the user to simply cast their <code>matrix</code> or <code>vector</code> and continue programming in R as normal.  For example:</p>
<pre>
ORDER = 1024

A = matrix(rnorm(ORDER^2), nrow=ORDER)
B = matrix(rnorm(ORDER^2), nrow=ORDER)
gpuA = gpuMatrix(A, type=&quot;double&quot;)
gpuB = gpuMatrix(B, type=&quot;double&quot;)

C = A %*% B
gpuC = gpuA %*% gpuB

all(C == gpuC[])
[1] TRUE
</pre>
<p>The <code>gpuMatrix</code> object points to a matrix in RAM which is then computed by the GPU when a desired function is called.  This avoids R’s habit of copying the memory of objects.  For example:</p>
<pre>

library(pryr)

# Initially points to same object
x = matrix(rnorm(16), 4)
y = x

address(x)
[1] &quot;0x16177f28&quot;

address(y)
[1] &quot;0x16177f28&quot;

# But once modify the second object it creates a copy
y[1,1] = 0

address(x)
[1] &quot;0x16177f28&quot;

address(y)
[1] &quot;0x15fbb1d8


</pre>
<p>In contrast, the same syntax for a <code>gpuMatrix</code> will modify the original object in-place without any copy.</p>
<pre>

library(pryr)
library(gpuR)

# Initially points to same object
x = gpuMatrix(rnorm(16), 4, 4)
y = x

x@address
[1] &lt;pointer: 0x6baa040>

y@address
[1] &lt;pointer: 0x6baa040>

# Modification affects both objects without copy
y[1,1] = 0

x@address
[1] &lt;pointer: 0x6baa040>

y@address
[1] &lt;pointer: 0x6baa040>
</pre>
<p> </p>
<p>Each new variable assigned to this object will only copy the pointer thereby making the program more memory efficient.  However, the <code>gpuMatrix></code> class does still require allocating GPU memory and copying data to device for each function call. The most commonly used methods have been overloaded such as  %*%, +, -, *, /, crossprod, tcrossprod, and trig functions among others.  In this way, an R user can create these objects and leverage GPU resources without the need to know a bunch more functions that would break existing algorithms.</p>
<h3><strong>Distinct Copy/Compute Functionality:</strong></h3>
<p>For the <code>gpuMatix</code> and <code>gpuVector</code> classes there are companion <code>vclMatrix</code> and <code>vclVector</code> class that point to objects that persist in the GPU RAM.  In this way, the user explicitly decides when data needs to be moved back to the host.  By avoiding unnecessary data transfer between host and device performance can significantly improve.  For example:</p>
<pre>

vclA = vclMatrix(rnorm(10000), nrow = 100)
vclB = vclMatrix(rnorm(10000), nrow = 100)
vclC = vclMatrix(rnorm(10000), nrow = 100)

# GEMM
vclD = vclA %*% vclB

# Element-wise addition
vclD = vclD + vclC

</pre>
<p>In this code, the three initial matrices already exist in the GPU memory so no data transfer takes place in the GEMM call.  Furthermore, the returned matrix remains in the GPU memory.  In this case, the ‘vclD’ object is still in GPU RAM. As such, the element-wise addition call also happens directly on the GPU with no data transfers. It is worth also noting that the user can still modify elements, rows, or columns with the exact same syntax as a normal R matrix.</p>
<pre>
vclD[1,1] = 42
vclD[,2] = rep(12, 100)
vclD[3,] = rep(23, 100)
</pre>
<p>These operations simply copy the <em>new</em> elements to the GPU and modify the object in-place within the GPU memory. The ‘vclD’ object is never copied to the host.</p>
<h3>Benchmarks:</h3>
<p>With all that in mind, how does gpuR perform?  Here are some general benchmarks of the popular GEMM operation.  I currently only have access to a single NVIDIA GeForce GTX 970 for these simulations.  Users should expect to see differences with high performance GPUs (e.g. AMD FirePro, NVIDIA Tesla, etc.). Speedup relative to CPU will also vary depending upon user hardware.</p>
<h4>(1) Default dGEMM vs Base R</h4>
<p>R is known to only support two numeric types (integer and double).  As such, Figure 1 shows the fold speedup achieved by using the <code>gpuMatrix</code> and <code>vclMatrix</code> classes.  Since R is already known to not be the fastest language, an implementation with the OpenBLAS backend is included as well for reference using a 4 core Intel i5-2500 CPU @ 3.30GHz.  As can be seen there is a dramatic speedup from just using OpenBLAS or the gpuMatrix class (essentially equivalent).  Of interest is the impact of the transfer time from host-device-host that is typical in many GPU implementations.  This cost is eliminated by using the <code>vclMatrix</code> class which continues to scale with matrix size.</p>
<div class="wp-caption"><img sizes="(max-width: 640px) 100vw, 640px" srcset="http://i2.wp.com/www.parallelr.com/wp-content/uploads/2016/05/dgemm-1024x588.png?resize=456%2C262 1024w, http://www.parallelr.com/wp-content/uploads/2016/05/dgemm-300x172.png 300w, http://www.parallelr.com/wp-content/uploads/2016/05/dgemm-768x441.png 768w" alt="dgemm" class="wp-image-768" src="http://i2.wp.com/www.parallelr.com/wp-content/uploads/2016/05/dgemm-1024x588.png?resize=456%2C262" data-recalc-dims="1">
<p class="wp-caption-text">Figure 1 – Fold speedup achieved using openblas (CPU) as well as the gpuMatrix/vclMatrix (GPU) classes provided in gpuR.</p>
</div>
<p> </p>
<h4>(2) sGEMM vs Base R</h4>
<p>In many GPU benchmarks there is often float operations measured as well.  As noted above, R does not provide this by default.  One way to go around this is to use the <a target="_blank" href="https://cran.r-project.org/web/packages/RcppArmadillo/index.html">RcppArmadillo</a> package and explicitly casting R objects as float types.  The armadillo library will also default to using the BLAS backend provided (i.e. OpenBLAS).  Figure 2 shows the impact of using float data types.  OpenBLAS continues to provide a noticeable speedup but <code>gpuMatrix</code> begins to outperform once matrix order exceeds 1500.  The <code>vclMatrix</code> continues to demonstrate the value in retaining objects in GPU memory and avoiding memory transfers.</p>
<p> </p>
<div class="wp-caption"><img sizes="(max-width: 640px) 100vw, 640px" srcset="http://i1.wp.com/www.parallelr.com/wp-content/uploads/2016/05/sgemm-1024x588.png?resize=456%2C262 1024w, http://www.parallelr.com/wp-content/uploads/2016/05/sgemm-300x172.png 300w, http://www.parallelr.com/wp-content/uploads/2016/05/sgemm-768x441.png 768w" alt="sgemm" class="wp-image-769" src="http://i1.wp.com/www.parallelr.com/wp-content/uploads/2016/05/sgemm-1024x588.png?resize=456%2C262" data-recalc-dims="1">
<p class="wp-caption-text">Figure 2 – Float type GEMM comparisons. Fold speedup achieved using openblas (via RcppArmadillo) as well as the gpuMatrix/vclMatrix (GPU) classes provided in gpuR.</p>
</div>
<p> </p>
<p>To give an additional view on the performance achieved by <code>gpuMatrix</code> and <code>vclMatrix</code> is comparing directly against the OpenBLAS performance.  The <code>gpuMatrix</code> reaches ~2-3 fold speedup over OpenBLAS whereas <code>vclMatrix</code> scales to over 100 fold speedup!  It is curious as to why the performance with vcl<code>Matrix</code> is so much faster (only differing in host-device-host transfers).  Further optimization with <code>gpuMatrix</code> will need to be explored (fresh eyes are welcome) accepting limitations in the BUS transfer speed.  Performance will certainly improve with improved hardware capabilities such as NVIDIA’s NVLink.</p>
<div class="wp-caption"><img sizes="(max-width: 737px) 100vw, 737px" srcset="http://i1.wp.com/www.parallelr.com/wp-content/uploads/2016/05/sgemm_openblas-1024x528.png?resize=456%2C235 1024w, http://www.parallelr.com/wp-content/uploads/2016/05/sgemm_openblas-300x155.png 300w, http://www.parallelr.com/wp-content/uploads/2016/05/sgemm_openblas-768x396.png 768w" alt="sgemm_openblas" class="wp-image-831" src="http://i1.wp.com/www.parallelr.com/wp-content/uploads/2016/05/sgemm_openblas-1024x528.png?resize=456%2C235" data-recalc-dims="1">
<p class="wp-caption-text">Figure 3 – Fold speedup achieved over openblas (via RcppArmadillo) float type GEMM comparisons vs the gpuMatrix/vclMatrix (GPU) classes provided in gpuR.</p>
</div>
<h3>Conclusion</h3>
<p>The <code>gpuR</code> package has been created to bring GPU computing to as many R users as possible.  It is the intention to use <code>gpuR</code> to more easily supplement current and future algorithms that could benefit from GPU acceleration.  The <code>gpuR</code> package is currently available on <a target="_blank" href="https://cran.r-project.org/web/packages/gpuR/index.html">CRAN</a>.  The development version can be found on my <a target="_blank" href="https://github.com/cdeterman/gpuR">github</a> in addition to existing issues and wiki pages (assisting primarily in installation).  Future developments include solvers (e.g. QR, SVD, cholesky, etc.), scaling across multiple GPUs,  ‘sparse’ class objects, and custom OpenCL kernels.</p>
<p>As noted above, this package is intended to be used with a multitude of hardware and operating systems (it has been tested on Windows, Mac, and multiple Linux flavors).  I only have access to a limited set of hardware (I can’t access every GPU, let along the most expensive).  As such, the development of <code>gpuR</code> depends upon the R user community.  Volunteers who possess different hardware are always welcomed and encouraged to submit issues regarding any discovered bugs.  I have begun a gitter account for users to report on successful usage with alternate hardware.  Suggestions and general conversation about gpuR is welcome.</p>
<p> </p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/r-gpu-programming-for-all-with-gpur/"></div></div><div><a data-url="http://www.r-bloggers.com/r-gpu-programming-for-all-with-gpur/" data-counturl="http://www.r-bloggers.com/r-gpu-programming-for-all-with-gpur/" data-count="vertical" data-text="R – GPU Programming for All with ‘gpuR’" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.parallelr.com/r-gpu-programming-for-all-with-gpur/"> Blog – ParallelR</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/JnDxCncYSro"></a></div>