[Rd] Fw: Calling a LAPACK subroutine from R

[Serguei Sokol]

On 11/09/2019 21:38, Berend Hasselman wrote:
> 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>
another way is to use directly dyn.load():

lapack.path <- paste0(file.path(R.home(), ifelse(.Platform$OS.type == 
"windows",
          file.path("bin", .Platform$r_arch, "Rlapack"), 
file.path("lib", "libRlapack"))),
          .Platform$dynlib.ext)
dyn.load(lapack.path)

followed by your code.

Best,
Serguei.

>
> 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
> ______________________________________________
> R-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>