[Rd] [ALTREP] What is the meaning of the return value of Is_sorted and No_NA function?

[Wang Jiefei]

Hi,



I would like to figure out the meaning of the return value of these two
functions. Here are the default definitions I find from R source code:



static int altreal_Is_sorted_default(SEXP x) { return UNKNOWN_SORTEDNESS; }

static int altreal_No_NA_default(SEXP x) { return 0; }

I guess the macro *UNKNOWN_SORTEDNESS *in *Is_sorted* and 0 in *No_NA
*simply means
unknown sorted/NA status of the vector, so R will loop over the vector and
find the answer. However, what should we return in these functions to
indicate whether the vector has been sorted/ contains NA? My initial guess
is 0/1 but since *NA_NA *uses 0 as its default value so it will be
ambiguous. Are there any macros to define yes/no return values for these
functions? I would appreciate any thought here.



Best,

Jiefei

[[alternative HTML version deleted]]

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

Re: [Rd] '==' operator: inconsistency in data.frame(...) == NULL

[Martin Maechler]

> Hilmar Berger 
> on Wed, 4 Sep 2019 15:25:46 +0200 writes:

> Dear all,

> I just stumbled upon some behavior of the == operator which is at least 
> somewhat inconsistent.

> R version 3.6.1 (2019-07-05) -- "Action of the Toes"
> Copyright (C) 2019 The R Foundation for Statistical Computing
> Platform: x86_64-w64-mingw32/x64 (64-bit)

>> list(a=1:3, b=LETTERS[1:3]) == NULL
> logical(0)
>> matrix(1:6, 2,3) == NULL
> logical(0)
>> data.frame(a=1:3, b=LETTERS[1:3]) == NULL # same for == logical(0)
> Error in matrix(if (is.null(value)) logical() else value, nrow = nr, 
> dimnames = list(rn,  :
>   length of 'dimnames' [2] not equal to array extent

>> data.frame(NULL) == 1
> <0 x 0 matrix>
>> data.frame(NULL) == NULL
> <0 x 0 matrix>
>> data.frame(NULL) == logical(0)
> <0 x 0 matrix>

> I wonder if data.frame() == NULL should also return 
> a value instead of an error. R help reads:

> "At least one of |x| and |y| must be an atomic vector, but
>  if the other is a list R attempts to coerce it to the
>  type of the atomic vector: this will succeed if the list
>  is made up of elements of length one that can be coerced
>  to the correct type.

>  If the two arguments are atomic vectors of different
>  types, one is coerced to the type of the other, the
>  (decreasing) order of precedence being character, complex,
>  numeric, integer, logical and raw."

> It is not clear from the help what to expect for NULL or
> empty atomic vectors. 

Well, strictly speaking an error would be expected for NULL,
as it is *not* an atomic vector, and your main issue

 " data.frame(..) == NULL "

would already be settled by the first half sentence from the
doc, and strictly speaking, even  data.frame(NULL) == NULL
"should" return an error ((Note: I'm not saying it really
 should, but at least the reference does not say it should work at all))

Now,  logical(0)  on the other hand *is* an atomic vector ... 


> It is also strange that for list()
> there is no error but for data.frame() with the same data
> an error is thrown. I can see that there might be reasons
> to return logical(0) instead of FALSE, but I do not fully
> understand why there should be differences between
> e.g. matrix() and data.frame().

Well, a [regular base R] matrix() is atomic  and a data frame is not.

> Also, It is at least somewhat strange that
> data.frame(NULL) == NULL and similar expressions return an
> empty matrix, while comparing a normal filled matrix to
> NULL returns logical(0).

> Even if this behavior is expected, the error message shown
> by data.frame(...) == NULL is not very informative.

I'm not at all sure there's any need for a change here.

I would say the following general thinking should be applied

1. The general rule that '==' should be used only for comparing 
  atomic objects (as it returns an atomic object, a 'logical' with
  corresponding attributes), is really principal
  and using '==' for anything else has never been "the idea".

2. There are (two) "semi-exceptions" to the above:
2a) Sometimes it has been convenient to treat NULL as if it was
 a zero-length atomic object (of "arbitrary" type/mode).
2b) data.frame()s "should typically" behave like matrices in
many situations, notably when indexed {and that rule is
violated (on purpose) by tibbles .. ("drop=FALSE" etc, but
that's another story)} 

So because of these exceptions, you and possibly others may
think  '=='  should "work" with data.frame()s and/or NULL, but
I would not tend to agree.

> Thanks and best regards,
> Hilmar

You are welcome!
Martin

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

Re: [Rd] '==' operator: inconsistency in data.frame(...) == NULL

[Hilmar Berger]

Dear Martin,

On 11/09/2019 09:56, Martin Maechler wrote:
>
>  > I wonder if data.frame() == NULL should also 
> return
>  > a value instead of an error. R help reads:
>
>  > "At least one of |x| and |y| must be an atomic vector, but
>   >  if the other is a list R attempts to coerce it to the
>   >  type of the atomic vector: this will succeed if the list
>   >  is made up of elements of length one that can be coerced
>   >  to the correct type.
>
>   >  If the two arguments are atomic vectors of different
>   >  types, one is coerced to the type of the other, the
>   >  (decreasing) order of precedence being character, complex,
>   >  numeric, integer, logical and raw."
>
>  > It is not clear from the help what to expect for NULL or
>  > empty atomic vectors.
>
> Well, strictly speaking an error would be expected for NULL,
> as it is *not* an atomic vector, and your main issue
>
>   " data.frame(..) == NULL "
>
> would already be settled by the first half sentence from the
> doc, and strictly speaking, even  data.frame(NULL) == NULL
> "should" return an error ((Note: I'm not saying it really
>   should, but at least the reference does not say it should work at all))
Thanks, this explanation makes total sense to me. I did not consider 
that NULL might be non-atomic. Strangely, is.atomic(NULL) returns TRUE. 
On the other hand, I understand that one would not like to treat it like 
atomic in ==.

However, in this case one might expect that the error message would be 
more like that for S4 objects (which always seem to report an 
informative error message for ==):

 > Pos <- setClass("Pos", slots = c(latitude = "numeric", longitude = 
"numeric", altitude = "numeric"))
 > p = Pos()
 > p == NULL
Error in p == NULL :
   comparison (1) is possible only for atomic and list types
 > p == "FOO"
Error in p == "FOO" :
   comparison (1) is possible only for atomic and list types

In the data.frame()==NULL cases I have the impression that the fact that 
both sides are non-atomic is not properly detected and therefore R tries 
to go on with the == method for data.frames.

 From a cursory check in Ops.data.frame() and some debugging I have the 
impression that the case of the second argument being non-atomic or 
empty is not handled at all and the function progresses until the end, 
where it fails in the last step on an empty value:

matrix(unlist(value, recursive = FALSE, use.names = FALSE),
     nrow = nr, dimnames = list(rn, cn))

Best regards,
Hilmar

-- 
Dr. Hilmar Berger, MD
Max Planck Institute for Infection Biology
Charitéplatz 1
D-10117 Berlin
GERMANY

Phone:  + 49 30 28460 430
Fax:+ 49 30 28460 401
  
E-Mail: ber...@mpiib-berlin.mpg.de
Web   : www.mpiib-berlin.mpg.de


[[alternative HTML version deleted]]

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

Re: [Rd] '==' operator: inconsistency in data.frame(...) == NULL

[Hilmar Berger]

Another example where a data.frame is compared to (here non-null, 
non-empty) non-atomic values in Ops.data.frame, resulting in an error 
message:


setClass("FOOCLASS2",
 slots = c(M="matrix")
)
ma = new("FOOCLASS2", M=matrix(rnorm(300), 30,10))

> isS4(ma)
[1] TRUE
> ma == data.frame(a=1:3)
Error in eval(f) : dims [product 1] do not match the length of object [3]

As for the NULL/logical(0) cases I would suggest to explicitly test for 
invalid conditions in Ops.data.frame and generate a comprehensible 
message (e.g. "comparison is possible only for atomic and list types") 
if appropriate.


Best regards,
Hilmar


On 11/09/2019 11:55, Hilmar Berger wrote:


In the data.frame()==NULL cases I have the impression that the fact 
that both sides are non-atomic is not properly detected and therefore 
R tries to go on with the == method for data.frames.


From a cursory check in Ops.data.frame() and some debugging I have the 
impression that the case of the second argument being non-atomic or 
empty is not handled at all and the function progresses until the end, 
where it fails in the last step on an empty value:


matrix(unlist(value, recursive = FALSE, use.names = FALSE),
    nrow = nr, dimnames = list(rn, cn)) 


--
Dr. Hilmar Berger, MD
Max Planck Institute for Infection Biology
Charitéplatz 1
D-10117 Berlin
GERMANY

Phone:  + 49 30 28460 430
Fax:+ 49 30 28460 401
 
E-Mail: ber...@mpiib-berlin.mpg.de

Web   : www.mpiib-berlin.mpg.de

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

Re: [Rd] '==' operator: inconsistency in data.frame(...) == NULL

[Hilmar Berger]

Sorry, I can't reproduce the example below even on the same machine. 
However, the following example produces the same error as NULL values in 
prior examples:


> setClass("FOOCLASS",
+  representation("list")
+ )
> ma = new("FOOCLASS", list(M=matrix(rnorm(300), 30,10)))
> isS4(ma)
[1] TRUE
> data.frame(a=1:3) == ma
Error in matrix(unlist(value, recursive = FALSE, use.names = FALSE), 
nrow = nr,  :

  length of 'dimnames' [2] not equal to array extent

Best,
Hilmar


On 11/09/2019 12:24, Hilmar Berger wrote:
Another example where a data.frame is compared to (here non-null, 
non-empty) non-atomic values in Ops.data.frame, resulting in an error 
message:


setClass("FOOCLASS2",
 slots = c(M="matrix")
)
ma = new("FOOCLASS2", M=matrix(rnorm(300), 30,10))

> isS4(ma)
[1] TRUE
> ma == data.frame(a=1:3)
Error in eval(f) : dims [product 1] do not match the length of object [3]

As for the NULL/logical(0) cases I would suggest to explicitly test 
for invalid conditions in Ops.data.frame and generate a comprehensible 
message (e.g. "comparison is possible only for atomic and list types") 
if appropriate.


Best regards,
Hilmar


On 11/09/2019 11:55, Hilmar Berger wrote:


In the data.frame()==NULL cases I have the impression that the fact 
that both sides are non-atomic is not properly detected and therefore 
R tries to go on with the == method for data.frames.


From a cursory check in Ops.data.frame() and some debugging I have 
the impression that the case of the second argument being non-atomic 
or empty is not handled at all and the function progresses until the 
end, where it fails in the last step on an empty value:


matrix(unlist(value, recursive = FALSE, use.names = FALSE),
    nrow = nr, dimnames = list(rn, cn)) 




--
Dr. Hilmar Berger, MD
Max Planck Institute for Infection Biology
Charitéplatz 1
D-10117 Berlin
GERMANY

Phone:  + 49 30 28460 430
Fax:+ 49 30 28460 401
 
E-Mail: ber...@mpiib-berlin.mpg.de

Web   : www.mpiib-berlin.mpg.de

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

[Rd] Fw: Calling a LAPACK subroutine from R

[Giovanni Petris]

Sorry for cross-posting, but I realized my question might be more appropriate 
for r-devel...

Thank you,
Giovanni


From: R-help  on behalf of Giovanni Petris 

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.411
[2,] -0.64 -0.64 -0.640
>
> ### 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

Re: [Rd] [ALTREP] What is the meaning of the return value of Is_sorted and No_NA function?

[Gabriel Becker]

Hi Jiefei,

The meanings of the return values for sortedness can be found in
RInternals.h, and are as follows:

/* ALTREP sorting support */
enum {SORTED_DECR_NA_1ST = -2,
  SORTED_DECR = -1,
  UNKNOWN_SORTEDNESS = INT_MIN, /*INT_MIN is NA_INTEGER! */
  SORTED_INCR = 1,
  SORTED_INCR_NA_1ST = 2,
  KNOWN_UNSORTED = 0};

The default value there is NA_INTEGER (ie INT_MIN), indicating that there
is no sortedness information.

Currently, *_NO_NA  effectively return a boolean, (even though the actual
return value is int). This can be seen in the method we provide for compact
sequences in altclasses.c:


static int compact_intseq_No_NA(SEXP x)
{
#ifdef COMPACT_INTSEQ_MUTABLE
/* If the vector has been expanded it may have been modified. */
if (COMPACT_SEQ_EXPANDED(x) != R_NilValue)
return FALSE;
#endif
return TRUE;
}

(FALSE is a macro for 0, TRUE is a macro for 1).

Think of the meaning of the return value to No_NA methods as the object's
answer to the following question

"Are you sure there are zero NAs in your data?"

When it is sure of that, it  says "yes" (returning 1, ie TRUE). When it
either is sure there are NAs *OR* doesn't have any information about
whether there are NAs, it says "no" (returning 0, ie FALSE).

Also please note, it is possible there may be another API point in the
future which asks the object *how many NAs it has.∫ˆ* If that materializes,
No_NA would just  consume the answer to thatto get the binarized version,
but again there is nothing like that in there now.

Hope that helps.

Best,
~G

On Wed, Sep 11, 2019 at 12:04 AM Wang Jiefei  wrote:

> Hi,
>
>
>
> I would like to figure out the meaning of the return value of these two
> functions. Here are the default definitions I find from R source code:
>
>
>
> static int altreal_Is_sorted_default(SEXP x) { return UNKNOWN_SORTEDNESS; }
>
> static int altreal_No_NA_default(SEXP x) { return 0; }
>
> I guess the macro *UNKNOWN_SORTEDNESS *in *Is_sorted* and 0 in *No_NA
> *simply means
> unknown sorted/NA status of the vector, so R will loop over the vector and
> find the answer. However, what should we return in these functions to
> indicate whether the vector has been sorted/ contains NA? My initial guess
> is 0/1 but since *NA_NA *uses 0 as its default value so it will be
> ambiguous. Are there any macros to define yes/no return values for these
> functions? I would appreciate any thought here.
>
>
>
> Best,
>
> Jiefei
>
> [[alternative HTML version deleted]]
>
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

[[alternative HTML version deleted]]

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

Re: [Rd] [ALTREP] What is the meaning of the return value of Is_sorted and No_NA function?

[Wang Jiefei]

Hi Gabriel,

Thanks for your answer and future update plan. Somehow this email has been
delayed for a week, so there might be a wired reply from me saying that I
have found the answer from the R source code, it was sent from me last
week. Hopefully, this reply will not cost another week to post:)

As a side note, I like the idea that defining a macro for sortedness, and I
can see why we can only have a binary answer for NO_NA (since the return
value is actually bool). For making the code more readable, and for
possibly working with the future R release, is it possible to define a
macro for NO_NA function in RInternal.h? So if there is any change in NO_NA
function, there is no need to modify the code. Also, the code can be more
readable by doing that.

Best,
Jiefei

On Wed, Sep 11, 2019 at 1:58 PM Gabriel Becker 
wrote:

> Hi Jiefei,
>
> The meanings of the return values for sortedness can be found in
> RInternals.h, and are as follows:
>
> /* ALTREP sorting support */
> enum {SORTED_DECR_NA_1ST = -2,
>   SORTED_DECR = -1,
>   UNKNOWN_SORTEDNESS = INT_MIN, /*INT_MIN is NA_INTEGER! */
>   SORTED_INCR = 1,
>   SORTED_INCR_NA_1ST = 2,
>   KNOWN_UNSORTED = 0};
>
> The default value there is NA_INTEGER (ie INT_MIN), indicating that there
> is no sortedness information.
>
> Currently, *_NO_NA  effectively return a boolean, (even though the actual
> return value is int). This can be seen in the method we provide for compact
> sequences in altclasses.c:
>
>
> static int compact_intseq_No_NA(SEXP x)
> {
> #ifdef COMPACT_INTSEQ_MUTABLE
> /* If the vector has been expanded it may have been modified. */
> if (COMPACT_SEQ_EXPANDED(x) != R_NilValue)
> return FALSE;
> #endif
> return TRUE;
> }
>
> (FALSE is a macro for 0, TRUE is a macro for 1).
>
> Think of the meaning of the return value to No_NA methods as the object's
> answer to the following question
>
> "Are you sure there are zero NAs in your data?"
>
> When it is sure of that, it  says "yes" (returning 1, ie TRUE). When it
> either is sure there are NAs *OR* doesn't have any information about
> whether there are NAs, it says "no" (returning 0, ie FALSE).
>
> Also please note, it is possible there may be another API point in the
> future which asks the object *how many NAs it has.∫ˆ* If that
> materializes, No_NA would just  consume the answer to thatto get the
> binarized version, but again there is nothing like that in there now.
>
> Hope that helps.
>
> Best,
> ~G
>
> On Wed, Sep 11, 2019 at 12:04 AM Wang Jiefei  wrote:
>
>> Hi,
>>
>>
>>
>> I would like to figure out the meaning of the return value of these two
>> functions. Here are the default definitions I find from R source code:
>>
>>
>>
>> static int altreal_Is_sorted_default(SEXP x) { return UNKNOWN_SORTEDNESS;
>> }
>>
>> static int altreal_No_NA_default(SEXP x) { return 0; }
>>
>> I guess the macro *UNKNOWN_SORTEDNESS *in *Is_sorted* and 0 in *No_NA
>> *simply means
>> unknown sorted/NA status of the vector, so R will loop over the vector and
>> find the answer. However, what should we return in these functions to
>> indicate whether the vector has been sorted/ contains NA? My initial guess
>> is 0/1 but since *NA_NA *uses 0 as its default value so it will be
>> ambiguous. Are there any macros to define yes/no return values for these
>> functions? I would appreciate any thought here.
>>
>>
>>
>> Best,
>>
>> Jiefei
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>

[[alternative HTML version deleted]]

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

Re: [Rd] Fw: Calling a LAPACK subroutine from R

[Berend Hasselman]

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:


# dummy code to get LAPACK library loaded
X1 <- diag(2,2)
x1 <- rep(2,2)
# X1;x1
z <- solve(X1,x1)


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


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



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:


# 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


To load the dynamic library xdpbtrf.so  change your script into this


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



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  wrote:
> 
> Sorry for cross-posting, but I realized my question might be more appropriate 
> for r-devel...
> 
> Thank you,
> Giovanni
> 
> 
> From: R-help  on behalf of Giovanni Petris 
> 
> 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.411
> [2,] -0.64 -0.64 -0.640
>> 
>> ### 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 c

Re: [Rd] R-intro: Appendix A: attach position

[Kurt Hornik]

> Suharto Anggono Suharto Anggono via R-devel writes:

Thanks: fixed now in the trunk.

Best
-k

> In "An Introduction to R", in "Appendix A  A sample session", in the part on 
> Michelson data, information for
> attach(mm)
> is
> Make the data frame visible at position 3 (the default).

> In fact, the 'attach' is at position 2.

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

Re: [Rd] Fw: Calling a LAPACK subroutine from R

[Göran Broström]

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:


# dummy code to get LAPACK library loaded
X1 <- diag(2,2)
x1 <- rep(2,2)
# X1;x1
z <- solve(X1,x1)


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


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



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:


# 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


To load the dynamic library xdpbtrf.so  change your script into this


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



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

Sorry for cross-posting, but I realized my question might be more appropriate 
for r-devel...

Thank you,
Giovanni


From: R-help  on behalf of Giovanni Petris 

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.411
[2,] -0.64 -0.64 -0.640


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

Re: [Rd] Fw: Calling a LAPACK subroutine from R

[Avraham Adler]

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.

Avi

On Wed, Sep 11, 2019 at 4:08 PM Göran Broström 
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:
> >
> > 
> > # dummy code to get LAPACK library loaded
> > X1 <- diag(2,2)
> > x1 <- rep(2,2)
> > # X1;x1
> > z <- solve(X1,x1)
> > 
> >
> > 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
> >
> > 
> > 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
> > 
> >
> >
> > 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:
> >
> > 
> > # 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
> > 
> >
> > To load the dynamic library xdpbtrf.so  change your script into this
> >
> > 
> > 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
> >
> > 
> >
> > 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  wrote:
> >>
> >> Sorry for cross-posting, but I realized my question might be more
> appropriate for r-devel...
> >>
> >> Thank you,
> >> Giovanni
> >>
> >> 
> >> From: R-help  on behalf of Giovanni
> Petris 
> >> 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

Re: [Rd] Fw: Calling a LAPACK subroutine from R

[Göran Broström]

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 > 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:
 >
 > 
 > # dummy code to get LAPACK library loaded
 > X1 <- diag(2,2)
 > x1 <- rep(2,2)
 > # X1;x1
 > z <- solve(X1,x1)
 > 
 >
 > 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
 >
 > 
 > 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
 > 
 >
 >
 > 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:
 >
 > 
 > # 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
 > 
 >
 > To load the dynamic library xdpbtrf.so  change your script into this
 >
 > 
 > 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
 >
 > 
 >
 > and you are good to go.
 >
 > You should always do something  as descri

Re: [Rd] Fw: Calling a LAPACK subroutine from R

[Avraham Adler]

I think better stated that it is subset that relied on a “bug” that was
never trapped for until now. We’re it written “properly” this never would
have arisen. At least to the best of my understanding.

Avi

On Wed, Sep 11, 2019 at 4:34 PM Göran Broström 
wrote:

>
>
> 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  > > 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:
> >  >
> >  > 
> >  > # dummy code to get LAPACK library loaded
> >  > X1 <- diag(2,2)
> >  > x1 <- rep(2,2)
> >  > # X1;x1
> >  > z <- solve(X1,x1)
> >  > 
> >  >
> >  > 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
> >  >
> >  > 
> >  > 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
> >  > 
> >  >
> >  >
> >  > 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:
> >  >
> >  > 
> >  > # Build a binary dynamic library for accessing Lapack dpbtrf
> >  >
> >  > # syntax checking
> >  >
> >  > SONAME=xdpbtrf.so
> >  >
> >  > echo Strict syntax checking
> >  > echo --
> >  > gfortran -c -fsyntax-only -fimplicit-no

Re: [Rd] Calling a LAPACK subroutine from R

[Berend Hasselman]

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  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:
>> 
>> # dummy code to get LAPACK library loaded
>> X1 <- diag(2,2)
>> x1 <- rep(2,2)
>> # X1;x1
>> z <- solve(X1,x1)
>> 
>> 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
>> 
>> 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
>> 
>> 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:
>> 
>> # 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
>> 
>> To load the dynamic library xdpbtrf.so  change your script into this
>> 
>> 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
>> 
>> 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  wrote:
>>> 
>>> Sorry for cross-posting, but I realized my question might be more 
>>> appropriate for r-devel...
>>> 
>>> Thank you,
>>> Giovanni
>>> 
>>> 
>>> From: R-help  on behalf of Giovanni Petris 
>>> 
>>> 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.411
>>> [2,] -0.64 -0.64 -0.640
 
 ### Cholesk