Thank you very much, Tomas, now it's clear and I'll see what to do with that
knowledge!
Luc
________________________________________
Van: Tomas Kalibera <tomas.kalib...@gmail.com>
Verzonden: donderdag 5 december 2024 14:39
Aan: Luc De Wilde <luc.dewi...@ugent.be>; r-package-devel@r-project.org
<r-package-devel@r-project.org>
CC: Yves Rosseel <yves.ross...@ugent.be>
Onderwerp: Re: [R-pkg-devel] Cannot create C code with acceptable performance
with respect to internal R command.
On 12/5/24 14:21, Luc De Wilde wrote:
> Dear package developers,
>
> in creating a package lavaanC for use in lavaan, I need to perform some
> matrix computations involving matrix products and crossproducts. As far as I
> see I cannot directly call the C code in the R core. So I copied the code in
> the R core, but the same C/C++ code in a package is 2.5 à 3 times slower than
> executed directly in R :
>
> C code in package :
> SEXP prod0(SEXP mat1, SEXP mat2) {
> SEXP u1 = Rf_getAttrib(mat1, R_DimSymbol);
> int m1 = INTEGER(u1)[0];
> int n1 = INTEGER(u1)[1];
> SEXP u2 = Rf_getAttrib(mat2, R_DimSymbol);
> int m2 = INTEGER(u2)[0];
> int n2 = INTEGER(u2)[1];
> if (n1 != m2) Rf_error("matrices not conforming");
> SEXP retval = PROTECT(Rf_allocMatrix(REALSXP, m1, n2));
> double* left = REAL(mat1);
> double* right = REAL(mat2);
> double* ret = REAL(retval);
> double werk = 0.0;
> for (int j = 0; j < n2; j++) {
> for (int i = 0; i < m1; i++) {
> werk = 0.0;
> for (int k = 0; k < n1; k++) werk += (left[i + m1 * k] * right[k +
> m2 * j]);
> ret[j * m1 + i] = werk;
> }
> }
> UNPROTECT(1);
> return retval;
> }
>
> Test script :
> m1 <- matrix(rnorm(300000), nrow = 60)
> m2 <- matrix(rnorm(300000), ncol = 60)
> print(microbenchmark::microbenchmark(
> m1 %*% m2, .Call("prod0", m1, m2), times = 100
> ))
>
> Result on my pc:
> Unit: milliseconds
> expr min lq mean median uq max
> neval
> m1 %*% m2 10.5650 10.8967 11.13434 10.9449 11.02965 15.8397
> 100
> .Call("prod0", m1, m2) 29.3336 30.7868 32.05114 31.0408 33.85935 45.5321
> 100
>
>
> Can anyone explain why the compiled code in the package is so much slower
> than in R core?
By default, R would use BLAS, not the simple algorithm above. See
?options, look for "matprod" for more information on how to select an
algorithm. The code is then in array.c, function matprod().
> and
>
> Is there a way to improve the performance in R package?
One option is to use BLAS.
Best
Tomas
>
>
> Best regards,
>
> Luc De Wilde
>
>
>
> ______________________________________________
> R-package-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-package-devel
______________________________________________
R-package-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-package-devel
