[Rd] Fw: Calling a LAPACK subroutine from R
Göran Broström
gor@n@bro@trom @end|ng |rom umu@@e
Thu Sep 12 11:56:05 CEST 2019
Hi guys,
interestingly, my problem seems to be solved by writing a FORTRAN
wrapper for the Fortran code! (As long as the check doesn't get
smarter...). This is the relevant part of my Fortran code:
-----------------------------------------------------------
subroutine gmlfun(what,
& totevent, totrs, ns,
& antrs, antevents, size,
& totsize, eventset, riskset,
& nn, antcov, covar, offset,
& beta, gamma,
& loglik, h1, h2, h11, h21, h22,
& score)
......
call gmlfun1(what,
& totevent, totrs, ns,
& antrs, antevents, size,
& totsize, eventset, riskset,
& nn, antcov, covar, offset,
& beta, gamma,
& loglik, h1, h2, h11, h21, h22,
& score)
return
end
C ***
C
subroutine gmlfun1(what,
& totevent, totrs, ns,
& antrs, antevents, size,
& totsize, eventset, riskset,
& nn, antcov, covar, offset,
& beta, gamma,
& loglik, h1, h2, h11, h21, h22,
& score)
.....
C
call dcopy(nn, offset, ione, score, ione)
call dgemv(trans, nn, antcov, one, covar, nn, beta, ione, one,
& score, ione)
.....
return
end
-----------------------------------------------------------------
gmlfun is called directly by .Fortran, and in my original code gmlfun
was gmlfun1. Now gmlfun is just a wrapper that calls gmlfun1, which
calls dgemv, and no complaints!
This is apparently what Berend (and now myself) gets away with: Quite
bizarre in my mind.
Or: Is this an R-blessed solution?
Thanks, Göran
On 2019-09-12 11:25, Tomas Kalibera wrote:
> On 9/12/19 11:07 AM, Göran Broström wrote:
>> Kurt, see below:
>>
>> On 2019-09-12 08:42, Kurt Hornik wrote:
>>>>>>>> Göran Broström writes:
>>>
>>> Göran,
>>>
>>> Pls allow me to join the discussions on your pending CRAN submission.
>>>
>>> First, the current solution with the copy of the BLAS sources is not
>>> quite perfect. You now have
>>>
>>> Authors using R: c(person("Göran", "Broström", role = c("aut", "cre"),
>>> email = "goran.brostrom using umu.se"),
>>> person("Jack", "Dongarra", role = "aut"),
>>> person("Jeremy", "Du Croz", role = "aut"),
>>> person("Sven", "Hammarling", role = "aut"),
>>> person("Richard", "Hanson", role = "aut"),
>>> person("Jianming", "Jin", role = "aut"))
>>> Copyright: The blas Developers 2014-2018.
>>>
>>> I think the last 5 persons should get a ctb role instead of aut, and you
>>> should add
>>>
>>> person("The blas Developers", role = "cph", comment = "....")
>>>
>>> with a suitable comment explaining what the copyright is held for.
>>>
>>> However, I would still hope that you can do without the copy by
>>> following the recipe in section "Fortran character strings" in "Writing
>>> R Extensions". This says:
>>>
>>> ************************************************************************
>>> Alternatively, do as R does as from version 3.6.1-patched and pass the
>>> character length(s) from C to Fortran. A portable way to do this is
>>> // before any R headers, or define in PKG_CPPFLAGS
>>> #define USE_FC_LEN_T
>>> #include <Rconfig.h>
>>> #include <R_ext/BLAS.h>
>>> #ifndef FCONE
>>> # define FCONE
>>> #endif
>>> ...
>>> F77_CALL(dgemm)("N", "T", &nrx, &ncy, &ncx, &one, x,
>>> &nrx, y, &nry, &zero, z, &nrx FCONE
>>> FCONE);
>>> (Note there is no comma before or between the 'FCONE' invocations.) It
>>> is strongly recommended that packages which call from C/C++ BLAS/LAPACK
>>> routines with character arguments adopt this approach.
>>> ************************************************************************
>>>
>>> Does this really not work for you?
>>
>> Of course it would work: I call dgemv (and other BLAS routines)
>> frequently from C code in my package, that is not the problem. I am
>> just wondering why Berend gets away with direct calls to dgemv from
>> Fortran without using the 12th hidden parameter FCLEN (FCONE), when I
>> fail. I suspect an oversight in the R checking system.
>
> The hidden argument will be added by the Fortran compiler, at the call
> site. If BLAS is built with the same Fortran compiler as your Fortran
> code which calls it, there should be no problem (often it will work also
> when the compilers are different, but that depends on how different they
> are). If you can't find why it is not working for you, please create a
> minimum reproducible but complete example. Please try looking also at
> existing packages with Fortran code, including "stats" for instance.
> There are several people who are trying to help you on the mailing
> lists, it would be easier for them if they knew exactly what you were
> doing.
>
> Best
> Tomas
>
>
>>
>> I can write a C wrapper for my Fortran function that calls dgemv and
>> call it in R by .C, but is that really guaranteed to work, when the
>> checking routine gets sharper?
>
>>
>> Best, G,
>>
>>> Best
>>> -k
>>>
>>>
>>>> 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
>>>>>
>>>
>>>> ______________________________________________
>>>> R-devel using r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>
More information about the R-devel
mailing list