From gor@n@bro@trom @end|ng |rom umu@@e  Thu Sep 12 11:56:05 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 11:56:05 +0200
Subject: [Rd] Fw: Calling a LAPACK subroutine from R
In-Reply-To: <95db5dc3-19d4-4eca-fb32-53d4bef59f3a@gmail.com>
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>
 <23929.59590.434134.81889@hornik.net>
 <517a2e8d-10d7-ad58-d547-07b4b7411570@umu.se>
 <95db5dc3-19d4-4eca-fb32-53d4bef59f3a@gmail.com>
Message-ID: <3d293c58-bc55-7c84-4ed5-b6805dd40eb0@umu.se>

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 at R: c(person("G?ran", "Brostr?m", role = c("aut", "cre"),
>>> ???????????????????? email = "goran.brostrom at umu.se"),
>>> ????????????? person("Jack", "Dongarra", role = "aut"),
>>> ????????????? person("Jeremy", "Du Croz", role = "aut"),
>>> ????????????? person("Sven", "Hammarling", role = "aut"),
>>> ????????????? person("Richard", "Hanson", role = "aut"),
>>> ????????????? person("Jianming", "Jin", role = "aut"))
>>> Copyright: The blas Developers 2014-2018.
>>>
>>> I think the last 5 persons should get a ctb role instead of aut, and you
>>> should add
>>>
>>> ?? person("The blas Developers", role = "cph", comment = "....")
>>>
>>> with a suitable comment explaining what the copyright is held for.
>>>
>>> However, I would still hope that you can do without the copy by
>>> following the recipe in section "Fortran character strings" in "Writing
>>> R Extensions".? This says:
>>>
>>> ************************************************************************
>>> Alternatively, do as R does as from version 3.6.1-patched and pass the
>>> character length(s) from C to Fortran.? A portable way to do this is
>>> ????? // before any R headers, or define in PKG_CPPFLAGS
>>> ????? #define USE_FC_LEN_T
>>> ????? #include <Rconfig.h>
>>> ????? #include <R_ext/BLAS.h>
>>> ????? #ifndef FCONE
>>> ????? # define FCONE
>>> ????? #endif
>>> ????? ...
>>> ????????????? F77_CALL(dgemm)("N", "T", &nrx, &ncy, &ncx, &one, x,
>>> ????????????????????????????? &nrx, y, &nry, &zero, z, &nrx FCONE 
>>> FCONE);
>>> (Note there is no comma before or between the 'FCONE' invocations.)? It
>>> is strongly recommended that packages which call from C/C++ BLAS/LAPACK
>>> routines with character arguments adopt this approach.
>>> ************************************************************************
>>>
>>> Does this really not work for you?
>>
>> Of course it would work: I call dgemv (and other BLAS routines) 
>> frequently from C code in my package, that is not the problem. I am 
>> just wondering why Berend gets away with direct calls to dgemv from 
>> Fortran without using the 12th hidden parameter FCLEN (FCONE), when I 
>> fail. I suspect an oversight in the R checking system.
> 
> The hidden argument will be added by the Fortran compiler, at the call 
> site. If BLAS is built with the same Fortran compiler as your Fortran 
> code which calls it, there should be no problem (often it will work also 
> when the compilers are different, but that depends on how different they 
> are). If you can't find why it is not working for you, please create a 
> minimum reproducible but complete example. Please try looking also at 
> existing packages with Fortran code, including "stats" for instance. 
> There are several people who are trying to help you on the mailing 
> lists, it would be easier for them if they knew exactly what you were 
> doing.
> 
> Best
> Tomas
> 
> 
>>
>> I can write a C wrapper for my Fortran function that calls dgemv and 
>> call it in R by .C, but is that really guaranteed to work, when the 
>> checking routine gets sharper?
> 
>>
>> Best, G,
>>
>>> Best
>>> -k
>>>
>>>
>>>> Berend,
>>>> I do not think this works with gfortran 7+. I am calling the BLAS
>>>> subroutine dgemv from Fortran code in my package eha, and the check
>>>> (with R-devel) gives:
>>>
>>>> gmlfun.f:223:1: warning: type of ?dgemv? does not match original
>>>> declaration [-Wlto-type-mismatch]
>>>> ??????? &???? score, ione)
>>>> ?? ^
>>>> /home/gobr0002/R/src/R-devel/include/R_ext/BLAS.h:107:1: note: type
>>>> mismatch in parameter 12
>>>> ?? F77_NAME(dgemv)(const char *trans, const int *m, const int *n,
>>>
>>>> Type of a Fortran subroutine is matched against type of a C function?!
>>>> My conclusion is that it is impossible to call a BLAS subroutine with a
>>>> character parameter from Fortran code (nowadays). Calling from C 
>>>> code is
>>>> fine, on the other hand(!).
>>>
>>>> I have recently asked about this on R-pkg-devel, but not received any
>>>> useful answers, and my submission to CRAN is rejected. I solve it by
>>>> making a personal copy of dgemv and changing the character parameter to
>>>> integer, and adding Jack Dongarra, Jeremy Du Croz, Sven Hammarling, and
>>>> Richard Hanson as authors of eha. And a Copyright note, all in the
>>>> DESCRIPTION file. Ugly but what can I do (except rewriting the Fortran
>>>> code in C with f2c)?
>>>
>>>> G?ran
>>>
>>>> On 2019-09-11 21:38, Berend Hasselman wrote:
>>>>>
>>>>> The Lapack library is loaded automatically by R itself when it 
>>>>> needs it? for doing some calculation.
>>>>> You can force it to do that with a (dummy) solve for example.
>>>>> Put this at start of your script:
>>>>>
>>>>> <code>
>>>>> # dummy code to get LAPACK library loaded
>>>>> X1 <- diag(2,2)
>>>>> x1 <- rep(2,2)
>>>>> # X1;x1
>>>>> z <- solve(X1,x1)
>>>>> </code>
>>>>>
>>>>> followed by the rest of your script.
>>>>> You will get a warning (I do) that? "passing a character vector? to 
>>>>> .Fortran is not portable".
>>>>> On other systems this may gave fatal errors. This is quick and very 
>>>>> dirty. Don't do it.
>>>>>
>>>>> I believe there is a better and much safer way to achieve what you 
>>>>> want.
>>>>> Here goes.
>>>>>
>>>>> Create a folder (directory) src in the directory where your script 
>>>>> resides.
>>>>> Create a wrapper for "dpbtrf" file in a file xdpbtrf.f that takes 
>>>>> an integer instead of character
>>>>>
>>>>> <xdpbtrf.f>
>>>>> c intermediate for dpbtrf
>>>>>
>>>>> SUBROUTINE xDPBTRF( kUPLO, N, KD, AB, LDAB, INFO )
>>>>>
>>>>> c????? .. Scalar Arguments ..
>>>>> integer???????? kUPLO
>>>>> INTEGER???????? INFO, KD, LDAB, N
>>>>>
>>>>> c? .. Array Arguments ..
>>>>> DOUBLE PRECISION?? AB( LDAB, * )
>>>>>
>>>>> character UPLO
>>>>> c???? convert integer argument to character
>>>>> if(kUPLO .eq. 1 ) then
>>>>> UPLO = 'L'
>>>>> else
>>>>> UPLO = 'U'
>>>>> endif
>>>>>
>>>>> call dpbtrf(UPLO,N,KD,AB,LDAB,INFO)
>>>>> return
>>>>> end
>>>>> </xdpbtrf.f>
>>>>>
>>>>>
>>>>> Instead of a character argument UPLO it takes an integer argument 
>>>>> kUPLO.
>>>>> The meaning should be obvious from the code.
>>>>>
>>>>> Now create a shell script in the folder of your script to generate 
>>>>> a dynamic library to be loaded in your script:
>>>>>
>>>>> <mkso.sh>
>>>>> # Build a binary dynamic library for accessing Lapack dpbtrf
>>>>>
>>>>> # syntax checking
>>>>>
>>>>> SONAME=xdpbtrf.so
>>>>>
>>>>> echo Strict syntax checking
>>>>> echo ----------------------
>>>>> gfortran -c -fsyntax-only -fimplicit-none -Wall src/*.f || exit 1
>>>>>
>>>>> LAPACK=$(R CMD config LAPACK_LIBS)
>>>>> R CMD SHLIB --output=${SONAME} src/*.f ${LAPACK} || exit 1
>>>>> </mkso.sh>
>>>>>
>>>>> To load the dynamic library xdpbtrf.so? change your script into this
>>>>>
>>>>> <yourscript>
>>>>> dyn.load("xdpbtrf.so")
>>>>> n <- 4L
>>>>> phi <- 0.64
>>>>> AB <- matrix(0, 2, n)
>>>>> AB[1, ] <- c(1, rep(1 + phi^2, n-2), 1)
>>>>> AB[2, -n] <- -phi
>>>>> round(AB, 3)
>>>>>
>>>>> AB.ch <- .Fortran("xdpbtrf", kUPLO=1L, N = as.integer(n),
>>>>> KD = 1L, AB = AB, LDAB = 2L, INFO = as.integer(0))$AB
>>>>> AB.ch
>>>>>
>>>>> </yourscript>
>>>>>
>>>>> and you are good to go.
>>>>>
>>>>> You should always do something? as described above when you need to 
>>>>> pass character arguments to Fortran code.
>>>>>
>>>>> All of this was tested and run on macOS using the CRAN version of R.
>>>>>
>>>>> Berend Hasselman
>>>>>
>>>>>> On 11 Sep 2019, at 15:47, Giovanni Petris <gpetris 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] Fw: Calling a LAPACK subroutine from R
=?UTF-8?Q?G=c3=b6ran_Brostr=c3=b6m?=