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

    ______________________________________________
    R-devel@r-project.org <mailto:R-devel@r-project.org> mailing list
    https://stat.ethz.ch/mailman/listinfo/r-devel

--
Sent from Gmail Mobile

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Reply via email to