From gor@n@bro@trom @end|ng |rom umu@@e  Thu Sep 12 10:41:01 2019
From: gor@n@bro@trom @end|ng |rom umu@@e (=?UTF-8?Q?G=c3=b6ran_Brostr=c3=b6m?=)
Date: Thu, 12 Sep 2019 10:41:01 +0200
Subject: [Rd] Calling a LAPACK subroutine from R
In-Reply-To: <66C8836B-379D-42FE-B6CF-CBD4D2336EF9@xs4all.nl>
References: <DM6PR04MB39154A2AC724DC1A1ED7ABFEC3B60@DM6PR04MB3915.namprd04.prod.outlook.com>
 <DM6PR04MB3915AAF7C4F24F2D854AD8FBC3B10@DM6PR04MB3915.namprd04.prod.outlook.com>
 <329F5311-5041-4297-B5C3-C05C258B1AF7@xs4all.nl>
 <ae8c6025-0a61-331e-c3e6-ce5fa1391da1@umu.se>
 <129C0E17-CBFC-4A53-A001-B768587B2BCC@xs4all.nl>
 <66C8836B-379D-42FE-B6CF-CBD4D2336EF9@xs4all.nl>
Message-ID: <da5a9a32-1860-c323-182c-071c5d803ba0@umu.se>

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 <bhh at xs4all.nl> 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 <goran.brostrom at umu.se> wrote:
>>>
>>> Berend,
>>>
>>> I do not think this works with gfortran 7+. I am calling the BLAS subroutine dgemv from Fortran code in my package eha, and the check (with R-devel) gives:
>>>
>>> gmlfun.f:223:1: warning: type of ?dgemv? does not match original declaration [-Wlto-type-mismatch]
>>>      &     score, ione)
>>> ^
>>> /home/gobr0002/R/src/R-devel/include/R_ext/BLAS.h:107:1: note: type mismatch in parameter 12
>>> F77_NAME(dgemv)(const char *trans, const int *m, const int *n,
>>>
>>> Type of a Fortran subroutine is matched against type of a C function?! My conclusion is that it is impossible to call a BLAS subroutine with a character parameter from Fortran code (nowadays). Calling from C code is fine, on the other hand(!).
>>>
>>> I have recently asked about this on R-pkg-devel, but not received any useful answers, and my submission to CRAN is rejected. I solve it by making a personal copy of dgemv and changing the character parameter to integer, and adding Jack Dongarra, Jeremy Du Croz, Sven Hammarling, and Richard Hanson as authors of eha. And a Copyright note, all in the DESCRIPTION file. Ugly but what can I do (except rewriting the Fortran code in C with f2c)?
>>>
>>> G?ran
>>>
>>> On 2019-09-11 21:38, Berend Hasselman wrote:
>>>> The Lapack library is loaded automatically by R itself when it needs it  for doing some calculation.
>>>> You can force it to do that with a (dummy) solve for example.
>>>> Put this at start of your script:
>>>> <code>
>>>> # dummy code to get LAPACK library loaded
>>>> X1 <- diag(2,2)
>>>> x1 <- rep(2,2)
>>>> # X1;x1
>>>> z <- solve(X1,x1)
>>>> </code>
>>>> followed by the rest of your script.
>>>> You will get a warning (I do) that  "passing a character vector  to .Fortran is not portable".
>>>> On other systems this may gave fatal errors. This is quick and very dirty. Don't do it.
>>>> I believe there is a better and much safer way to achieve what you want.
>>>> Here goes.
>>>> Create a folder (directory) src in the directory where your script resides.
>>>> Create a wrapper for "dpbtrf" file in a file xdpbtrf.f that takes an integer instead of character
>>>> <xdpbtrf.f>
>>>> c intermediate for dpbtrf
>>>>       SUBROUTINE xDPBTRF( kUPLO, N, KD, AB, LDAB, INFO )
>>>> c      .. Scalar Arguments ..
>>>>       integer         kUPLO
>>>>       INTEGER         INFO, KD, LDAB, N
>>>> c  .. Array Arguments ..
>>>>       DOUBLE PRECISION   AB( LDAB, * )
>>>>       character UPLO
>>>> c     convert integer argument to character
>>>>       if(kUPLO .eq. 1 ) then
>>>>           UPLO = 'L'
>>>>       else
>>>>           UPLO = 'U'
>>>>       endif
>>>>       call dpbtrf(UPLO,N,KD,AB,LDAB,INFO)
>>>>       return
>>>>       end
>>>> </xdpbtrf.f>
>>>> Instead of a character argument UPLO it takes an integer argument kUPLO.
>>>> The meaning should be obvious from the code.
>>>> Now create a shell script in the folder of your script to generate a dynamic library to be loaded in your script:
>>>> <mkso.sh>
>>>> # Build a binary dynamic library for accessing Lapack dpbtrf
>>>> # syntax checking
>>>> SONAME=xdpbtrf.so
>>>> echo Strict syntax checking
>>>> echo ----------------------
>>>> gfortran -c -fsyntax-only -fimplicit-none -Wall src/*.f || exit 1
>>>> LAPACK=$(R CMD config LAPACK_LIBS)
>>>> R CMD SHLIB --output=${SONAME} src/*.f ${LAPACK} || exit 1
>>>> </mkso.sh>
>>>> To load the dynamic library xdpbtrf.so  change your script into this
>>>> <yourscript>
>>>> dyn.load("xdpbtrf.so")
>>>> n <- 4L
>>>> phi <- 0.64
>>>> AB <- matrix(0, 2, n)
>>>> AB[1, ] <- c(1, rep(1 + phi^2, n-2), 1)
>>>> AB[2, -n] <- -phi
>>>> round(AB, 3)
>>>> AB.ch <- .Fortran("xdpbtrf", kUPLO=1L, N = as.integer(n),
>>>>                             KD = 1L, AB = AB, LDAB = 2L, INFO = as.integer(0))$AB
>>>> AB.ch
>>>> </yourscript>
>>>> and you are good to go.
>>>> You should always do something  as described above when you need to pass character arguments to Fortran code.
>>>> All of this was tested and run on macOS using the CRAN version of R.
>>>> Berend Hasselman
>>>>> On 11 Sep 2019, at 15:47, Giovanni Petris <gpetris at 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 at r-project.org> on behalf of Giovanni Petris <gpetris at uark.edu>
>>>>> Sent: Tuesday, September 10, 2019 16:44
>>>>> To: r-help at 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 at 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 at r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>> ______________________________________________
>>>> R-devel at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>



[Rd] Calling a LAPACK subroutine from R
=?UTF-8?Q?G=c3=b6ran_Brostr=c3=b6m?=