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
_______________________________________________
R-SIG-Debian mailing list
R-SIG-Debian@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-debian
