Re: [Rd] F77_CALL(dgetrs) C++ call in R-devel
Thanks I got it to work!
I will try to migrate to (Rcpp)Armadillo asap.
Lars
Den tir. 20. dec. 2022 kl. 14.09 skrev Dirk Eddelbuettel :
>
> On 20 December 2022 at 12:33, Lars Relund wrote:
> | In my package, I have the method:
> |
> | /** Solve equations transpose(P)w = r. */
> | int LASolveT(MatSimple &P, MatSimple &w, const
> | MatSimple &r) {
> | int rows = P.rows;
> | int nrhs = 1;
> | int lda = rows;
> | int ldb = rows;
> | int info = -1;
> | MatSimple ipivot(1,rows);
> | w.Inject(r);// copy r to w;
> | F77_CALL(dgetrf)(&rows, &rows, &P(0,0), &lda, &ipivot(0,0),
> &info);
> | if (info!=0) {
> | cout << "Error in LASolve (dgetrf). Info=" << info << endl;
> | return 1;
> | }
> | F77_CALL(dgetrs)("T", &rows, &nrhs, &P(0,0), &lda, &ipivot(0,0),
> | &w(0,0), &ldb, &info);("T", &rows, &nrhs, &P(0,0), &lda, &ipivot(0,0),
> | &w(0,0), &ldb, &info);
> | if (info!=0) {
> | cout << "Error in LASolve (dgetrs). Info=" << info << endl;
> | return 1;
> | }
> | return 0;
> | }
> |
> | When compiling the package on using R-devel the error for
> F77_CALL(dgetrs)
> | occur:
> |
> | matalg.h:67:25: error: too few arguments to function ‘void dgetrs_(const
> | char*, const int*, const int*, const double*, const int*, const int*,
> | double*, const int*, int*, size_t)’
> |67 | F77_CALL(dgetrs)("T", &rows, &nrhs, &P(0,0), &lda,
> | &ipivot(0,0), &w(0,0), &ldb, &info);
> |
> | It works in R-release and I guess it have something to do with
> | https://cran.r-project.org/doc/manuals/r-devel/NEWS.html and LAPACK.
> |
> | Any hints on how to get it to work for both R-release and R-devel.
>
> I can offer you two answers. The first, and narrower, is in Writing R
> Extensions and concerns FC_LEN. R now 'automagically' injects additional
> parameters for a better, more stringent, control of character variable
> length. See eg
>
>
> https://rstudio.github.io/r-manuals/r-exts/The-R-API.html#fortran-character-strings
>
> and related. (And this isn't new per se, those of use with packages with
> Fortran interfaces have been keeping this up.)
>
> The second, more pragmatic answer, is of course 'to not do that' but to
> rely
> on the decade of tuning and bazillion of test and runs a higher-end Linear
> Algebra package like (Rcpp)Armadillo offers by wrapping around LAPACK and
> BLAS for you. You already are in C++, so there os essentially no switching
> cost. And (Rcpp)Armadillo is header-only and hence free of added
> complications.
>
> Dirk
>
> --
> dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org
>
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] F77_CALL(dgetrs) C++ call in R-devel
In my package, I have the method:
/** Solve equations transpose(P)w = r. */
int LASolveT(MatSimple &P, MatSimple &w, const
MatSimple &r) {
int rows = P.rows;
int nrhs = 1;
int lda = rows;
int ldb = rows;
int info = -1;
MatSimple ipivot(1,rows);
w.Inject(r);// copy r to w;
F77_CALL(dgetrf)(&rows, &rows, &P(0,0), &lda, &ipivot(0,0), &info);
if (info!=0) {
cout << "Error in LASolve (dgetrf). Info=" << info << endl;
return 1;
}
F77_CALL(dgetrs)("T", &rows, &nrhs, &P(0,0), &lda, &ipivot(0,0),
&w(0,0), &ldb, &info);("T", &rows, &nrhs, &P(0,0), &lda, &ipivot(0,0),
&w(0,0), &ldb, &info);
if (info!=0) {
cout << "Error in LASolve (dgetrs). Info=" << info << endl;
return 1;
}
return 0;
}
When compiling the package on using R-devel the error for F77_CALL(dgetrs)
occur:
matalg.h:67:25: error: too few arguments to function ‘void dgetrs_(const
char*, const int*, const int*, const double*, const int*, const int*,
double*, const int*, int*, size_t)’
67 | F77_CALL(dgetrs)("T", &rows, &nrhs, &P(0,0), &lda,
&ipivot(0,0), &w(0,0), &ldb, &info);
It works in R-release and I guess it have something to do with
https://cran.r-project.org/doc/manuals/r-devel/NEWS.html and LAPACK.
Any hints on how to get it to work for both R-release and R-devel.
Best
Lars Relund
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] .C and .Call: convolve example not working
Thanks a lot for the help Thomas and William, that solved it. I am
used to programming Java and didn't think of checking my print
function.
Regards Lars
2011/8/11 Thomas Lumley :
> On Fri, Aug 12, 2011 at 1:06 AM, Lars Wißler wrote:
>> Dear R users,
>>
>> I want to call C code via the .C or .Call interface. This works fine
>> with integer values but using doubles the array received in C will be
>> set to zeros.
>>
>> I have tried the convolve examples (Writing R extensions, chapter 5.2)
>> and still the resulting array consists of zeros.
>>
>> My code (shortened for my purposes. Original did not work either):
>>
>> convolve.r
>>
>> a=array(1, c(4))
>> b=array(3, c(4))
>>
>> print(a)
>> print(b)
>>
>> .C("convolve",
>> as.double(a),
>> as.double(b))
>>
>>
>> convolve.c
>>
>> void convolve(double* a, double* b){
>> int i;
>>
>> for(i=0;i<4;i++) Rprintf("a: %d", a[i]);
>> for(i=0;i<4;i++) Rprintf("b: %d", b[i]);
>> }
>
>
> The C code here is wrong for two reasons. Firstly, there's no
> declaration for Rprintf(), because you don't include the header file.
>
> Secondly, you're using %d to print, which means you are telling the
> compiler you're passing ints to Rprintf(), but you are actually
> passing doubles.
>
> When I fix these problems the code works for me.
>
> -thomas
>
> Thomas Lumley
> Professor of Biostatistics
> University of Auckland
>
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] .C and .Call: convolve example not working
Dear R users,
I want to call C code via the .C or .Call interface. This works fine
with integer values but using doubles the array received in C will be
set to zeros.
I have tried the convolve examples (Writing R extensions, chapter 5.2)
and still the resulting array consists of zeros.
My code (shortened for my purposes. Original did not work either):
convolve.r
a=array(1, c(4))
b=array(3, c(4))
print(a)
print(b)
.C("convolve",
as.double(a),
as.double(b))
convolve.c
void convolve(double* a, double* b){
int i;
for(i=0;i<4;i++) Rprintf("a: %d", a[i]);
for(i=0;i<4;i++) Rprintf("b: %d", b[i]);
}
--
ouput:
[1] "starting C code..."
[1] 1 1 1 1
[1] 3 3 3 3
a: 30467528
a: 0
a: 0
a: 0
b: 0
b: 0
b: 0
b: 0
Any suggestions as to why this is happening and what I am doing wrong
would be much appreciated. I have tried .Call with conversion SEXP to
double with REAL(a)/REAL(b) with the same results (the entry first
entry of a has a different number, but is huge as well). I am quite
astonished with the results I am getting.
Thanks and regards,
Lars Wissler
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Adressing Problems: R with Fortran and OpenMP
Hello,
I am programming an R program with nested Fortran calls for
calculations and OpenMP for parallelization. I am getting a changing
error corresponding to memory addressing problems, when using a 64-bit
system. Using a 32-bit System the application runs without problems.
The errors on 64-bit range from null-pointer failures, over
segmentation faults, over stack imbalances (changing differences and I
am not using C/C++) to finishing without exception but with wrong
values. Sometimes it even works correctly on 64-bit, mostly when
executing a second time within the same R session. Sometimes an
endless loop "Error: bad target context--should NEVER happen; please
bug.report() [R_run_onexits]" appears.
The problem seems to be platform independent. I have tried windows 7,
windows vista and open suse 11.3. (x86-64). Evaluation with valgrid
reveals a major possible memory leak, though the leak appears on
32-bit systems as well, just no errors. I am using a gfortran 4.5.0
x86-64 compiler and R version 2.12.
valgrid log extract:
==22989== 25,559,200 bytes in 4 blocks are possibly lost in loss
record 5,678 of 5,678
==22989==at 0x4C26C3A: malloc (in
/usr/lib64/valgrind/vgpreload_memcheck-amd64-linux.so)
==22989==by 0x4F39907: Rf_allocVector (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4EDAF96: duplicate1 (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4FC204E: R_subassign3_dflt (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4FC24A2: do_subassign3 (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4F00B4A: Rf_eval (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4F0346F: do_set (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4F00B4A: Rf_eval (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4F02EB1: applydefine (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4F00B4A: Rf_eval (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4F035EB: do_begin (in /usr/lib64/R/lib/libR.so)
==22989==by 0x4F00B4A: Rf_eval (in /usr/lib64/R/lib/libR.so)
==22989==
==22989== LEAK SUMMARY:
==22989==definitely lost: 82 bytes in 1 blocks
==22989==indirectly lost: 0 bytes in 0 blocks
==22989== possibly lost: 109,720,966 bytes in 26,330 blocks
==22989==still reachable: 23,101,045 bytes in 5,105 blocks
==22989== suppressed: 0 bytes in 0 blocks
All pointers in Fortran are explicitly defined with integer*4 and
real*8 as double.
I am really lost in this, because i just dont know where to start and
stop looking. It is obvious to me, that there is some kind of memory
adressing problem related to 64-bit architecture but since I dont know
if its related to R or Fortran or OpenMp or a combination of those, it
is very hard to find. Also the program is part of a library with 40+
files which interact, so I it would be really hard and time consuming
to cut the program down to a size, where the error will be reproduced
and still managable.
Any help, ideas, suggestions as to what to do, where to look and what
to try would be very welcome. I have been trying to solve this problem
for nearly two weeks and read everything I could find regarding
x86-64, R, Fortran, OpenMP and memory issues. I could post more and
more specific information regarding the errors, but then the
description would get even bigger. So if I need to supply more
information, please tell me and I will do so.
Regards
Lars
Following are the code snippets for the Fortran call and the entrance
to the Fortran program with OpenMp definition. If the program fails
with an statement about where it failed (i.e. segmentation fault),
then it gives this call as place. But since I only get R errors and
not Fortran errors, the error might actually occur anywhere in
Fortran.
z <- .Fortran("nlrdtirg",
as.integer(si),
as.integer(ngrad),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.logical(mask),
as.double(object@btb),
as.double(sdcoef),
th0=as.double(s0),
D=double(6*prod(ddim)),
as.integer(200),
as.double(1e-6),
res=double(ngrad*prod(ddim)),
rss=double(prod(ddim)),
double(ngrad*num_threads),
as.integer(num_threads),
PACKAGE="dti",DUP=TRUE)
subroutine nlrdtirg(s,nb,n1,n2,n3,mask,b,sdcoef,th0,D,niter,eps,
1res,rss,varinv,nt)
use omp_lib
implicit logical*4 (a-z)
integer*4 nb,n1,n2,n3,s(nb,n1,n2,n3),niter,nt,tid
logical mask(n1,n2,n3)
real*8 D(6,n1,n2,n3),b(6,nb),res(nb,n1,n2,n3),
1th0(n1,n2,n3),eps,rss(n1,n2,n3),sdcoef(4),varinv(nt*nb)
integer*4 i1,i2,i3,j
DO i3=1,n3
DO i2=1,n2
C$OMP PARALLEL DEFAULT(NONE)
C$OMP& SHARED(mask,s,b,sdcoef,th0,D,res,rss,varinv,nb,niter,eps)
C$OMP& FIRSTPRIVATE(i2,i3,n1)
C$OMP&
[Rd] unexpected behavior of rpart 3.1-43, loss matrix
Hi,
I just noticed that rpart behaves unexpectecly, when performing
classification learning and specifying a loss matrix.
if the response variable y is a factor and if not all levels of the
factor occur in the observations, rpart exits with an error:
> df=data.frame(attr=1:5,class=factor(c(2,3,1,5,3),levels=1:6))
> rpart(class~attr,df,parms=list(loss=matrix(0,6,6)))
Error in (get(paste("rpart", method, sep = ".")))(Y, offset, parms, wt)
: Wrong length for loss matrix
note that while the levels of the factor range from 1:6, for the
concrete obseration data, only levels 1, 2, 3, 5 do occur.
the error is caused by the code of rpart.class:
fy <- as.factor(y)
y <- as.integer(fy)
numclass <- max(y[!is.na(y)])
...
temp2 <- parms$loss
if (length(temp2) != numclass^2)
stop("Wrong length for loss matrix")
for the example, numclass is set to 5 instead of 6.
while for that small example, it may be discussable whether or not
numclass should be 6, consider a set of data for that the response
variable has a certain range. Then, it may be the case that for some
data, not all levels of the response variable do occur. at the same
time, it is desirable to use the same loss matrix when training a
deicision tree from the data.
having said that, i am very happy with the rpart package and with its
high configurability.
best regards
lars
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Using R from Java
Hi, Take a look at the rJava package. It includes JRI that lets you call R from Java. From the README: This package contains code that is necessary to run R as a single thread of a Java application. It provides callback that make it possible to run R in REPL mode thus giving the Java application full access to the console. Currently the API is very, very low-level, comparable to the C level interface to R. Convenience methods for mid to high-level are planned, but not implemented yet. Good luck, Lars Marzio Sala wrote: Hello, I am interesting in using R from a web application, for basic statistics and plots. The server is Java-based (tomcat). The simplest solution is a system call that generates the text or the image, then the servlet forwards the output. This can be done from any language, but it is quite inelegant and slow for the initialization time. Is there any package or approach for accessing R from a Java servlet you can suggest? Thanks in advance for any suggestion. Regards, -Marzio __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] R-Project logo in SVG format
Hi, I'm looking for a version of the R-Project logo in SVG format. I've found the bitmapped versions, http://developer.r-project.org/Logo/ but would prefer a scalable version as it usually looks better when printed. Where may I find one? -Lars Lars Nooden ([EMAIL PROTECTED]) On the Internet, nobody knows you're a dog ... ... until you start barking. __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] (PR#7951) DispatchOrEval missing in do_isfinite and do_isinfinite
Hi,
I ran into another internal function that is missing S4 dispatch. It is
the binary operator ":". Looking at the code, I see that it is actually
a common problem. Other candidates are operators like "~", "&&", "||"
and functions like: "length<-", "row", "col", "unlist", "cbind", etc. It
would for instance be nice to be able to write a matrix class that has
the same operators and functions as the built-in class. In general, I
think that all the operators and functions associates with built-in
types like vectors, lists, matrices and data frames should have S4 dispatch.
Thanks,
Lars
lars wrote:
> Hi,
>
> OK, if you try to explicitly make them generic, you are told that they
> are implicitly already generic:
>
> > setGeneric("is.finite", function(from, ...)
> standardGeneric("is.finite"))
> Error in setGeneric("is.finite", function(from, ...)
> standardGeneric("is.finite")) :
>"is.finite" is a primitive function; methods can be defined, but
> the generic function is implicit, and can't be changed.
>
> If you query about its genericness before you define you own generic,
> you get:
>
> > isGeneric("is.finite")
> [1] FALSE
>
> But after you define you own generic, you get:
>
> > setMethod("is.finite", signature(x="TS"),
> + function(x) {
> + Data(x) = callNextMethod()
> + x
> + })
> [1] "is.finite"
>
> > isGeneric("is.finite")
> [1] TRUE
>
> This all makes some sense, but I am not familiar enough with he
> internals to explain exactly why it is done this way. I think you will
> fine that 'is.nan' behave exactly the same way.
>
> Thanks,
> Lars
>
>
> Prof Brian Ripley wrote:
>
>> These functions are not generic according to the help page.
>> The same page says explicitly that is.nan is generic.
>>
>> Where did you get the (false) idea that they were generic?
>>
>> On Thu, 16 Jun 2005 [EMAIL PROTECTED] wrote:
>>
>>> Full_Name: Lars Hansen
>>> Version: 2.1.0
>>> OS: SunOS 5.8
>>> Submission from: (NULL) (207.66.36.189)
>>>
>>>
>>> Hi,
>>>
>>> S4 method displacth does not work for the two generic functions
>>> 'is.finite' and 'is.infinite'. It turns out that the C functions
>>> 'do_isfinite' and 'do_isinfinite' in src/main/coerce.c are missing a
>>> call to 'DispatchOrEval' (see do_isnan). Added in the call fixed the
>>> problem. My functions no look like this:
>>>
>>> Form coerce.c:
>>>
>>> SEXP do_isfinite(SEXP call, SEXP op, SEXP args, SEXP rho)
>>> {
>>>SEXP ans, x, names, dims;
>>>int i, n;
>>>
>>>if (DispatchOrEval(call, op, "is.finite", args, rho, &ans, 1, 1))
>>>return(ans);
>>>
>>>checkArity(op, args);
>>>...
>>>
>>> SEXP do_isinfinite(SEXP call, SEXP op, SEXP args, SEXP rho)
>>> {
>>>SEXP ans, x, names, dims;
>>>double xr, xi;
>>>int i, n;
>>>
>>>if (DispatchOrEval(call, op, "is.infinite", args, rho, &ans, 1, 1))
>>>return(ans);
>>>
>>>checkArity(op, args);
>>>...
>>
>>
>>
>
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] (PR#7951) DispatchOrEval missing in do_isfinite and do_isinfinite
Hi,
OK, if you try to explicitly make them generic, you are told that they
are implicitly already generic:
> setGeneric("is.finite", function(from, ...) standardGeneric("is.finite"))
Error in setGeneric("is.finite", function(from, ...)
standardGeneric("is.finite")) :
"is.finite" is a primitive function; methods can be defined, but
the generic function is implicit, and can't be changed.
If you query about its genericness before you define you own generic,
you get:
> isGeneric("is.finite")
[1] FALSE
But after you define you own generic, you get:
> setMethod("is.finite", signature(x="TS"),
+ function(x) {
+ Data(x) = callNextMethod()
+ x
+ })
[1] "is.finite"
> isGeneric("is.finite")
[1] TRUE
This all makes some sense, but I am not familiar enough with he
internals to explain exactly why it is done this way. I think you will
fine that 'is.nan' behave exactly the same way.
Thanks,
Lars
Prof Brian Ripley wrote:
> These functions are not generic according to the help page.
> The same page says explicitly that is.nan is generic.
>
> Where did you get the (false) idea that they were generic?
>
> On Thu, 16 Jun 2005 [EMAIL PROTECTED] wrote:
>
>> Full_Name: Lars Hansen
>> Version: 2.1.0
>> OS: SunOS 5.8
>> Submission from: (NULL) (207.66.36.189)
>>
>>
>> Hi,
>>
>> S4 method displacth does not work for the two generic functions
>> 'is.finite' and 'is.infinite'. It turns out that the C functions
>> 'do_isfinite' and 'do_isinfinite' in src/main/coerce.c are missing a
>> call to 'DispatchOrEval' (see do_isnan). Added in the call fixed the
>> problem. My functions no look like this:
>>
>> Form coerce.c:
>>
>> SEXP do_isfinite(SEXP call, SEXP op, SEXP args, SEXP rho)
>> {
>>SEXP ans, x, names, dims;
>>int i, n;
>>
>>if (DispatchOrEval(call, op, "is.finite", args, rho, &ans, 1, 1))
>>return(ans);
>>
>>checkArity(op, args);
>>...
>>
>> SEXP do_isinfinite(SEXP call, SEXP op, SEXP args, SEXP rho)
>> {
>>SEXP ans, x, names, dims;
>>double xr, xi;
>>int i, n;
>>
>>if (DispatchOrEval(call, op, "is.infinite", args, rho, &ans, 1, 1))
>>return(ans);
>>
>>checkArity(op, args);
>>...
>
>
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] DispatchOrEval missing in do_isfinite and do_isinfinite (PR#7951)
Full_Name: Lars Hansen
Version: 2.1.0
OS: SunOS 5.8
Submission from: (NULL) (207.66.36.189)
Hi,
S4 method displacth does not work for the two generic functions 'is.finite' and
'is.infinite'. It turns out that the C functions 'do_isfinite' and
'do_isinfinite' in src/main/coerce.c are missing a call to 'DispatchOrEval' (see
do_isnan). Added in the call fixed the problem. My functions no look like this:
Form coerce.c:
SEXP do_isfinite(SEXP call, SEXP op, SEXP args, SEXP rho)
{
SEXP ans, x, names, dims;
int i, n;
if (DispatchOrEval(call, op, "is.finite", args, rho, &ans, 1, 1))
return(ans);
checkArity(op, args);
...
SEXP do_isinfinite(SEXP call, SEXP op, SEXP args, SEXP rho)
{
SEXP ans, x, names, dims;
double xr, xi;
int i, n;
if (DispatchOrEval(call, op, "is.infinite", args, rho, &ans, 1, 1))
return(ans);
checkArity(op, args);
...
Thanks you,
Lars
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
