[R-sig-Debian] Openblas?

[Göran Broström]

Here is the reproducible example (very educational ...):
----------------------------------------------------------
$ wget http://ehar.se/data/reex_0.1.1.tar.gz
$ R CMD INSTALL reex_0.1.1.tar.gz

{with "status = 1: blas/libblas.so.3, and lapack..}
$ R

 > library("reex")
 > system.time(res <- coxfunk(beta = 1, X = X, rs = rs, what = 2))
#   user  system elapsed
#  0.093   0.000   0.093

 > q()

{Change to "status 2; openmp/libblas.so.3", and lapack..}

$ R
 > library(reex)
 > system.time(res <- coxfunk(beta = 1, X = X, rs = rs, what = 2))
#   user  system elapsed
# 72.050   1.006   6.123
 >
 > system.time(res <- coxfunk(beta = 1, X = X, rs = rs, what = 1))
#   user  system elapsed
#  0.088   0.000   0.088
-----------------------------------------------------------
Comment:
## what = 1 calculates loglik and score, what = 2 in addition hessian
## "Extra" code when what = 2:

if (*what >= 2){ /* Second derivatives: */
     F77_CALL(dsyr)(&up, p, (wsc + i),
               (x + (*p) * i), &ione,
	      sumd2score, p FCONE);
}


and

if (*what >= 2){
     alpha = -alpha;
     F77_CALL(daxpy)(&p2, &alpha, sumd2score, &ione,
		    d2loglik, &ione);
     alpha = -alpha / sumscore;
     F77_CALL(dsyr)(&up, p, &alpha, sumdscore, &ione,
		   d2loglik, p FCONE);
}

Full C and R code in package.

Göran

On 2020-07-15 16:32, Dirk Eddelbuettel wrote:
> 
> On 15 July 2020 at 16:13, Göran Broström wrote:
> | On 2020-07-15 14:36, Dirk Eddelbuettel wrote:
> | >
> | > Göran,
> | >
> | > This is not an easy email to reply to because it _contains nothing
> | > reproducible_.
> |
> | Thanks Dirk,
> |
> | Sorry about that, but my real question was (see below): "Is the problem
> | that openblas uses C versions of blas?" That is, do I need to change
> |
> | F77_CALL(name)(...);
> |
> | to
> |
> | cblas_name(...);
> |
> | everywhere? And if so, is this really a good idea with old code?
> 
> I don't think so. At the end of the day it comes from "higher up" (say,
> crossprod()) and is just passed down.
> 
> Remember, at the end it's all assembler :)
> 
> | I'll try to extract a reproducible example from the package (eha) where
> | I run it.
> 
> +1
> 
> Dirk
>