CRAN Package Check Results for Package expm

Last updated on 2019-03-20 11:52:33 CET.

Flavor Version Tinstall Tcheck Ttotal Status Flags
r-devel-linux-x86_64-debian-clang 0.999-3 9.55 91.42 100.97 OK
r-devel-linux-x86_64-debian-gcc 0.999-3 8.01 71.94 79.95 OK
r-devel-linux-x86_64-fedora-clang 0.999-3 120.62 OK
r-devel-linux-x86_64-fedora-gcc 0.999-3 113.69 OK
r-devel-windows-ix86+x86_64 0.999-3 28.00 181.00 209.00 OK
r-patched-linux-x86_64 0.999-3 7.01 86.04 93.05 OK
r-patched-solaris-x86 0.999-3 176.20 OK
r-release-linux-x86_64 0.999-3 6.64 85.62 92.26 OK
r-release-windows-ix86+x86_64 0.999-3 27.00 182.00 209.00 OK
r-release-osx-x86_64 0.999-3 OK
r-oldrel-windows-ix86+x86_64 0.999-3 21.00 155.00 176.00 ERROR
r-oldrel-osx-x86_64 0.999-3 OK

Check Details

Version: 0.999-3
Check: running examples for arch ‘i386’
Result: ERROR
    Running examples in 'expm-Ex.R' failed
    The error most likely occurred in:
    
    > ### Name: expAtv
    > ### Title: Compute Matrix Exponential exp(A t) * v directly
    > ### Aliases: expAtv expmv
    > ### Keywords: algebra math
    >
    > ### ** Examples
    >
    > source(system.file("demo", "exact-fn.R", package = "expm"))
    > ##-> rnilMat() ; xct10
    > set.seed(1)
    > (s5 <- Matrix(m5 <- rnilMat(5))); v <- c(1,6:9)
    5 x 5 sparse Matrix of class "dtCMatrix"
    Error in isFALSE(suppRows) : could not find function "isFALSE"
    Calls: <Anonymous> -> <Anonymous> -> printSpMatrix2
    Execution halted
Flavor: r-oldrel-windows-ix86+x86_64

Version: 0.999-3
Check: running examples for arch ‘x64’
Result: ERROR
    Running examples in 'expm-Ex.R' failed
    The error most likely occurred in:
    
    > ### Name: expAtv
    > ### Title: Compute Matrix Exponential exp(A t) * v directly
    > ### Aliases: expAtv expmv
    > ### Keywords: algebra math
    >
    > ### ** Examples
    >
    > source(system.file("demo", "exact-fn.R", package = "expm"))
    > ##-> rnilMat() ; xct10
    > set.seed(1)
    > (s5 <- Matrix(m5 <- rnilMat(5))); v <- c(1,6:9)
    5 x 5 sparse Matrix of class "dtCMatrix"
    Error in isFALSE(suppRows) : could not find function "isFALSE"
    Calls: <Anonymous> -> <Anonymous> -> printSpMatrix2
    Execution halted
Flavor: r-oldrel-windows-ix86+x86_64

Version: 0.999-3
Check: running tests for arch ‘i386’
Result: ERROR
     Running 'Frechet-test.R' [2s]
     Running 'bal-ex.R' [3s]
     Running 'ex.R' [3s]
     Running 'ex2.R' [4s]
     Running 'exact-ex.R' [6s]
    Running the tests in 'tests/exact-ex.R' failed.
    Complete output:
     > #### Examples where we know the result "exactly"
     >
     > library(expm)
     Loading required package: Matrix
    
     Attaching package: 'expm'
    
     The following object is masked from 'package:Matrix':
    
     expm
    
     >
     > options(digits = 4, width = 90, keep.source = FALSE)
     >
     > mSource <- function(file, ...)
     + source(system.file(file, ..., package = "expm", mustWork=TRUE))
     > mSource("test-tools.R")## -> assertError(), rMat(), .. doExtras
     Loading required package: tools
     > mSource("demo", "exact-fn.R")
     > doExtras
     [1] FALSE
     >
     > re.nilA3 <- function(xyz, EXPMlist)
     + {
     + stopifnot(is.list(EXPMlist))
     + r <- do.call(nilA3, as.list(xyz))
     + sapply(EXPMlist, function(Efn) relErr(r$expA, Efn(r$A)))
     + }
     >
     > re.facMat <- function(n, EXPMlist, rFUN = rnorm, ...)
     + {
     + stopifnot(is.list(EXPMlist))
     + r <- facMat(n, rFUN, ...)
     + vapply(EXPMlist, function(EXPM) {
     + ct <- system.time(E <- EXPM(r$A), gc = FALSE)[[1]]
     + c(relErr = relErr(r$expA, E), c.time = ct)
     + }, double(2))
     + }
     >
     > re.m2ex3 <- function(eps, EXPMlist)
     + {
     + stopifnot(is.list(EXPMlist))
     + r <- m2ex3(eps)
     + sapply(EXPMlist, function(EXPM) relErr(r$expA, EXPM(r$A)))
     + }
     >
     > ## check 1x1 matrices
     > stopifnot(
     + ## these had failed before 2017-03-28 (= Liselotte's 58-th birthday):
     + all.equal(as.matrix(sqrtm(matrix(4))), matrix(2)) ,
     + all.equal(as.matrix(logm (matrix(pi))), matrix(log(pi))) ,
     + ## these had "always" worked :
     + all.equal(as.matrix(expm (matrix(0))), matrix(1)) ,
     + all.equal(as.matrix(expm (matrix(1))), matrix(exp(1)))
     + )
     >
     >
     > set.seed(321)
     > re <- replicate(1000,
     + c(re.nilA3(rlnorm(3),list(function(x)expm(x,"Pade"))),
     + re.nilA3(rnorm(3), list(function(x)expm(x,"Pade")))))
     >
     > summary(t(re))
     V1 V2
     Min. :0.00e+00 Min. :0.00e+00
     1st Qu.:0.00e+00 1st Qu.:0.00e+00
     Median :2.60e-17 Median :1.60e-17
     Mean :4.40e-16 Mean :1.90e-16
     3rd Qu.:6.20e-17 3rd Qu.:4.10e-17
     Max. :2.12e-14 Max. :4.68e-15
     > stopifnot(rowMeans(re) < 1e-15,
     + apply(re, 1, quantile, 0.80) < 1e-16,
     + apply(re, 1, quantile, 0.90) < 2e-15,
     + apply(re, 1, max) < c(4e-14, 6e-15))
     >
     > showProc.time()
     Time elapsed: 1.46 0 1.46
     >
     > ## Check *many* random nilpotent 3 x 3 matrices:
     > set.seed(321)
     > RE <- replicate(1000,
     + c(re.nilA3(rlnorm(3), list(function(x) expm(x, "Ward77"))),
     + re.nilA3(rnorm(3), list(function(x) expm(x, "Ward77")))))
     > stopifnot(rowMeans(RE) < 1e-15,
     + apply(RE, 1, quantile, 0.80) < 1e-16,
     + apply(RE, 1, quantile, 0.90) < 2e-15,
     + apply(RE, 1, max) < c(4e-14, 6e-15))
     >
     > print(summary(t(RE)))
     V1 V2
     Min. :0.00e+00 Min. :0.00e+00
     1st Qu.:0.00e+00 1st Qu.:0.00e+00
     Median :0.00e+00 Median :0.00e+00
     Mean :1.56e-17 Mean :8.48e-18
     3rd Qu.:2.88e-17 3rd Qu.:1.13e-17
     Max. :1.85e-16 Max. :1.08e-16
     > epsC <- .Machine$double.eps
     > cat("relErr(expm(.,Pade)) - relErr(expm(.,Ward77)) in Machine_eps units:\n")
     relErr(expm(.,Pade)) - relErr(expm(.,Ward77)) in Machine_eps units:
     > print(summary(c(re - RE)) / epsC)
     Min. 1st Qu. Median Mean 3rd Qu. Max.
     -0.62 0.00 0.00 1.36 0.14 94.98
     > ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     > ## -0.6183442 0.0000000 0.0000000 1.3650410 0.1399719 94.9809161
     > ## nb-mm3; ditto lynne (both x64), 2014-09-11:
     > ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     > ## -0.8422 0.0000 0.0000 0.0725 0.1067 1.2205
     > ## 32-bit [i686, florence, Linx 3.14.8-100.fc19..]:
     > ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     > ## -0.62 0.00 0.00 1.36 0.14 95.93
     >
     >
     > showProc.time()
     Time elapsed: 0.6 0 0.6
     >
     > ###--- A second group --- where we know the diagonalization of A ---
     >
     > if(!require("Matrix"))
     + q('no')
     > ## ------ the rest really uses 'Matrix'
     > ##---> now use expm::expm() since Matrix has its own may mask the expm one
     > ## ^^^^^^^^^^^^
     >
     > ## rMat() relies on Matrix::rcond():
     > ## Now with the change default rcondMin, this "works"
     > R40 <- rMat(40)
     > R80 <- rMat(80)
     > showProc.time()
     Time elapsed: 0.07 0 0.07
     >
     > expm.safe.Eigen <- function(x, silent = FALSE) {
     + r <- try(expm::expm(x, "R_Eigen"), silent = silent)
     + if(inherits(r, "try-error")) NA else r
     + }
     >
     > ## the S4 generic
     > Matrix::expm
     standardGeneric for "expm" defined from package "Matrix"
    
     function (x)
     standardGeneric("expm")
     <environment: 0x05009cd4>
     Methods may be defined for arguments: x
     Use showMethods("expm") for currently available ones.
     > ## the dgeMatrix method:
     > expm.Matr.dge <- function(x) getDataPart(getMethod("expm", "dgeMatrix"))(Matrix::..2dge(x))
     > expmList <-
     + list(Matr = Matrix::expm,
     + Matr.d = expm.Matr.dge,
     + Ward = function(x) expm::expm(x, "Ward77"),
     + s.P.s = function(x) expm::expm(x, "Pade"),
     + s.P.sO= function(x) expm::expm(x, "PadeO"),
     + s.P.sRBS= function(x) expm::expm(x, "PadeRBS"),
     + sPs.H08.= function(x) expm:: expm.Higham08(x, balancing=FALSE),
     + sPs.H08b= function(x) expm:: expm.Higham08(x, balancing= TRUE),
     + AmHi09.06= function(x) expm:::expm.AlMoHi09(x, p = 6),
     + AmHi09.08= function(x) expm:::expm.AlMoHi09(x, p = 8),
     + AmHi09.10= function(x) expm:::expm.AlMoHi09(x, p = 10),
     + AmHi09.12= function(x) expm:::expm.AlMoHi09(x, p = 12),
     + AmHi09.13= function(x) expm:::expm.AlMoHi09(x, p = 13),
     + s.T.s = function(x) expm::expm(x, "Taylor"),
     + s.T.sO= function(x) expm::expm(x, "TaylorO"),
     + Eigen = expm.safe.Eigen,
     + hybrid= function(x) expm::expm(x, "hybrid")
     + )
     >
     >
     > set.seed(12)
     > fRE <- replicate(if(doExtras) 100 else 20,
     + re.facMat(20, expmList))
     > cat("Number of correct decimal digits for facMat(20, rnorm):\n")
     Number of correct decimal digits for facMat(20, rnorm):
     > summary(-log10(t(fRE["relErr",,])))
     Matr Matr.d Ward s.P.s s.P.sO
     Min. :13.7 Min. :13.7 Min. :13.7 Min. :12.8 Min. :12.9
     1st Qu.:14.2 1st Qu.:14.2 1st Qu.:14.2 1st Qu.:13.3 1st Qu.:13.3
     Median :14.3 Median :14.3 Median :14.3 Median :13.4 Median :13.4
     Mean :14.3 Mean :14.3 Mean :14.3 Mean :13.4 Mean :13.4
     3rd Qu.:14.5 3rd Qu.:14.5 3rd Qu.:14.5 3rd Qu.:13.5 3rd Qu.:13.6
     Max. :14.7 Max. :14.7 Max. :14.7 Max. :13.6 Max. :13.7
     s.P.sRBS sPs.H08. sPs.H08b AmHi09.06 AmHi09.08
     Min. :14.0 Min. :14.4 Min. :14.4 Min. : 6.61 Min. : 9.78
     1st Qu.:14.3 1st Qu.:14.8 1st Qu.:14.9 1st Qu.: 7.69 1st Qu.:11.27
     Median :14.4 Median :15.0 Median :15.0 Median : 8.24 Median :11.72
     Mean :14.4 Mean :14.9 Mean :14.9 Mean : 8.47 Mean :12.19
     3rd Qu.:14.6 3rd Qu.:15.0 3rd Qu.:15.0 3rd Qu.: 9.41 3rd Qu.:13.55
     Max. :14.8 Max. :15.2 Max. :15.2 Max. :10.35 Max. :14.44
     AmHi09.10 AmHi09.12 AmHi09.13 s.T.s s.T.sO
     Min. :13.4 Min. :14.4 Min. :14.4 Min. :13.0 Min. :13.3
     1st Qu.:14.6 1st Qu.:14.7 1st Qu.:14.7 1st Qu.:13.4 1st Qu.:13.6
     Median :14.8 Median :14.9 Median :14.9 Median :13.5 Median :13.7
     Mean :14.7 Mean :14.8 Mean :14.8 Mean :13.5 Mean :13.8
     3rd Qu.:14.9 3rd Qu.:15.0 3rd Qu.:15.0 3rd Qu.:13.7 3rd Qu.:13.9
     Max. :15.0 Max. :15.1 Max. :15.1 Max. :13.9 Max. :14.1
     Eigen hybrid
     Min. :14.5 Min. :14.4
     1st Qu.:14.6 1st Qu.:14.6
     Median :14.7 Median :14.7
     Mean :14.7 Mean :14.7
     3rd Qu.:14.8 3rd Qu.:14.8
     Max. :15.0 Max. :15.0
     >
     >
     > ## Now look at that:
     > boxplot(t(fRE["relErr",,]), log="y", notch=TRUE, ylim = c(8e-16, 1e-8),
     + main = "relative errors for 'random' eigen-ok 20 x 20 matrix")
     Warning message:
     In bxp(list(stats = c(1.8789451450914e-15, 2.89597166164401e-15, :
     some notches went outside hinges ('box'): maybe set notch=FALSE
     >
     > showProc.time()
     Time elapsed: 0.75 0 0.76
     >
     > if(doExtras) {
     + str(rf100 <- replicate(20, re.facMat(100, expmList)))
     + print(1000*t(apply(rf100["c.time",,], 1, summary)))
     + ## lynne {Linux 2.6.34.7-56.fc13.x86_64 --- AMD Phenom II X4 925}:
     + ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     + ## Ward 23 24 24.5 24.4 25.0 25
     + ## s.P.s 107 109 109.0 109.0 109.0 112
     + ## s.P.sO 188 190 191.0 192.0 193.0 198
     + ## s.P.sRBS 17 18 19.0 18.9 19.2 21
     + ## sPs.H08. 15 17 18.0 17.6 18.0 19
     + ## sPs.H08b 18 18 19.0 23.4 20.0 107
     + ## s.T.s 44 45 45.0 45.6 46.0 48
     + ## s.T.sO 96 98 99.0 100.0 100.0 116
     + ## Eigen 18 19 20.0 24.4 21.0 109
     + ## hybrid 40 42 42.0 47.1 44.0 133
     +
     + ##--> take out the real slow ones for the subsequent tests:
     + `%w/o%` <- function(x, y) x[!x %in% y] #-- x without y
     + print(nms.swift <- names(expmList) %w/o%
     + c("s.P.s", "s.P.sO", "s.T.s", "s.T.sO"))
     + expmL.swift <- expmList[nms.swift]
     +
     + set.seed(18) ## 12 replicates is too small .. but then it's too slow otherwise:
     + rf400 <- replicate(12, re.facMat(400, expmL.swift))
     + print(1000*t(apply(rf400["c.time",,], 1, summary)))
     + ## lynne:
     + ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     + ## Ward 1740 1790 1830 1820 1860 1900
     + ## s.P.sRBS 1350 1420 1440 1430 1450 1460
     + ## sPs.H08. 1020 1030 1130 1140 1210 1290
     + ## sPs.H08b 1120 1130 1220 1220 1300 1390
     + ## Eigen 962 977 989 992 1000 1030
     + ## hybrid 2740 2800 2840 2840 2890 2910
     +
     + showProc.time()
     +
     + }## if(doExtras) only
     >
     > ## Now try an example with badly conditioned "random" M matrix...
     > ## ...
     > ## ... (not yet)
     >
     >
     > ### m2ex3() --- The 2x2 example with bad condition , see A3 in ./ex2.R
     >
     > RE <- re.m2ex3(1e-8, expmList)
     > sort(RE)# Ward + both sps.H08 are best; s.P.s fair, Eigen (and hybrid): ~1e-9
     sPs.H08. sPs.H08b s.P.sRBS Matr Matr.d Ward AmHi09.13 AmHi09.12 s.P.s
     0.000e+00 0.000e+00 5.584e-33 5.030e-17 5.030e-17 5.030e-17 1.509e-16 3.018e-16 2.163e-15
     s.T.s s.P.sO s.T.sO AmHi09.10 hybrid Eigen AmHi09.08 AmHi09.06
     2.163e-15 2.213e-15 2.213e-15 2.325e-12 1.402e-09 1.836e-09 4.819e-09 3.156e-06
     >
     > eps <- 10^-(1:18)
     > t.m2 <- t(sapply(eps, re.m2ex3, EXPMlist = expmList))
     Error in solve.default(V) :
     system is computationally singular: reciprocal condition number = 1.11022e-16
     Error in solve.default(V) :
     system is computationally singular: reciprocal condition number = 1.11022e-16
     Error in solve.default(V) :
     system is computationally singular: reciprocal condition number = 1.11022e-16
     > ## --> 3 error messages from solve(V), 5 error messages from try(. "R_Eigen" ...)
     >
     > showProc.time()
     Time elapsed: 0.34 0 0.34
     >
     > cbind(sort(apply(log(t.m2),2, median, na.rm=TRUE)))
     [,1]
     sPs.H08. -90.55
     sPs.H08b -90.55
     s.P.sRBS -83.28
     Matr -37.53
     Matr.d -37.53
     Ward -37.53
     AmHi09.13 -36.43
     AmHi09.12 -35.74
     s.P.s -33.79
     s.P.sO -33.79
     s.T.s -33.79
     s.T.sO -33.79
     AmHi09.10 -26.79
     hybrid -24.44
     Eigen -20.12
     AmHi09.08 -19.15
     AmHi09.06 -12.67
     > ## 'na.rm=TRUE' needed for Eigen which blows up for the last 3 eps
     > t.m2.ranks <- sort(rowMeans(apply(t.m2, 1, rank)))
     > cbind(signif(t.m2.ranks, 3))
     [,1]
     sPs.H08. 2.67
     sPs.H08b 2.67
     s.P.sRBS 3.08
     AmHi09.13 5.33
     Matr 5.56
     Matr.d 5.56
     Ward 5.56
     AmHi09.12 7.11
     s.T.s 10.60
     s.P.s 10.90
     s.T.sO 10.90
     s.P.sO 11.00
     hybrid 12.80
     AmHi09.10 13.70
     Eigen 14.10
     AmHi09.08 15.10
     AmHi09.06 16.40
     > ## lynne (x86_64, Linux 3.14.4-100; Intel i7-4765T), 2014-09:
     > ## sPs.H08. 2.67
     > ## sPs.H08b 2.67
     > ## s.P.sRBS 3.06
     > ## Ward 4.03
     > ## AmHi09.13 4.33 <<- still not close to H08 !
     > ## AmHi09.12 5.86
     > ## s.T.s 8.33
     > ## s.T.sO 8.33
     > ## s.P.s 9.11
     > ## s.P.sO 9.11
     > ## hybrid 10.80
     > ## AmHi09.10 11.70 << astonishingly bad
     > ## Eigen 12.60
     > ## AmHi09.08 13.10
     > ## AmHi09.06 14.40
     >
     > print(t.m2[, names(t.m2.ranks)[1:8]], digits = 3)
     sPs.H08. sPs.H08b s.P.sRBS AmHi09.13 Matr Matr.d Ward AmHi09.12
     [1,] 9.99e-17 9.99e-17 2.00e-16 3.51e-16 6.51e-16 6.51e-16 6.51e-16 5.52e-16
     [2,] 1.01e-16 1.01e-16 2.01e-16 1.01e-16 5.53e-16 5.53e-16 5.53e-16 2.51e-16
     [3,] 5.03e-17 5.03e-17 2.51e-16 1.01e-16 2.51e-16 2.51e-16 2.51e-16 2.01e-16
     [4,] 2.51e-16 2.51e-16 3.75e-25 1.51e-16 2.01e-16 2.01e-16 2.01e-16 2.51e-16
     [5,] 3.52e-16 3.52e-16 2.01e-16 1.01e-16 1.51e-16 1.51e-16 1.51e-16 1.01e-16
     [6,] 2.01e-16 2.01e-16 2.01e-16 2.01e-16 2.01e-16 2.01e-16 2.01e-16 3.52e-16
     [7,] 3.52e-16 3.52e-16 2.51e-16 7.15e-31 4.53e-16 4.53e-16 4.53e-16 2.01e-16
     [8,] 0.00e+00 0.00e+00 5.58e-33 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [9,] 4.36e-35 4.36e-35 8.73e-35 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [10,] 0.00e+00 0.00e+00 0.00e+00 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [11,] 5.33e-39 5.33e-39 5.33e-39 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [12,] 4.16e-41 4.16e-41 0.00e+00 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [13,] 0.00e+00 0.00e+00 6.50e-43 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [14,] 0.00e+00 0.00e+00 5.08e-45 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [15,] 3.97e-47 3.97e-47 3.97e-47 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [16,] 0.00e+00 0.00e+00 6.20e-49 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [17,] 0.00e+00 0.00e+00 4.84e-51 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     [18,] 3.78e-53 3.78e-53 0.00e+00 1.51e-16 5.03e-17 5.03e-17 5.03e-17 3.02e-16
     > ## ==> 1st class: H08 (both) and (but slightly better than) Ward
     > ## 2nd class s.T.s and s.P.s
     > ## "bad" : hybrid and Eigen
     >
     > ## ??? AmHi09 - methods, up to order = 10 are worse !
     > if(require(RColorBrewer)) {
     + ## Bcol <- brewer.pal(ncol(t.m2),"Dark2")
     + Bcol <- brewer.pal(min(9, ncol(t.m2)), "Set1")
     + Bcol <- Bcol[sqrt(colSums(col2rgb(Bcol)^2)) < 340]
     + ## FIXME: more colors ==> ~/R/MM/GRAPHICS/color-palettes.R
     + } else {
     + ## 7 from Dark2
     + ## Bcol <- c("#1B9E77", "#D95F02", "#7570B3", "#E7298A",
     + ## "#66A61E", "#E6AB02", "#A6761D")
     + ## Rather: those from "Set1"
     + Bcol <- c("#E41A1C", "#377EB8", "#4DAF4A",
     + "#984EA3", "#FF7F00", # too bright: "#FFFF33",
     + "#A65628", "#F781BF", "#999999")
     + }
     Loading required package: RColorBrewer
     >
     > matplot(eps, t.m2, type = "b", log = "xy", col=Bcol, lty = 1:9, pch=1:15,
     + axes=FALSE, frame = TRUE,
     + xlab = expression(epsilon), ylab = "relative error",
     + main = expression(expm(A, method == "*") *" relative errors for " *
     + A == bgroup("[", atop({-1} *" "* 1, {epsilon^2} *" "*{-1}), "]")))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     15 y values <= 0 omitted from logarithmic plot
     > legend("bottomright",colnames(t.m2), col=Bcol, lty = 1:9, pch=1:15,
     + inset = 0.02)
     > if(require("sfsmisc")) {
     + sfsmisc::eaxis(1, labels=FALSE)
     + sfsmisc::eaxis(1, at = eps[c(TRUE,FALSE)])
     + sfsmisc::eaxis(2)
     + ## sfsmisc::eaxis(2, labels=FALSE)
     + ## op <- par(las=2)
     + ## sfsmisc::eaxis(2, at = axTicks(2,log=TRUE)[c(TRUE,FALSE,FALSE)])
     + ## par(op)
     + } else {
     + axis(1)
     + axis(2)
     + }
     Loading required package: sfsmisc
     >
     > ## typical case:
     > ep <- 1e-10
     > (me <- m2ex3(ep))
     $A
     [,1] [,2]
     [1,] -1e+00 1
     [2,] 1e-20 -1
    
     $expA
     [,1] [,2]
     [1,] 3.679e-01 0.3679
     [2,] 3.679e-21 0.3679
    
     > me$expA * exp(1) ## the correct value ; numerically identical to simple matrix:
     [,1] [,2]
     [1,] 1e+00 1
     [2,] 1e-20 1
     > ## identical() not fulfilled e.g. on Solaris
     > stopifnot(all.equal(me$expA * exp(1),
     + rbind(c( 1, 1),
     + c(ep^2, 1)),
     + tolerance = 1e-14))
     > ## The relative error (matrices):
     > lapply(expmList, function(EXPM) 1 - EXPM(me$A)/me$expA)
     $Matr
     2 x 2 Matrix of class "dgeMatrix"
     [,1] [,2]
     [1,] 0 1.11e-16
     [2,] 0 0.00e+00
    
     $Matr.d
     2 x 2 Matrix of class "dgeMatrix"
     [,1] [,2]
     [1,] 0 1.11e-16
     [2,] 0 0.00e+00
    
     $Ward
     [,1] [,2]
     [1,] 0 1.11e-16
     [2,] 0 0.00e+00
    
     $s.P.s
     [,1] [,2]
     [1,] -2.22e-15 -2.22e-15
     [2,] -2.22e-15 -2.22e-15
     attr(,"accuracy")
     [1] 7.523e-37
    
     $s.P.sO
     [,1] [,2]
     [1,] -2.22e-15 -2.22e-15
     [2,] -2.22e-15 -2.22e-15
     attr(,"accuracy")
     [1] 7.523e-37
    
     $s.P.sRBS
     [,1] [,2]
     [1,] 0 0
     [2,] 0 0
    
     $sPs.H08.
     [,1] [,2]
     [1,] 0 0
     [2,] 0 0
    
     $sPs.H08b
     [,1] [,2]
     [1,] 0 0
     [2,] 0 0
    
     $AmHi09.06
     [,1] [,2]
     [1,] -1.153e-06 7.160e-06
     [2,] 7.160e-06 -1.153e-06
    
     $AmHi09.08
     [,1] [,2]
     [1,] -1.414e-09 1.163e-08
     [2,] 1.163e-08 -1.414e-09
    
     $AmHi09.10
     [,1] [,2]
     [1,] -5.702e-13 5.836e-12
     [2,] 5.836e-12 -5.702e-13
    
     $AmHi09.12
     [,1] [,2]
     [1,] 1.110e-16 5.551e-16
     [2,] 6.661e-16 1.110e-16
    
     $AmHi09.13
     [,1] [,2]
     [1,] 1.11e-16 1.11e-16
     [2,] 2.22e-16 1.11e-16
    
     $s.T.s
     [,1] [,2]
     [1,] -2.22e-15 -2.22e-15
     [2,] -2.22e-15 -2.22e-15
     attr(,"accuracy")
     [1] 0
    
     $s.T.sO
     [,1] [,2]
     [1,] -2.22e-15 -2.22e-15
     [2,] -2.22e-15 -2.22e-15
     attr(,"accuracy")
     [1] 0
    
     $Eigen
     [,1] [,2]
     [1,] 0.000e+00 -5.109e-07
     [2,] -1.048e-06 0.000e+00
    
     $hybrid
     [,1] [,2]
     [1,] 0.000e+00 -2.394e-07
     [2,] -1.561e-07 -2.220e-16
    
     >
     > ## Average number of correct digits [less "extreme" than plot above]
     > nDig <- sapply(expmList, function(EXPM) -log10(mean(abs(1 - EXPM(me$A)/me$expA))))
     > round(nDig, 2)
     Matr Matr.d Ward s.P.s s.P.sO s.P.sRBS sPs.H08. sPs.H08b AmHi09.06
     16.56 16.56 16.56 14.65 14.65 Inf Inf Inf 5.38
     AmHi09.08 AmHi09.10 AmHi09.12 AmHi09.13 s.T.s s.T.sO Eigen hybrid
     8.19 11.49 15.44 15.86 14.65 14.65 6.41 7.00
     > ## Ward s.P.s s.P.sO s.T.s s.T.sO Eigen hybrid
     > ## 16.26 14.65 14.65 14.65 14.65 6.20 6.39 [AMD Opteron 64-bit]
     > ## Inf 14.65 14.65 14.65 14.65 6.74 6.33 [Pentium-M (32-bit)]
     >
     > ###--- rnilMat() : random upper triangular (zero-diagonal) nilpotent n x n matrix
     >
     > set.seed(17)
     > m <- rnilMat(10)
     > (m. <- as(m,"sparseMatrix"))# for nicer printing - and more *below*
     10 x 10 sparse Matrix of class "dgCMatrix"
     Error in isFALSE(suppRows) : could not find function "isFALSE"
     Calls: <Anonymous> -> <Anonymous> -> printSpMatrix2
     Execution halted
Flavor: r-oldrel-windows-ix86+x86_64

Version: 0.999-3
Check: running tests for arch ‘x64’
Result: ERROR
     Running 'Frechet-test.R' [3s]
     Running 'bal-ex.R' [3s]
     Running 'ex.R' [3s]
     Running 'ex2.R' [4s]
     Running 'exact-ex.R' [7s]
    Running the tests in 'tests/exact-ex.R' failed.
    Complete output:
     > #### Examples where we know the result "exactly"
     >
     > library(expm)
     Loading required package: Matrix
    
     Attaching package: 'expm'
    
     The following object is masked from 'package:Matrix':
    
     expm
    
     >
     > options(digits = 4, width = 90, keep.source = FALSE)
     >
     > mSource <- function(file, ...)
     + source(system.file(file, ..., package = "expm", mustWork=TRUE))
     > mSource("test-tools.R")## -> assertError(), rMat(), .. doExtras
     Loading required package: tools
     > mSource("demo", "exact-fn.R")
     > doExtras
     [1] FALSE
     >
     > re.nilA3 <- function(xyz, EXPMlist)
     + {
     + stopifnot(is.list(EXPMlist))
     + r <- do.call(nilA3, as.list(xyz))
     + sapply(EXPMlist, function(Efn) relErr(r$expA, Efn(r$A)))
     + }
     >
     > re.facMat <- function(n, EXPMlist, rFUN = rnorm, ...)
     + {
     + stopifnot(is.list(EXPMlist))
     + r <- facMat(n, rFUN, ...)
     + vapply(EXPMlist, function(EXPM) {
     + ct <- system.time(E <- EXPM(r$A), gc = FALSE)[[1]]
     + c(relErr = relErr(r$expA, E), c.time = ct)
     + }, double(2))
     + }
     >
     > re.m2ex3 <- function(eps, EXPMlist)
     + {
     + stopifnot(is.list(EXPMlist))
     + r <- m2ex3(eps)
     + sapply(EXPMlist, function(EXPM) relErr(r$expA, EXPM(r$A)))
     + }
     >
     > ## check 1x1 matrices
     > stopifnot(
     + ## these had failed before 2017-03-28 (= Liselotte's 58-th birthday):
     + all.equal(as.matrix(sqrtm(matrix(4))), matrix(2)) ,
     + all.equal(as.matrix(logm (matrix(pi))), matrix(log(pi))) ,
     + ## these had "always" worked :
     + all.equal(as.matrix(expm (matrix(0))), matrix(1)) ,
     + all.equal(as.matrix(expm (matrix(1))), matrix(exp(1)))
     + )
     >
     >
     > set.seed(321)
     > re <- replicate(1000,
     + c(re.nilA3(rlnorm(3),list(function(x)expm(x,"Pade"))),
     + re.nilA3(rnorm(3), list(function(x)expm(x,"Pade")))))
     >
     > summary(t(re))
     V1 V2
     Min. :0.00e+00 Min. :0.00e+00
     1st Qu.:0.00e+00 1st Qu.:0.00e+00
     Median :2.15e-17 Median :1.47e-17
     Mean :3.15e-17 Mean :2.27e-17
     3rd Qu.:4.98e-17 3rd Qu.:3.61e-17
     Max. :2.71e-16 Max. :1.29e-16
     > stopifnot(rowMeans(re) < 1e-15,
     + apply(re, 1, quantile, 0.80) < 1e-16,
     + apply(re, 1, quantile, 0.90) < 2e-15,
     + apply(re, 1, max) < c(4e-14, 6e-15))
     >
     > showProc.time()
     Time elapsed: 1.85 0 1.86
     >
     > ## Check *many* random nilpotent 3 x 3 matrices:
     > set.seed(321)
     > RE <- replicate(1000,
     + c(re.nilA3(rlnorm(3), list(function(x) expm(x, "Ward77"))),
     + re.nilA3(rnorm(3), list(function(x) expm(x, "Ward77")))))
     > stopifnot(rowMeans(RE) < 1e-15,
     + apply(RE, 1, quantile, 0.80) < 1e-16,
     + apply(RE, 1, quantile, 0.90) < 2e-15,
     + apply(RE, 1, max) < c(4e-14, 6e-15))
     >
     > print(summary(t(RE)))
     V1 V2
     Min. :0.00e+00 Min. :0.00e+00
     1st Qu.:0.00e+00 1st Qu.:0.00e+00
     Median :0.00e+00 Median :0.00e+00
     Mean :1.49e-17 Mean :7.03e-18
     3rd Qu.:2.70e-17 3rd Qu.:5.20e-18
     Max. :1.87e-16 Max. :7.52e-17
     > epsC <- .Machine$double.eps
     > cat("relErr(expm(.,Pade)) - relErr(expm(.,Ward77)) in Machine_eps units:\n")
     relErr(expm(.,Pade)) - relErr(expm(.,Ward77)) in Machine_eps units:
     > print(summary(c(re - RE)) / epsC)
     Min. 1st Qu. Median Mean 3rd Qu. Max.
     -0.8429 0.0000 0.0000 0.0727 0.1065 1.2215
     > ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     > ## -0.6183442 0.0000000 0.0000000 1.3650410 0.1399719 94.9809161
     > ## nb-mm3; ditto lynne (both x64), 2014-09-11:
     > ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     > ## -0.8422 0.0000 0.0000 0.0725 0.1067 1.2205
     > ## 32-bit [i686, florence, Linx 3.14.8-100.fc19..]:
     > ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     > ## -0.62 0.00 0.00 1.36 0.14 95.93
     >
     >
     > showProc.time()
     Time elapsed: 0.9 0 0.9
     >
     > ###--- A second group --- where we know the diagonalization of A ---
     >
     > if(!require("Matrix"))
     + q('no')
     > ## ------ the rest really uses 'Matrix'
     > ##---> now use expm::expm() since Matrix has its own may mask the expm one
     > ## ^^^^^^^^^^^^
     >
     > ## rMat() relies on Matrix::rcond():
     > ## Now with the change default rcondMin, this "works"
     > R40 <- rMat(40)
     > R80 <- rMat(80)
     > showProc.time()
     Time elapsed: 0.06 0 0.06
     >
     > expm.safe.Eigen <- function(x, silent = FALSE) {
     + r <- try(expm::expm(x, "R_Eigen"), silent = silent)
     + if(inherits(r, "try-error")) NA else r
     + }
     >
     > ## the S4 generic
     > Matrix::expm
     standardGeneric for "expm" defined from package "Matrix"
    
     function (x)
     standardGeneric("expm")
     <environment: 0x000000000b6d27e0>
     Methods may be defined for arguments: x
     Use showMethods("expm") for currently available ones.
     > ## the dgeMatrix method:
     > expm.Matr.dge <- function(x) getDataPart(getMethod("expm", "dgeMatrix"))(Matrix::..2dge(x))
     > expmList <-
     + list(Matr = Matrix::expm,
     + Matr.d = expm.Matr.dge,
     + Ward = function(x) expm::expm(x, "Ward77"),
     + s.P.s = function(x) expm::expm(x, "Pade"),
     + s.P.sO= function(x) expm::expm(x, "PadeO"),
     + s.P.sRBS= function(x) expm::expm(x, "PadeRBS"),
     + sPs.H08.= function(x) expm:: expm.Higham08(x, balancing=FALSE),
     + sPs.H08b= function(x) expm:: expm.Higham08(x, balancing= TRUE),
     + AmHi09.06= function(x) expm:::expm.AlMoHi09(x, p = 6),
     + AmHi09.08= function(x) expm:::expm.AlMoHi09(x, p = 8),
     + AmHi09.10= function(x) expm:::expm.AlMoHi09(x, p = 10),
     + AmHi09.12= function(x) expm:::expm.AlMoHi09(x, p = 12),
     + AmHi09.13= function(x) expm:::expm.AlMoHi09(x, p = 13),
     + s.T.s = function(x) expm::expm(x, "Taylor"),
     + s.T.sO= function(x) expm::expm(x, "TaylorO"),
     + Eigen = expm.safe.Eigen,
     + hybrid= function(x) expm::expm(x, "hybrid")
     + )
     >
     >
     > set.seed(12)
     > fRE <- replicate(if(doExtras) 100 else 20,
     + re.facMat(20, expmList))
     > cat("Number of correct decimal digits for facMat(20, rnorm):\n")
     Number of correct decimal digits for facMat(20, rnorm):
     > summary(-log10(t(fRE["relErr",,])))
     Matr Matr.d Ward s.P.s s.P.sO
     Min. :13.7 Min. :13.7 Min. :13.7 Min. :12.8 Min. :12.8
     1st Qu.:14.2 1st Qu.:14.2 1st Qu.:14.2 1st Qu.:13.1 1st Qu.:13.1
     Median :14.3 Median :14.3 Median :14.3 Median :13.2 Median :13.2
     Mean :14.3 Mean :14.3 Mean :14.3 Mean :13.2 Mean :13.2
     3rd Qu.:14.5 3rd Qu.:14.5 3rd Qu.:14.5 3rd Qu.:13.3 3rd Qu.:13.3
     Max. :14.7 Max. :14.7 Max. :14.7 Max. :13.6 Max. :13.6
     s.P.sRBS sPs.H08. sPs.H08b AmHi09.06 AmHi09.08
     Min. :13.9 Min. :14.4 Min. :14.4 Min. : 6.61 Min. : 9.78
     1st Qu.:14.3 1st Qu.:14.8 1st Qu.:14.8 1st Qu.: 7.69 1st Qu.:11.27
     Median :14.4 Median :14.9 Median :14.9 Median : 8.24 Median :11.72
     Mean :14.4 Mean :14.9 Mean :14.9 Mean : 8.47 Mean :12.19
     3rd Qu.:14.6 3rd Qu.:15.0 3rd Qu.:15.0 3rd Qu.: 9.41 3rd Qu.:13.56
     Max. :14.8 Max. :15.2 Max. :15.2 Max. :10.35 Max. :14.40
     AmHi09.10 AmHi09.12 AmHi09.13 s.T.s s.T.sO
     Min. :13.4 Min. :14.5 Min. :14.5 Min. :13.0 Min. :13.0
     1st Qu.:14.6 1st Qu.:14.7 1st Qu.:14.7 1st Qu.:13.4 1st Qu.:13.4
     Median :14.7 Median :14.8 Median :14.8 Median :13.5 Median :13.5
     Mean :14.6 Mean :14.8 Mean :14.8 Mean :13.5 Mean :13.5
     3rd Qu.:14.8 3rd Qu.:14.9 3rd Qu.:14.9 3rd Qu.:13.7 3rd Qu.:13.7
     Max. :15.0 Max. :15.2 Max. :15.2 Max. :13.9 Max. :13.9
     Eigen hybrid
     Min. :14.2 Min. :14.1
     1st Qu.:14.3 1st Qu.:14.3
     Median :14.5 Median :14.5
     Mean :14.5 Mean :14.4
     3rd Qu.:14.6 3rd Qu.:14.6
     Max. :14.6 Max. :14.6
     >
     >
     > ## Now look at that:
     > boxplot(t(fRE["relErr",,]), log="y", notch=TRUE, ylim = c(8e-16, 1e-8),
     + main = "relative errors for 'random' eigen-ok 20 x 20 matrix")
     Warning message:
     In bxp(list(stats = c(1.8711533705969e-15, 2.97283678100251e-15, :
     some notches went outside hinges ('box'): maybe set notch=FALSE
     >
     > showProc.time()
     Time elapsed: 0.88 0 0.89
     >
     > if(doExtras) {
     + str(rf100 <- replicate(20, re.facMat(100, expmList)))
     + print(1000*t(apply(rf100["c.time",,], 1, summary)))
     + ## lynne {Linux 2.6.34.7-56.fc13.x86_64 --- AMD Phenom II X4 925}:
     + ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     + ## Ward 23 24 24.5 24.4 25.0 25
     + ## s.P.s 107 109 109.0 109.0 109.0 112
     + ## s.P.sO 188 190 191.0 192.0 193.0 198
     + ## s.P.sRBS 17 18 19.0 18.9 19.2 21
     + ## sPs.H08. 15 17 18.0 17.6 18.0 19
     + ## sPs.H08b 18 18 19.0 23.4 20.0 107
     + ## s.T.s 44 45 45.0 45.6 46.0 48
     + ## s.T.sO 96 98 99.0 100.0 100.0 116
     + ## Eigen 18 19 20.0 24.4 21.0 109
     + ## hybrid 40 42 42.0 47.1 44.0 133
     +
     + ##--> take out the real slow ones for the subsequent tests:
     + `%w/o%` <- function(x, y) x[!x %in% y] #-- x without y
     + print(nms.swift <- names(expmList) %w/o%
     + c("s.P.s", "s.P.sO", "s.T.s", "s.T.sO"))
     + expmL.swift <- expmList[nms.swift]
     +
     + set.seed(18) ## 12 replicates is too small .. but then it's too slow otherwise:
     + rf400 <- replicate(12, re.facMat(400, expmL.swift))
     + print(1000*t(apply(rf400["c.time",,], 1, summary)))
     + ## lynne:
     + ## Min. 1st Qu. Median Mean 3rd Qu. Max.
     + ## Ward 1740 1790 1830 1820 1860 1900
     + ## s.P.sRBS 1350 1420 1440 1430 1450 1460
     + ## sPs.H08. 1020 1030 1130 1140 1210 1290
     + ## sPs.H08b 1120 1130 1220 1220 1300 1390
     + ## Eigen 962 977 989 992 1000 1030
     + ## hybrid 2740 2800 2840 2840 2890 2910
     +
     + showProc.time()
     +
     + }## if(doExtras) only
     >
     > ## Now try an example with badly conditioned "random" M matrix...
     > ## ...
     > ## ... (not yet)
     >
     >
     > ### m2ex3() --- The 2x2 example with bad condition , see A3 in ./ex2.R
     >
     > RE <- re.m2ex3(1e-8, expmList)
     > sort(RE)# Ward + both sps.H08 are best; s.P.s fair, Eigen (and hybrid): ~1e-9
     Matr Matr.d Ward s.P.sRBS sPs.H08. sPs.H08b AmHi09.13 AmHi09.12 s.P.s
     5.584e-33 5.584e-33 5.584e-33 5.584e-33 5.584e-33 5.584e-33 1.509e-16 2.515e-16 2.113e-15
     s.P.sO s.T.s s.T.sO AmHi09.10 Eigen hybrid AmHi09.08 AmHi09.06
     2.113e-15 2.113e-15 2.113e-15 2.325e-12 1.836e-09 1.836e-09 4.819e-09 3.156e-06
     >
     > eps <- 10^-(1:18)
     > t.m2 <- t(sapply(eps, re.m2ex3, EXPMlist = expmList))
     Error in solve.default(V) :
     system is computationally singular: reciprocal condition number = 1.11022e-16
     Error in solve.default(V) :
     system is computationally singular: reciprocal condition number = 1.11022e-16
     Error in solve.default(V) :
     system is computationally singular: reciprocal condition number = 1.11022e-16
     > ## --> 3 error messages from solve(V), 5 error messages from try(. "R_Eigen" ...)
     >
     > showProc.time()
     Time elapsed: 0.37 0 0.38
     >
     > cbind(sort(apply(log(t.m2),2, median, na.rm=TRUE)))
     [,1]
     sPs.H08. -86.05
     sPs.H08b -86.05
     Matr -83.62
     Matr.d -83.62
     Ward -83.62
     s.P.sRBS -81.20
     AmHi09.13 -36.43
     AmHi09.12 -35.92
     s.P.s -33.79
     s.P.sO -33.79
     s.T.s -33.79
     s.T.sO -33.79
     AmHi09.10 -26.79
     hybrid -23.74
     Eigen -20.12
     AmHi09.08 -19.15
     AmHi09.06 -12.67
     > ## 'na.rm=TRUE' needed for Eigen which blows up for the last 3 eps
     > t.m2.ranks <- sort(rowMeans(apply(t.m2, 1, rank)))
     > cbind(signif(t.m2.ranks, 3))
     [,1]
     sPs.H08. 3.19
     sPs.H08b 3.19
     s.P.sRBS 3.92
     Matr 5.00
     Matr.d 5.00
     Ward 5.00
     AmHi09.13 5.67
     AmHi09.12 7.36
     s.T.s 10.20
     s.T.sO 10.20
     s.P.s 11.10
     s.P.sO 11.10
     hybrid 12.50
     AmHi09.10 13.70
     Eigen 14.40
     AmHi09.08 15.10
     AmHi09.06 16.40
     > ## lynne (x86_64, Linux 3.14.4-100; Intel i7-4765T), 2014-09:
     > ## sPs.H08. 2.67
     > ## sPs.H08b 2.67
     > ## s.P.sRBS 3.06
     > ## Ward 4.03
     > ## AmHi09.13 4.33 <<- still not close to H08 !
     > ## AmHi09.12 5.86
     > ## s.T.s 8.33
     > ## s.T.sO 8.33
     > ## s.P.s 9.11
     > ## s.P.sO 9.11
     > ## hybrid 10.80
     > ## AmHi09.10 11.70 << astonishingly bad
     > ## Eigen 12.60
     > ## AmHi09.08 13.10
     > ## AmHi09.06 14.40
     >
     > print(t.m2[, names(t.m2.ranks)[1:8]], digits = 3)
     sPs.H08. sPs.H08b s.P.sRBS Matr Matr.d Ward AmHi09.13 AmHi09.12
     [1,] 9.99e-17 9.99e-17 2.00e-16 6.02e-16 6.02e-16 6.02e-16 3.50e-16 5.52e-16
     [2,] 1.01e-16 1.01e-16 2.01e-16 5.53e-16 5.53e-16 5.53e-16 1.01e-16 2.51e-16
     [3,] 5.03e-17 5.03e-17 2.51e-16 3.02e-16 3.02e-16 3.02e-16 2.51e-16 3.02e-16
     [4,] 2.51e-16 2.51e-16 3.75e-25 1.51e-16 1.51e-16 1.51e-16 1.51e-16 2.51e-16
     [5,] 3.52e-16 3.52e-16 2.01e-16 1.51e-16 1.51e-16 1.51e-16 1.01e-16 1.01e-16
     [6,] 2.01e-16 2.01e-16 3.52e-16 3.52e-16 3.52e-16 3.52e-16 2.01e-16 3.52e-16
     [7,] 3.52e-16 3.52e-16 3.52e-16 6.04e-16 6.04e-16 6.04e-16 7.15e-31 2.01e-16
     [8,] 5.58e-33 5.58e-33 5.58e-33 5.58e-33 5.58e-33 5.58e-33 1.51e-16 2.51e-16
     [9,] 4.36e-35 4.36e-35 4.36e-35 4.36e-35 4.36e-35 4.36e-35 1.51e-16 2.51e-16
     [10,] 0.00e+00 0.00e+00 6.82e-37 0.00e+00 0.00e+00 0.00e+00 1.51e-16 2.51e-16
     [11,] 0.00e+00 0.00e+00 5.33e-39 5.33e-39 5.33e-39 5.33e-39 1.51e-16 2.51e-16
     [12,] 4.16e-41 4.16e-41 0.00e+00 4.16e-41 4.16e-41 4.16e-41 1.51e-16 2.51e-16
     [13,] 0.00e+00 0.00e+00 6.50e-43 6.50e-43 6.50e-43 6.50e-43 1.51e-16 2.51e-16
     [14,] 0.00e+00 0.00e+00 5.08e-45 5.08e-45 5.08e-45 5.08e-45 1.51e-16 2.51e-16
     [15,] 3.97e-47 3.97e-47 7.94e-47 7.94e-47 7.94e-47 7.94e-47 1.51e-16 2.51e-16
     [16,] 0.00e+00 0.00e+00 6.20e-49 6.20e-49 6.20e-49 6.20e-49 1.51e-16 2.51e-16
     [17,] 0.00e+00 0.00e+00 4.84e-51 0.00e+00 0.00e+00 0.00e+00 1.51e-16 2.51e-16
     [18,] 3.78e-53 3.78e-53 0.00e+00 3.78e-53 3.78e-53 3.78e-53 1.51e-16 2.51e-16
     > ## ==> 1st class: H08 (both) and (but slightly better than) Ward
     > ## 2nd class s.T.s and s.P.s
     > ## "bad" : hybrid and Eigen
     >
     > ## ??? AmHi09 - methods, up to order = 10 are worse !
     > if(require(RColorBrewer)) {
     + ## Bcol <- brewer.pal(ncol(t.m2),"Dark2")
     + Bcol <- brewer.pal(min(9, ncol(t.m2)), "Set1")
     + Bcol <- Bcol[sqrt(colSums(col2rgb(Bcol)^2)) < 340]
     + ## FIXME: more colors ==> ~/R/MM/GRAPHICS/color-palettes.R
     + } else {
     + ## 7 from Dark2
     + ## Bcol <- c("#1B9E77", "#D95F02", "#7570B3", "#E7298A",
     + ## "#66A61E", "#E6AB02", "#A6761D")
     + ## Rather: those from "Set1"
     + Bcol <- c("#E41A1C", "#377EB8", "#4DAF4A",
     + "#984EA3", "#FF7F00", # too bright: "#FFFF33",
     + "#A65628", "#F781BF", "#999999")
     + }
     Loading required package: RColorBrewer
     >
     > matplot(eps, t.m2, type = "b", log = "xy", col=Bcol, lty = 1:9, pch=1:15,
     + axes=FALSE, frame = TRUE,
     + xlab = expression(epsilon), ylab = "relative error",
     + main = expression(expm(A, method == "*") *" relative errors for " *
     + A == bgroup("[", atop({-1} *" "* 1, {epsilon^2} *" "*{-1}), "]")))
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     21 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     2 y values <= 0 omitted from logarithmic plot
     > legend("bottomright",colnames(t.m2), col=Bcol, lty = 1:9, pch=1:15,
     + inset = 0.02)
     > if(require("sfsmisc")) {
     + sfsmisc::eaxis(1, labels=FALSE)
     + sfsmisc::eaxis(1, at = eps[c(TRUE,FALSE)])
     + sfsmisc::eaxis(2)
     + ## sfsmisc::eaxis(2, labels=FALSE)
     + ## op <- par(las=2)
     + ## sfsmisc::eaxis(2, at = axTicks(2,log=TRUE)[c(TRUE,FALSE,FALSE)])
     + ## par(op)
     + } else {
     + axis(1)
     + axis(2)
     + }
     Loading required package: sfsmisc
     >
     > ## typical case:
     > ep <- 1e-10
     > (me <- m2ex3(ep))
     $A
     [,1] [,2]
     [1,] -1e+00 1
     [2,] 1e-20 -1
    
     $expA
     [,1] [,2]
     [1,] 3.679e-01 0.3679
     [2,] 3.679e-21 0.3679
    
     > me$expA * exp(1) ## the correct value ; numerically identical to simple matrix:
     [,1] [,2]
     [1,] 1e+00 1
     [2,] 1e-20 1
     > ## identical() not fulfilled e.g. on Solaris
     > stopifnot(all.equal(me$expA * exp(1),
     + rbind(c( 1, 1),
     + c(ep^2, 1)),
     + tolerance = 1e-14))
     > ## The relative error (matrices):
     > lapply(expmList, function(EXPM) 1 - EXPM(me$A)/me$expA)
     $Matr
     2 x 2 Matrix of class "dgeMatrix"
     [,1] [,2]
     [1,] 0 0
     [2,] 0 0
    
     $Matr.d
     2 x 2 Matrix of class "dgeMatrix"
     [,1] [,2]
     [1,] 0 0
     [2,] 0 0
    
     $Ward
     [,1] [,2]
     [1,] 0 0
     [2,] 0 0
    
     $s.P.s
     [,1] [,2]
     [1,] -2.22e-15 -2.22e-15
     [2,] -2.22e-15 -2.22e-15
     attr(,"accuracy")
     [1] 1.665e-16
    
     $s.P.sO
     [,1] [,2]
     [1,] -2.22e-15 -2.22e-15
     [2,] -2.22e-15 -2.22e-15
     attr(,"accuracy")
     [1] 1.665e-16
    
     $s.P.sRBS
     [,1] [,2]
     [1,] 0.00e+00 0
     [2,] -2.22e-16 0
    
     $sPs.H08.
     [,1] [,2]
     [1,] 0 0
     [2,] 0 0
    
     $sPs.H08b
     [,1] [,2]
     [1,] 0 0
     [2,] 0 0
    
     $AmHi09.06
     [,1] [,2]
     [1,] -1.153e-06 7.160e-06
     [2,] 7.160e-06 -1.153e-06
    
     $AmHi09.08
     [,1] [,2]
     [1,] -1.414e-09 1.163e-08
     [2,] 1.163e-08 -1.414e-09
    
     $AmHi09.10
     [,1] [,2]
     [1,] -5.702e-13 5.836e-12
     [2,] 5.836e-12 -5.702e-13
    
     $AmHi09.12
     [,1] [,2]
     [1,] 1.110e-16 4.441e-16
     [2,] 7.772e-16 1.110e-16
    
     $AmHi09.13
     [,1] [,2]
     [1,] 1.11e-16 1.11e-16
     [2,] 2.22e-16 1.11e-16
    
     $s.T.s
     [,1] [,2]
     [1,] -2.22e-15 -2.22e-15
     [2,] -2.22e-15 -2.22e-15
     attr(,"accuracy")
     [1] 0
    
     $s.T.sO
     [,1] [,2]
     [1,] -2.22e-15 -2.22e-15
     [2,] -2.22e-15 -2.22e-15
     attr(,"accuracy")
     [1] 0
    
     $Eigen
     [,1] [,2]
     [1,] 0.00e+00 -5.109e-07
     [2,] -1.12e-06 0.000e+00
    
     $hybrid
     [,1] [,2]
     [1,] 0.000e+00 -5.109e-07
     [2,] -2.419e-07 -2.220e-16
    
     >
     > ## Average number of correct digits [less "extreme" than plot above]
     > nDig <- sapply(expmList, function(EXPM) -log10(mean(abs(1 - EXPM(me$A)/me$expA))))
     > round(nDig, 2)
     Matr Matr.d Ward s.P.s s.P.sO s.P.sRBS sPs.H08. sPs.H08b AmHi09.06
     Inf Inf Inf 14.65 14.65 16.26 Inf Inf 5.38
     AmHi09.08 AmHi09.10 AmHi09.12 AmHi09.13 s.T.s s.T.sO Eigen hybrid
     8.19 11.49 15.44 15.86 14.65 14.65 6.39 6.73
     > ## Ward s.P.s s.P.sO s.T.s s.T.sO Eigen hybrid
     > ## 16.26 14.65 14.65 14.65 14.65 6.20 6.39 [AMD Opteron 64-bit]
     > ## Inf 14.65 14.65 14.65 14.65 6.74 6.33 [Pentium-M (32-bit)]
     >
     > ###--- rnilMat() : random upper triangular (zero-diagonal) nilpotent n x n matrix
     >
     > set.seed(17)
     > m <- rnilMat(10)
     > (m. <- as(m,"sparseMatrix"))# for nicer printing - and more *below*
     10 x 10 sparse Matrix of class "dgCMatrix"
     Error in isFALSE(suppRows) : could not find function "isFALSE"
     Calls: <Anonymous> -> <Anonymous> -> printSpMatrix2
     Execution halted
Flavor: r-oldrel-windows-ix86+x86_64