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 = "goran.brost...@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 <gpet...@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-boun...@r-project.org> on behalf of Giovanni
Petris <gpet...@uark.edu>
Sent: Tuesday, September 10, 2019 16:44
To: r-h...@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-h...@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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
