Fw: Calling a LAPACK subroutine from R

classic Classic list List threaded Threaded
13 messages Options
Reply | Threaded
Open this post in threaded view
|

Fw: Calling a LAPACK subroutine from R

Giovanni Petris
Sorry for cross-posting, but I realized my question might be more appropriate for r-devel...

Thank you,
Giovanni

________________________________________
From: R-help <[hidden email]> on behalf of Giovanni Petris <[hidden email]>
Sent: Tuesday, September 10, 2019 16:44
To: [hidden email]
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


______________________________________________
[hidden email] 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.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[hidden email]> wrote:
>
> Sorry for cross-posting, but I realized my question might be more appropriate for r-devel...
>
> Thank you,
> Giovanni
>
> ________________________________________
> From: R-help <[hidden email]> on behalf of Giovanni Petris <[hidden email]>
> Sent: Tuesday, September 10, 2019 16:44
> To: [hidden email]
> 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
>
>
> ______________________________________________
> [hidden email] 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.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Calling a LAPACK subroutine from R

Göran Broström-3
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 <[hidden email]> wrote:
>>
>> Sorry for cross-posting, but I realized my question might be more appropriate for r-devel...
>>
>> Thank you,
>> Giovanni
>>
>> ________________________________________
>> From: R-help <[hidden email]> on behalf of Giovanni Petris <[hidden email]>
>> Sent: Tuesday, September 10, 2019 16:44
>> To: [hidden email]
>> 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
>>
>>
>> ______________________________________________
>> [hidden email] 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.
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Calling a LAPACK subroutine from R

Avraham Adler
Can you write a small C function that calls LAPACK call that fro your
Fortran code? Yes, an extra step but maybe less traumatic than rewriting
parts of LAPACK directly.

Avi

On Wed, Sep 11, 2019 at 4:08 PM Göran Broström <[hidden email]>
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 <[hidden email]> wrote:
> >>
> >> Sorry for cross-posting, but I realized my question might be more
> appropriate for r-devel...
> >>
> >> Thank you,
> >> Giovanni
> >>
> >> ________________________________________
> >> From: R-help <[hidden email]> on behalf of Giovanni
> Petris <[hidden email]>
> >> Sent: Tuesday, September 10, 2019 16:44
> >> To: [hidden email]
> >> 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
> >>
> >>
> >> ______________________________________________
> >> [hidden email] 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.
> >>
> >> ______________________________________________
> >> [hidden email] mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-devel
> >
> > ______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
> >
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
--
Sent from Gmail Mobile

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Calling a LAPACK subroutine from R

Göran Broström-3


On 2019-09-11 22:16, Avraham Adler wrote:
> Can you write a small C function that calls LAPACK call that fro your
> Fortran code? Yes, an extra step but maybe less traumatic than rewriting
> parts of LAPACK directly.

Yes, I know how to do that, but I find it somewhat bizarre that it is
impossible to call a Fortran subroutine from Fortran. And rewriting
'dgemv' was simple: Just change character to integer and 'N' to 1. And
rename the subroutine. The hard (tedious) part was to include all the
LAPACK authors in my DESCRIPTION file.

My guess is that the root cause is that BLAS/LAPACK is written in
FORTRAN 77, which is said to be a subset of the current Fortran version
but obviously isn't.

Thanks, Göran

>
> Avi
>
> On Wed, Sep 11, 2019 at 4:08 PM Göran Broström <[hidden email]
> <mailto:[hidden email]>> 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 <[hidden email]
>     <mailto:[hidden email]>> wrote:
>      >>
>      >> Sorry for cross-posting, but I realized my question might be
>     more appropriate for r-devel...
>      >>
>      >> Thank you,
>      >> Giovanni
>      >>
>      >> ________________________________________
>      >> From: R-help <[hidden email]
>     <mailto:[hidden email]>> on behalf of Giovanni Petris
>     <[hidden email] <mailto:[hidden email]>>
>      >> Sent: Tuesday, September 10, 2019 16:44
>      >> To: [hidden email] <mailto:[hidden email]>
>      >> 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
>      >>
>      >>
>      >> ______________________________________________
>      >> [hidden email] <mailto:[hidden email]> 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.
>      >>
>      >> ______________________________________________
>      >> [hidden email] <mailto:[hidden email]> mailing list
>      >> https://stat.ethz.ch/mailman/listinfo/r-devel
>      >
>      > ______________________________________________
>      > [hidden email] <mailto:[hidden email]> mailing list
>      > https://stat.ethz.ch/mailman/listinfo/r-devel
>      >
>
>     ______________________________________________
>     [hidden email] <mailto:[hidden email]> mailing list
>     https://stat.ethz.ch/mailman/listinfo/r-devel
>
> --
> Sent from Gmail Mobile

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Calling a LAPACK subroutine from R

Avraham Adler
I think better stated that it is subset that relied on a “bug” that was
never trapped for until now. We’re it written “properly” this never would
have arisen. At least to the best of my understanding.

Avi

On Wed, Sep 11, 2019 at 4:34 PM Göran Broström <[hidden email]>
wrote:

>
>
> On 2019-09-11 22:16, Avraham Adler wrote:
> > Can you write a small C function that calls LAPACK call that fro your
> > Fortran code? Yes, an extra step but maybe less traumatic than rewriting
> > parts of LAPACK directly.
>
> Yes, I know how to do that, but I find it somewhat bizarre that it is
> impossible to call a Fortran subroutine from Fortran. And rewriting
> 'dgemv' was simple: Just change character to integer and 'N' to 1. And
> rename the subroutine. The hard (tedious) part was to include all the
> LAPACK authors in my DESCRIPTION file.
>
> My guess is that the root cause is that BLAS/LAPACK is written in
> FORTRAN 77, which is said to be a subset of the current Fortran version
> but obviously isn't.
>
> Thanks, Göran
>
> >
> > Avi
> >
> > On Wed, Sep 11, 2019 at 4:08 PM Göran Broström <[hidden email]
> > <mailto:[hidden email]>> 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 <[hidden email]
> >     <mailto:[hidden email]>> wrote:
> >      >>
> >      >> Sorry for cross-posting, but I realized my question might be
> >     more appropriate for r-devel...
> >      >>
> >      >> Thank you,
> >      >> Giovanni
> >      >>
> >      >> ________________________________________
> >      >> From: R-help <[hidden email]
> >     <mailto:[hidden email]>> on behalf of Giovanni Petris
> >     <[hidden email] <mailto:[hidden email]>>
> >      >> Sent: Tuesday, September 10, 2019 16:44
> >      >> To: [hidden email] <mailto:[hidden email]>
> >      >> 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
> >      >>
> >      >>
> >      >> ______________________________________________
> >      >> [hidden email] <mailto:[hidden email]> 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.
> >      >>
> >      >> ______________________________________________
> >      >> [hidden email] <mailto:[hidden email]> mailing
> list
> >      >> https://stat.ethz.ch/mailman/listinfo/r-devel
> >      >
> >      > ______________________________________________
> >      > [hidden email] <mailto:[hidden email]> mailing list
> >      > https://stat.ethz.ch/mailman/listinfo/r-devel
> >      >
> >
> >     ______________________________________________
> >     [hidden email] <mailto:[hidden email]> mailing list
> >     https://stat.ethz.ch/mailman/listinfo/r-devel
> >
> > --
> > Sent from Gmail Mobile
>
--
Sent from Gmail Mobile

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Calling a LAPACK subroutine from R

Berend Hasselman
In reply to this post by Göran Broström-3

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 <[hidden email]> 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 <[hidden email]> wrote:
>>>
>>> Sorry for cross-posting, but I realized my question might be more appropriate for r-devel...
>>>
>>> Thank you,
>>> Giovanni
>>>
>>> ________________________________________
>>> From: R-help <[hidden email]> on behalf of Giovanni Petris <[hidden email]>
>>> Sent: Tuesday, September 10, 2019 16:44
>>> To: [hidden email]
>>> 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
>>>
>>>
>>> ______________________________________________
>>> [hidden email] 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.
>>>
>>> ______________________________________________
>>> [hidden email] mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Calling a LAPACK subroutine from R

Berend Hasselman
Followup:

I have checked my package nleqslv which uses dgemv only from Fortran, on Kubuntu 18.04 with the development version of R.
No errors or problems.

Berend


> On 12 Sep 2019, at 08:57, Berend Hasselman <[hidden email]> wrote:
>
>
> 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 <[hidden email]> 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 <[hidden email]> wrote:
>>>>
>>>> Sorry for cross-posting, but I realized my question might be more appropriate for r-devel...
>>>>
>>>> Thank you,
>>>> Giovanni
>>>>
>>>> ________________________________________
>>>> From: R-help <[hidden email]> on behalf of Giovanni Petris <[hidden email]>
>>>> Sent: Tuesday, September 10, 2019 16:44
>>>> To: [hidden email]
>>>> 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
>>>>
>>>>
>>>> ______________________________________________
>>>> [hidden email] 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.
>>>>
>>>> ______________________________________________
>>>> [hidden email] mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>> ______________________________________________
>>> [hidden email] mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Calling a LAPACK subroutine from R

Serguei Sokol
In reply to this post by Berend Hasselman
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 <[hidden email]> wrote:
>>
>> Sorry for cross-posting, but I realized my question might be more appropriate for r-devel...
>>
>> Thank you,
>> Giovanni
>>
>> ________________________________________
>> From: R-help <[hidden email]> on behalf of Giovanni Petris <[hidden email]>
>> Sent: Tuesday, September 10, 2019 16:44
>> To: [hidden email]
>> 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
>>
>>
>> ______________________________________________
>> [hidden email] 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.
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Calling a LAPACK subroutine from R

Göran Broström-3
In reply to this post by Berend Hasselman
Thanks Berend, I got curious and checked your package, and no errors.
However, two open questions:

1. You invoke the call to dgemv from a Fortran subroutine  that is
called from a C function, which in turn is called in R by .Call. I go
directly via .Fortran, no C involved, except by "base R", see 2. below.

2. I am still wondering why your route gets away with avoiding "the 12th
hidden argument" FCLEN, see R-devel/include/R_ext/BLAS.h:

F77_NAME(dgemv)(const char *trans, const int *m, const int *n,
                const double *alpha, const double *a, const int *lda,
                const double *x, const int *incx, const double *beta,
                double *y, const int *incy FCLEN);

when my direct route doesn't?

Göran

On 2019-09-12 09:15, Berend Hasselman wrote:

> Followup:
>
> I have checked my package nleqslv which uses dgemv only from Fortran, on Kubuntu 18.04 with the development version of R.
> No errors or problems.
>
> Berend
>
>
>> On 12 Sep 2019, at 08:57, Berend Hasselman <[hidden email]> wrote:
>>
>>
>> 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 <[hidden email]> 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 <[hidden email]> wrote:
>>>>>
>>>>> Sorry for cross-posting, but I realized my question might be more appropriate for r-devel...
>>>>>
>>>>> Thank you,
>>>>> Giovanni
>>>>>
>>>>> ________________________________________
>>>>> From: R-help <[hidden email]> on behalf of Giovanni Petris <[hidden email]>
>>>>> Sent: Tuesday, September 10, 2019 16:44
>>>>> To: [hidden email]
>>>>> 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
>>>>>
>>>>>
>>>>> ______________________________________________
>>>>> [hidden email] 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.
>>>>>
>>>>> ______________________________________________
>>>>> [hidden email] mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>> ______________________________________________
>>>> [hidden email] mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Calling a LAPACK subroutine from R

Berend Hasselman
In reply to this post by Serguei Sokol


> On 12 Sep 2019, at 10:36, Serguei Sokol <[hidden email]> wrote:
>
> 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)


This will not work on macOS.
The extension for dynamic libraries is .dylib.
So you would need

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

See the help for .Platform and dyn.load for the details for macOS.

Berend

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Calling a LAPACK subroutine from R

Göran Broström-3
In reply to this post by Göran Broström-3
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@R: c(person("Göran", "Broström", role = c("aut", "cre"),
>>>                      email = "[hidden email]"),
>>>               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 <[hidden email]> wrote:
>>>>>>
>>>>>> Sorry for cross-posting, but I realized my question might be more
>>>>>> appropriate for r-devel...
>>>>>>
>>>>>> Thank you,
>>>>>> Giovanni
>>>>>>
>>>>>> ________________________________________
>>>>>> From: R-help <[hidden email]> on behalf of Giovanni
>>>>>> Petris <[hidden email]>
>>>>>> Sent: Tuesday, September 10, 2019 16:44
>>>>>> To: [hidden email]
>>>>>> 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
>>>>>>
>>>>>>
>>>>>> ______________________________________________
>>>>>> [hidden email] 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.
>>>>>>
>>>>>> ______________________________________________
>>>>>> [hidden email] mailing list
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>>>
>>>>> ______________________________________________
>>>>> [hidden email] mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>>>
>>>
>>>> ______________________________________________
>>>> [hidden email] mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Calling a LAPACK subroutine from R

Serguei Sokol
In reply to this post by Berend Hasselman
On 12/09/2019 11:07, Berend Hasselman wrote:

>> On 12 Sep 2019, at 10:36, Serguei Sokol<[hidden email]>  wrote:
>>
>> 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)
> This will not work on macOS.
> The extension for dynamic libraries is .dylib.
> So you would need
>
> lapack.path <- paste0(file.path(R.home(), ifelse(.Platform$OS.type == "windows",
>           file.path("bin", .Platform$r_arch, "Rlapack"), file.path("lib", "libRlapack"))),
>           ".dylib")
>
> See the help for .Platform and dyn.load for the details for macOS.
Indeed. I was surprised to discover that .Platform$dynlib.ext is set to
".so" on macos,
not to ".dylib". Thank you to point me to a special note about it in
?.Platform
Is there a R predefined variable set to ".dylib" on macos ?
Meanwhile, the code for lapack path detection will become a little bit
more complicated:

dynlib.ext=ifelse(Sys.info()[["sysname"]] == "Darwin", ".dylib",
.Platform$dynlib.ext)
lapack.path <- paste0(file.path(R.home(), ifelse(.Platform$OS.type ==
"windows",
          file.path("bin", .Platform$r_arch, "Rlapack"),
file.path("lib", "libRlapack"))),
          dynlib.ext)
dyn.load(lapack.path)
is.loaded("dgemv") # must be TRUE

Serguei.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel