title: Without root access, run R with tuned BLAS when it is linked with reference BLAS
body: <p>Can any one tell me why I can <strong>not successfully</strong> test OpenBLAS's <code>dgemm</code> performance (in GFLOPs) in R via the following way?</p>
  
  <ol>
  <li>link R with the "reference BLAS" <code>libblas.so</code></li>
  <li>compile my C program <code>mmperf.c</code> with OpenBLAS library <code>libopenblas.so</code></li>
  <li>load the resulting shared library <code>mmperf.so</code> into R, call the R wrapper function <code>mmperf</code> and report <code>dgemm</code> performance in GFLOPs.</li>
  </ol>
  
  <p>Point 1 looks strange, but I have no choice because I have no root access on machines I want to test, so actual linking to OpenBLAS is impossible. By <strong>"not successfully"</strong> I mean my program ends up reporting <code>dgemm</code> performance for reference BLAS instead of OpenBLAS. I hope someone can explain to me:</p>
  
  <ol>
  <li>why my way does not work;</li>
  <li>is it possible at all to make it work (this is important, because if it is impossible, I must write a C <code>main</code> function and do my job in a C program.)</li>
  </ol>
  
  <p>I've investigated into this issue for two days, here I will include various system output to assist you to make a diagnose. To make things reproducible, I will also include the code, makefile as well as shell command.</p>
  
  <p><strong>Part 1: system environment before testing</strong></p>
  
  <p>There are 2 ways to invoke R, either using <code>R</code> or <code>Rscript</code>. There are some differences in what is loaded when they are invoked:</p>
  
  <pre><code>~/Desktop/dgemm$ readelf -d $(R RHOME)/bin/exec/R | grep "NEEDED"
  0x00000001 (NEEDED)         Shared library: [libR.so]
  0x00000001 (NEEDED)         Shared library: [libpthread.so.0]
  0x00000001 (NEEDED)         Shared library: [libc.so.6]
  
  ~/Desktop/dgemm$ readelf -d $(R RHOME)/bin/Rscript | grep "NEEDED"
  0x00000001 (NEEDED)         Shared library: [libc.so.6]
  </code></pre>
  
  <p>Here we need to choose <code>Rscript</code>, because <code>R</code> loads <code>libR.so</code>, which will automatically load the reference BLAS <code>libblas.so.3</code>:</p>
  
  <pre><code>~/Desktop/dgemm$ readelf -d $(R RHOME)/lib/libR.so | grep blas
  0x00000001 (NEEDED)         Shared library: [libblas.so.3]
  
  ~/Desktop/dgemm$ ls -l /etc/alternatives/libblas.so.3
  ... 31 May /etc/alternatives/libblas.so.3 -&gt; /usr/lib/libblas/libblas.so.3.0
  
  ~/Desktop/dgemm$ readelf -d /usr/lib/libblas/libblas.so.3 | grep SONAME
  0x0000000e (SONAME)         Library soname: [libblas.so.3]
  </code></pre>
  
  <p>Comparatively, <code>Rscript</code> gives a cleaner environment.</p>
  
  <p><strong>Part 2: OpenBLAS</strong></p>
  
  <p>After downloading source file from <a href="http://www.openblas.net/" rel="noreferrer">OpenBLAS</a> and a simple <code>make</code> command, a shared library of the form <code>libopenblas-&lt;arch&gt;-&lt;release&gt;.so-&lt;version&gt;</code> can be generated. Note that we will not have root access to install it; instead, we copy this library into our working directory <code>~/Desktop/dgemm</code> and rename it simply to <code>libopenblas.so</code>. At the same time we have to make another copy with name <code>libopenblas.so.0</code>, as this is the <em>SONAME</em> which run time loader will seek for:</p>
  
  <pre><code>~/Desktop/dgemm$ readelf -d libopenblas.so | grep "RPATH\|SONAME"
  0x0000000e (SONAME)         Library soname: [libopenblas.so.0]
  </code></pre>
  
  <p>Note that the <code>RPATH</code> attribute is not given, which means this library is intended to be put in <code>/usr/lib</code> and we should call <code>ldconfig</code> to add it to <code>ld.so.cache</code>. But again we don't have root access to do this. In fact, if this can be done, then all the difficulties are gone. We could then use <code>update-alternatives --config libblas.so.3</code> to effectively link R to OpenBLAS.</p>
  
  <p><strong>Part 3: C code, Makefile, and R code</strong></p>
  
  <p>Here is a C script <code>mmperf.c</code> computing GFLOPs of multiplying 2 square matrices of size <code>N</code>:</p>
  
  <pre><code>#include &lt;R.h&gt;
  #include &lt;Rmath.h&gt;
  #include &lt;Rinternals.h&gt;
  #include &lt;R_ext/BLAS.h&gt;
  #include &lt;sys/time.h&gt;
  
  /* standard C subroutine */
  double mmperf (int n) {
    /* local vars */
    int n2 = n * n, tmp; double *A, *C, one = 1.0;
    struct timeval t1, t2; double elapsedTime, GFLOPs;
    /* simulate N-by-N matrix A */
    A = (double *)calloc(n2, sizeof(double));
    GetRNGstate();
    tmp = 0; while (tmp &lt; n2) {A[tmp] = runif(0.0, 1.0); tmp++;}
    PutRNGstate();
    /* generate N-by-N zero matrix C */
    C = (double *)calloc(n2, sizeof(double));
    /* time 'dgemm.f' for C &lt;- A * A + C */
    gettimeofday(&amp;t1, NULL);
    F77_CALL(dgemm) ("N", "N", &amp;n, &amp;n, &amp;n, &amp;one, A, &amp;n, A, &amp;n, &amp;one, C, &amp;n);
    gettimeofday(&amp;t2, NULL);
    /* free memory */
    free(A); free(C);
    /* compute and return elapsedTime in microseconds (usec or 1e-6 sec) */
    elapsedTime = (double)(t2.tv_sec - t1.tv_sec) * 1e+6;
    elapsedTime += (double)(t2.tv_usec - t1.tv_usec);
    /* convert microseconds to nanoseconds (1e-9 sec) */
    elapsedTime *= 1e+3;
    /* compute and return GFLOPs */
    GFLOPs = 2.0 * (double)n2 * (double)n / elapsedTime;
    return GFLOPs;
    }
  
  /* R wrapper */
  SEXP R_mmperf (SEXP n) {
    double GFLOPs = mmperf(asInteger(n));
    return ScalarReal(GFLOPs);
    }
  </code></pre>
  
  <p>Here is a simple R script <code>mmperf.R</code> to report GFLOPs for case <code>N = 2000</code></p>
  
  <pre><code>mmperf &lt;- function (n) {
    dyn.load("mmperf.so")
    GFLOPs &lt;- .Call("R_mmperf", n)
    dyn.unload("mmperf.so")
    return(GFLOPs)
    }
  
  GFLOPs &lt;- round(mmperf(2000), 2)
  cat(paste("GFLOPs =",GFLOPs, "\n"))
  </code></pre>
  
  <p>Finally there is a simple makefile to generate the shared library <code>mmperf.so</code>:</p>
  
  <pre><code>mmperf.so: mmperf.o
      gcc -shared -L$(shell pwd) -Wl,-rpath=$(shell pwd) -o mmperf.so mmperf.o -lopenblas
  
  mmperf.o: mmperf.c
      gcc -fpic -O2 -I$(shell Rscript --default-packages=base --vanilla -e 'cat(R.home("include"))') -c mmperf.c
  </code></pre>
  
  <p>Put all these files under working directory <code>~/Desktop/dgemm</code>, and compile it:</p>
  
  <pre><code>~/Desktop/dgemm$ make
  ~/Desktop/dgemm$ readelf -d mmperf.so | grep "NEEDED\|RPATH\|SONAME"
  0x00000001 (NEEDED)            Shared library: [libopenblas.so.0]
  0x00000001 (NEEDED)            Shared library: [libc.so.6]
  0x0000000f (RPATH)             Library rpath: [/home/zheyuan/Desktop/dgemm]
  </code></pre>
  
  <p>The output reassures us that OpenBLAS is correctly linked, and the run time load path is correctly set.</p>
  
  <p><strong>Part 4: testing OpenBLAS in R</strong></p>
  
  <p>Let's do</p>
  
  <pre><code>~/Desktop/dgemm$ Rscript --default-packages=base --vanilla mmperf.R
  </code></pre>
  
  <p>Note our script needs only the <code>base</code> package in R, and <code>--vanilla</code> is used to ignore all user settings on R start-up. On my laptop, my program returns:</p>
  
  <pre><code>GFLOPs = 1.11
  </code></pre>
  
  <p>Oops! This is truely reference BLAS performance not OpenBLAS (which is about 8-9 GFLOPs).</p>
  
  <p><strong>Part 5: Why?</strong></p>
  
  <p>To be honest, I don't know why this happens. Each step seems to work correctly. Does something subtle occurs when R is invoked? For example, any possibility that OpenBLAS library is overridden by reference BLAS at some point for some reason? Any explanations and solutions? Thanks!</p>