Running 'Frechet-test.R' [2s]
  Running 'bal-ex.R' [3s]
  Running 'ex.R' [2s]
  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.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.48 0 1.48 
  > 
  > ## 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.64 0 0.64 
  > 
  > ###--- 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.04 0 0.05 
  > 
  > 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: 0x000000000b0d1940>
  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.77 0.02 0.78 
  > 
  > 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.36 0 0.36 
  > 
  > 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