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 <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
