[Rd] Fw: Calling a LAPACK subroutine from R

[Berend Hasselman]

The Lapack library is loaded automatically by R itself when it needs it  for doing some calculation.
You can force it to do that with a (dummy) solve for example.
Put this at start of your script:

<code>
# dummy code to get LAPACK library loaded
X1 <- diag(2,2)
x1 <- rep(2,2)
# X1;x1
z <- solve(X1,x1)
</code>

followed by the rest of your script.
You will get a warning (I do) that  "passing a character vector  to .Fortran is not portable".
On other systems this may gave fatal errors. This is quick and very dirty. Don't do it.

I believe there is a better and much safer way to achieve what you want.
Here goes.

Create a folder (directory) src in the directory where your script resides.
Create a wrapper for "dpbtrf" file in a file xdpbtrf.f that takes an integer instead of character

<xdpbtrf.f>
c intermediate for dpbtrf

      SUBROUTINE xDPBTRF( kUPLO, N, KD, AB, LDAB, INFO )

c      .. Scalar Arguments ..
      integer         kUPLO
      INTEGER         INFO, KD, LDAB, N

c  .. Array Arguments ..
      DOUBLE PRECISION   AB( LDAB, * )

      character UPLO
c     convert integer argument to character
      if(kUPLO .eq. 1 ) then
          UPLO = 'L'
      else
          UPLO = 'U'
      endif

      call dpbtrf(UPLO,N,KD,AB,LDAB,INFO)
      return
      end
</xdpbtrf.f>


Instead of a character argument UPLO it takes an integer argument kUPLO.
The meaning should be obvious from the code.

Now create a shell script in the folder of your script to generate a dynamic library to be loaded in your script:

<mkso.sh>
# Build a binary dynamic library for accessing Lapack dpbtrf

# syntax checking
 
SONAME=xdpbtrf.so

echo Strict syntax checking
echo ----------------------
gfortran -c -fsyntax-only -fimplicit-none -Wall src/*.f || exit 1

LAPACK=$(R CMD config LAPACK_LIBS)
R CMD SHLIB --output=${SONAME} src/*.f ${LAPACK} || exit 1
</mkso.sh>

To load the dynamic library xdpbtrf.so  change your script into this

<yourscript>
dyn.load("xdpbtrf.so")
n <- 4L
phi <- 0.64
AB <- matrix(0, 2, n)
AB[1, ] <- c(1, rep(1 + phi^2, n-2), 1)
AB[2, -n] <- -phi
round(AB, 3)

AB.ch <- .Fortran("xdpbtrf", kUPLO=1L, N = as.integer(n),
                            KD = 1L, AB = AB, LDAB = 2L, INFO = as.integer(0))$AB
AB.ch

</yourscript>

and you are good to go.

You should always do something  as described above when you need to pass character arguments to Fortran code.

All of this was tested and run on macOS using the CRAN version of R.

Berend Hasselman

> On 11 Sep 2019, at 15:47, Giovanni Petris <gpetris using uark.edu> wrote:
> 
> Sorry for cross-posting, but I realized my question might be more appropriate for r-devel...
> 
> Thank you,
> Giovanni
> 
> ________________________________________
> From: R-help <r-help-bounces using r-project.org> on behalf of Giovanni Petris <gpetris using uark.edu>
> Sent: Tuesday, September 10, 2019 16:44
> To: r-help using r-project.org
> Subject: [R] Calling a LAPACK subroutine from R
> 
> Hello R-helpers!
> 
> I am trying to call a LAPACK subroutine directly from my R code using .Fortran(), but R cannot find the symbol name. How can I register/load the appropriate library?
> 
>> ### AR(1) Precision matrix
>> n <- 4L
>> phi <- 0.64
>> AB <- matrix(0, 2, n)
>> AB[1, ] <- c(1, rep(1 + phi^2, n-2), 1)
>> AB[2, -n] <- -phi
>> round(AB, 3)
>      [,1]  [,2]  [,3] [,4]
> [1,]  1.00  1.41  1.41    1
> [2,] -0.64 -0.64 -0.64    0
>> 
>> ### Cholesky factor
>> AB.ch <- .Fortran("dpbtrf", UPLO = 'L', N = as.integer(n),
> +                  KD = 1L, AB = AB, LDAB = 2L, INFO = as.integer(0))$AB
> Error in .Fortran("dpbtrf", UPLO = "L", N = as.integer(n), KD = 1L, AB = AB,  :
>  Fortran symbol name "dpbtrf" not in load table
>> sessionInfo()
> R version 3.6.0 (2019-04-26)
> Platform: x86_64-apple-darwin18.5.0 (64-bit)
> Running under: macOS Mojave 10.14.6
> 
> Matrix products: default
> BLAS/LAPACK: /usr/local/Cellar/openblas/0.3.6_1/lib/libopenblasp-r0.3.6.dylib
> 
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
> 
> loaded via a namespace (and not attached):
> [1] compiler_3.6.0 tools_3.6.0
> 
> Thank you in advance for your help!
> 
> Best,
> Giovanni Petris
> 
> 
> 
> --
> Giovanni Petris, PhD
> Professor
> Director of Statistics
> Department of Mathematical Sciences
> University of Arkansas - Fayetteville, AR 72701
> 
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Dhelp&d=DwICAg&c=7ypwAowFJ8v-mw8AB-SdSueVQgSDL4HiiSaLK01W8HA&r=C3DNvy_azplKSvJKgvsgjA&m=C-MwKl__0xz-98RBbu7QNXJjqWkRr4xp6c0cz9Dck7A&s=a1vAu3mcXKObTLwP19vOmRPq55h6oQTh_vnS6BEibF0&e=
> PLEASE do read the posting guide https://urldefense.proofpoint.com/v2/url?u=http-3A__www.R-2Dproject.org_posting-2Dguide.html&d=DwICAg&c=7ypwAowFJ8v-mw8AB-SdSueVQgSDL4HiiSaLK01W8HA&r=C3DNvy_azplKSvJKgvsgjA&m=C-MwKl__0xz-98RBbu7QNXJjqWkRr4xp6c0cz9Dck7A&s=qFGlplF9cOSmnDUvugsPRDn4iZS7v-LuWNAvfY69sbA&e=
> and provide commented, minimal, self-contained, reproducible code.
> 
> ______________________________________________
> R-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel