[Rd] Calling a LAPACK subroutine from R

Berend Hasselman bhh @end|ng |rom x@4@||@n|
Thu Sep 12 08:57:02 CEST 2019 

I have tried what I proposed in a virtual Kubuntu 18.04 which uses gfortran 7.4.
I used the latest development version of R.

It worked just as on macOS.

Berend


> On 11 Sep 2019, at 22:07, Göran Broström <goran.brostrom using umu.se> wrote:
> 
> Berend,
> 
> I do not think this works with gfortran 7+. I am calling the BLAS subroutine dgemv from Fortran code in my package eha, and the check (with R-devel) gives:
> 
> gmlfun.f:223:1: warning: type of ‘dgemv’ does not match original declaration [-Wlto-type-mismatch]
>      &     score, ione)
> ^
> /home/gobr0002/R/src/R-devel/include/R_ext/BLAS.h:107:1: note: type mismatch in parameter 12
> F77_NAME(dgemv)(const char *trans, const int *m, const int *n,
> 
> Type of a Fortran subroutine is matched against type of a C function?! My conclusion is that it is impossible to call a BLAS subroutine with a character parameter from Fortran code (nowadays). Calling from C code is fine, on the other hand(!).
> 
> I have recently asked about this on R-pkg-devel, but not received any useful answers, and my submission to CRAN is rejected. I solve it by making a personal copy of dgemv and changing the character parameter to integer, and adding Jack Dongarra, Jeremy Du Croz, Sven Hammarling, and Richard Hanson as authors of eha. And a Copyright note, all in the DESCRIPTION file. Ugly but what can I do (except rewriting the Fortran code in C with f2c)?
> 
> Göran
> 
> On 2019-09-11 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>
>> 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

More information about the R-devel mailing list