R Under development (unstable) (2024-08-17 r87027 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(expm)
Loading required package: Matrix

Attaching package: 'expm'

The following object is masked from 'package:Matrix':

    expm

> source(system.file("test-tools.R", package= "expm"), keep.source=FALSE)## -> assertError()...
Loading required package: tools
> 
> ## A matrix with 'Inf'
> mI <- rbind(0, c(-Inf, Inf, 0, 0), 0, 0)
> bal3 <-
+     list(dB = dgebal(mI, "B"), # = default
+          dP = dgebal(mI, "P"),
+          dN = dgebal(mI, "N"))
Warning messages:
1: In dgebal(mI, "B") : 'dgebal' is deprecated.
Use 'balance' instead.
See help("Deprecated")
2: In dgebal(mI, "P") : 'dgebal' is deprecated.
Use 'balance' instead.
See help("Deprecated")
3: In dgebal(mI, "N") : 'dgebal' is deprecated.
Use 'balance' instead.
See help("Deprecated")
> str(bal3)
List of 3
 $ dB:List of 4
  ..$ z    : num [1:4, 1:4] Inf 0 0 0 -Inf ...
  ..$ scale: num [1:4] 1 1 3 4
  ..$ i1   : int 1
  ..$ i2   : int 1
 $ dP:List of 4
  ..$ z    : num [1:4, 1:4] Inf 0 0 0 -Inf ...
  ..$ scale: num [1:4] 1 1 3 4
  ..$ i1   : int 1
  ..$ i2   : int 1
 $ dN:List of 4
  ..$ z    : num [1:4, 1:4] 0 -Inf 0 0 0 ...
  ..$ scale: num [1:4] 1 1 1 1
  ..$ i1   : int 1
  ..$ i2   : int 4
> stopifnot(identical(mI, bal3$dN$z),
+           with(bal3, all.equal(dB, dP, tol=1e-14)),
+           all.equal(bal3$dB$z, rbind(c(Inf,-Inf,0,0), 0,0,0), tol=1e-14),
+           all.equal(bal3$dB$scale, c(1,1,3,4)))
> assertError(dgebal(mI, "S"), verbose=TRUE)# gave infinite loop
Asserted error: R_dgebal(*, type="S"): Infinite matrix entry
> 
> 
> 
> ## Compare the two different "balance" pre-conditioning versions in Ward77:
> set.seed(1)
> mList <- lapply(integer(100), function(...) rSpMatrix(20, nnz=80))
> re20 <- sapply(mList, function(M)
+                relErr(expm(M, precond = "2bal"),
+                       expm(M, precond = "1bal")))
> re20 ## ahh.: zero or ~ 1e-13 ... good
  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> table(re20 == 0)

TRUE 
 100 
> summary(re20[re20 != 0])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
                                                
> ## Pentium M (ubuntu)
> ##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
> ## 2.593e-14 8.703e-14 1.282e-13 2.434e-13 4.177e-13 6.295e-13
> 
> demo(balanceTst) #-> the function definition and the first few examples


	demo(balanceTst)
	---- ~~~~~~~~~~

> balanceTst <- function(A) {
+ 
+     ## Purpose: Consistency checking of	 balance() {was "dgebal()"}
+     ## ----------------------------------------------------------------------
+     ## Arguments: a square matrix
+     ## ----------------------------------------------------------------------
+     ## Author: Martin Maechler, 20 Feb 2008 and on
+ 
+     n <- dim(A)[1]
+     ## do *the* three calls and look at result
+     P <- balance(A, "P")
+ 
+     doPerm <- function(A, pp, i1, i2) {
+         stopifnot(length(pp) == n, dim(A) == c(n,n),
+                   1 <= i1, i1 <= i2, i2 <= n)
+         A. <- A
+         if(i2 < n) { ## The upper part
+             for(i in n:(i2+1)) {    # 'p2' in *reverse* order
+                 ## swap	 i <-> pp[i]   both rows and columns
+                 tt <- A.[,i]; A.[,i] <- A.[,pp[i]]; A.[,pp[i]] <- tt
+                 tt <- A.[i,]; A.[i,] <- A.[pp[i],]; A.[pp[i],] <- tt
+             }
+         }
+         if(i1 > 1) { ## The lower part
+             for(i in 1:(i1-1)) {    # 'p1' in *forward* order
+                 tt <- A.[,i]; A.[,i] <- A.[,pp[i]]; A.[,pp[i]] <- tt
+                 tt <- A.[i,]; A.[i,] <- A.[pp[i],]; A.[pp[i],] <- tt
+             }
+         }
+         A.
+     }
+ 
+     checkPerm <- function(P, orig.A) {
+         didPerm <- ((leftP  <- (i1 <- P$i1) != 1L) |
+                     (rightP <- (i2 <- P$i2) != n))
+         if(didPerm) { ## *had* permutation -- now check my idea about it
+             pp <- as.integer(P$scale)
+             ## Permute A to become P$z :
+             A. <- doPerm(orig.A, pp = pp, i1=i1, i2=i2)
+             stopifnot(isTRUE(all.equal(A., P$z, tolerance = 1e-15)))
+ 
+             ## Now the reverse: Use pp[] and permute  A. "back to A":
+             if(leftP) { ## The lower part
+                 for(i in (i1-1):1) {    # 'p1' in *reverse* order
+                     tt <- A.[,i]; A.[,i] <- A.[,pp[i]]; A.[,pp[i]] <- tt
+                     tt <- A.[i,]; A.[i,] <- A.[pp[i],]; A.[pp[i],] <- tt
+                 }
+             }
+             if(rightP) { ## The upper part
+                 for(i in (i2+1):n) {    # 'p2' in *forward* order
+                     ## swap	 i <-> pp[i]   both rows and columns
+                     tt <- A.[,i]; A.[,i] <- A.[,pp[i]]; A.[,pp[i]] <- tt
+                     tt <- A.[i,]; A.[i,] <- A.[pp[i],]; A.[pp[i],] <- tt
+                 }
+             }
+             stopifnot(isTRUE(all.equal(A., orig.A, tolerance = 1e-15)))
+         }
+     }
+     checkPerm(P, orig.A = A)
+ 
+     S <- balance(P$z, "S")# "S" starting from result of "P"
+     stopifnot(S$i1 == 1, S$i2 == n)
+ 
+     ## Now check the scaling
+     checkScal <- function (d, A1, A2) {
+         stopifnot(length(d) == n, dim(A1) == dim(A2), dim(A2) == c(n,n))
+ 
+         ## A.scaled <- diag(1/d, n) \%*\% A1 \%*\% diag(d, n)
+         ## more efficiently:
+         A.scaled <- A1 * (rep(d, each = n) / d)
+         stopifnot(isTRUE(all.equal(A2, A.scaled, tolerance = 1e-15)))
+         ## Check the reverse:
+         S.rescaled <- A2 * (d * rep(1/d, each = n))
+         stopifnot(isTRUE(all.equal(A1, S.rescaled, tolerance = 1e-15)))
+     }
+     checkScal(d = S$scale, A1 = P$z, A2 = S$z)
+ 
+     B <- balance(A, "B")# "B" : B[oth]
+     stopifnot(P$i1 == B$i1, P$i2 == B$i2)
+     ## now check *both* permutation and scaling
+ 
+     A.perm <- doPerm(A, pp = as.integer(B$scale), i1=B$i1, i2=B$i2)
+     ## checkPerm(B, orig.A = A)
+ 
+     dB <- B$scale
+     dB[c(if(B$i1 > 1) 1:(B$i1-1),
+          if(B$i2 < n) (B$i2+1):n)] <- 1
+     checkScal(d = dB, A1 = A.perm, A2 = B$z)
+ 
+     ## return
+     list(P = P, S = S, B = B, Sz.eq.Bz = isTRUE(all.equal(S$z, B$z)))
+ }

> m4. <- rbind(c(-1,-2, 0, 0),
+              c( 0, 0,10,11),
+              c( 0, 0,12, 0),
+              c( 0,13, 0, 0))

> op <- options(str = strOptions(vec.len = 12))

> str(b4. <- balanceTst(m4.))
List of 4
 $ P       :List of 4
  ..$ z    : num [1:4, 1:4] -1 0 0 0 -2 0 13 0 0 11 0 0 0 10 0 12
  ..$ scale: num [1:4] 1 1 1 3
  ..$ i1   : int 2
  ..$ i2   : int 3
 $ S       :List of 4
  ..$ z    : num [1:4, 1:4] -1 0 0 0 -2 0 13 0 0 11 0 0 0 10 0 12
  ..$ scale: num [1:4] 1 1 1 1
  ..$ i1   : int 1
  ..$ i2   : int 4
 $ B       :List of 4
  ..$ z    : num [1:4, 1:4] -1 0 0 0 -2 0 13 0 0 11 0 0 0 10 0 12
  ..$ scale: num [1:4] 1 1 1 3
  ..$ i1   : int 2
  ..$ i2   : int 3
 $ Sz.eq.Bz: logi TRUE

> with(b4., all.equal(P, B)) # TRUE (everywhere?)
[1] TRUE

> ## better (?) example
> (m <- matrix(c(0,-1,0,-2,10, rep(0,11)), 4,4))
     [,1] [,2] [,3] [,4]
[1,]    0   10    0    0
[2,]   -1    0    0    0
[3,]    0    0    0    0
[4,]   -2    0    0    0

> str(ba <- balanceTst(m))
List of 4
 $ P       :List of 4
  ..$ z    : num [1:4, 1:4] 0 0 0 0 0 0 10 0 -2 -1 0 0 0 0 0 0
  ..$ scale: num [1:4] 3 1 1 3
  ..$ i1   : int 2
  ..$ i2   : int 3
 $ S       :List of 4
  ..$ z    : num [1:4, 1:4] 0 0 0 0 0 0 2.5 0 -2 -4 0 0 0 0 0 0
  ..$ scale: num [1:4] 1 0.25 1 1
  ..$ i1   : int 1
  ..$ i2   : int 4
 $ B       :List of 4
  ..$ z    : num [1:4, 1:4] 0 0 0 0 0 0 2.5 0 -2 -4 0 0 0 0 0 0
  ..$ scale: num [1:4] 3 0.25 1 3
  ..$ i1   : int 2
  ..$ i2   : int 3
 $ Sz.eq.Bz: logi TRUE

> (eq <- with(ba, all.equal(S$z, B$z))) # TRUE now (everywhere?)
[1] TRUE

> ba$Sz.eq.Bz # ditto
[1] TRUE

> ## a non-empty ``less-balanced'' example  ---
> 
> m4 <- matrix(outer(2^(0:7),c(-1,1)), 4,4)

> m4[lower.tri(m4)] <- 0 #--> upper triangular ==> will have many permutations

> ## now permute it; so balance() will find the permutation
> p <- c(4,2:1,3); m4 <- m4[p,p]

> m4
     [,1] [,2] [,3] [,4]
[1,]  128    0    0    0
[2,]   32  -32    0    2
[3,]   16  -16   -1    1
[4,]   64    0    0    4

> str(dm4 <- balanceTst(m4)) # much permutation!  i1 = i2 = 1 !
List of 4
 $ P       :List of 4
  ..$ z    : num [1:4, 1:4] -1 0 0 0 -16 -32 0 0 1 2 4 0 16 32 64 128
  ..$ scale: num [1:4] 1 2 1 1
  ..$ i1   : int 1
  ..$ i2   : int 1
 $ S       :List of 4
  ..$ z    : num [1:4, 1:4] -1 0 0 0 -1 -32 0 0 0.25 8 4 0 1 32 16 128
  ..$ scale: num [1:4] 16 1 4 1
  ..$ i1   : int 1
  ..$ i2   : int 4
 $ B       :List of 4
  ..$ z    : num [1:4, 1:4] -1 0 0 0 -16 -32 0 0 1 2 4 0 16 32 64 128
  ..$ scale: num [1:4] 1 2 1 1
  ..$ i1   : int 1
  ..$ i2   : int 1
 $ Sz.eq.Bz: logi FALSE

> ##----------- Complex examples
> zba4 <- balanceTst(m4 + 3i * m4)

> str(zba4)
List of 4
 $ P       :List of 4
  ..$ z    : cplx [1:4, 1:4] -1-3i 0+0i 0+0i 0+0i -16-48i -32-96i 0+0i 0+0i 1+3i ...
  ..$ scale: num [1:4] 1 2 1 1
  ..$ i1   : int 1
  ..$ i2   : int 1
 $ S       :List of 4
  ..$ z    : cplx [1:4, 1:4] -1-3i 0+0i 0+0i 0+0i -1-3i -32-96i 0+0i 0+0i 0.25+0.75i ...
  ..$ scale: num [1:4] 16 1 4 1
  ..$ i1   : int 1
  ..$ i2   : int 4
 $ B       :List of 4
  ..$ z    : cplx [1:4, 1:4] -1-3i 0+0i 0+0i 0+0i -16-48i -32-96i 0+0i 0+0i 1+3i ...
  ..$ scale: num [1:4] 1 2 1 1
  ..$ i1   : int 1
  ..$ i2   : int 1
 $ Sz.eq.Bz: logi FALSE

> zba <-  balanceTst(m*(1 + 1i))

> str(zba)
List of 4
 $ P       :List of 4
  ..$ z    : cplx [1:4, 1:4] 0+0i 0+0i 0+0i 0+0i 0+0i 0+0i 10+10i 0+0i -2-2i ...
  ..$ scale: num [1:4] 3 1 1 3
  ..$ i1   : int 2
  ..$ i2   : int 3
 $ S       :List of 4
  ..$ z    : cplx [1:4, 1:4] 0+0i 0+0i 0+0i 0+0i 0+0i 0+0i 2.5+2.5i 0+0i -2-2i ...
  ..$ scale: num [1:4] 1 0.25 1 1
  ..$ i1   : int 1
  ..$ i2   : int 4
 $ B       :List of 4
  ..$ z    : cplx [1:4, 1:4] 0+0i 0+0i 0+0i 0+0i 0+0i 0+0i 2.5+2.5i 0+0i -2-2i ...
  ..$ scale: num [1:4] 3 0.25 1 3
  ..$ i1   : int 2
  ..$ i2   : int 3
 $ Sz.eq.Bz: logi TRUE

> stopifnot(exprs = {
+     all.equal(ba$ S$z, Re(zba$ S$z))
+     all.equal(ba$ S$z, Im(zba$ S$z))
+     all.equal(dm4$ S$z,     Re(zba4$ S$z))
+     all.equal(dm4$ S$z * 3, Im(zba4$ S$z))
+ })

> options(op) # revert
> ##   ==========
> 
> dm4. <- dgebal(m4)
Warning message:
In dgebal(m4) : 'dgebal' is deprecated.
Use 'balance' instead.
See help("Deprecated")
> storage.mode(m4) <- "integer"
> stopifnot(identical(dm4., dgebal(m4)))
Warning message:
In dgebal(m4) : 'dgebal' is deprecated.
Use 'balance' instead.
See help("Deprecated")
> 
> expm(m)
             [,1]        [,2] [,3] [,4]
[1,] -0.999786073 -0.06540707    0    0
[2,]  0.006540707 -0.99978607    0    0
[3,]  0.000000000  0.00000000    1    0
[4,]  0.013081414 -3.99957215    0    1
> expm(m,"Pade") ## are different indeed {when bug still existed}
             [,1]        [,2] [,3] [,4]
[1,] -0.999786073 -0.06540707    0    0
[2,]  0.006540707 -0.99978607    0    0
[3,]  0.000000000  0.00000000    1    0
[4,]  0.013081414 -3.99957215    0    1
attr(,"accuracy")
[1] 7.918807e-14
> expm(m,"R_Pade")# same as Pade
             [,1]        [,2] [,3] [,4]
[1,] -0.999786073 -0.06540707    0    0
[2,]  0.006540707 -0.99978607    0    0
[3,]  0.000000000  0.00000000    1    0
[4,]  0.013081414 -3.99957215    0    1
> 
> 
> ## a non-empty ``non-balanced'' example  ---
> 
> expm.t.identity(m4, "Ward")
[1] TRUE
> 
> m6 <- zeroTrace(matrix(outer(2^(-8:9),c(-1,1)), 6,6)); m6
            [,1]      [,2]      [,3]         [,4]      [,5]     [,6]
[1,] -75.2584635  -0.25000  -16.0000   0.00390625   0.25000  16.0000
[2,]  -0.0078125 -75.75456  -32.0000   0.00781250   0.50000  32.0000
[3,]  -0.0156250  -1.00000 -139.2546   0.01562500   1.00000  64.0000
[4,]  -0.0312500  -2.00000 -128.0000 -75.22330729   2.00000 128.0000
[5,]  -0.0625000  -4.00000 -256.0000   0.06250000 -71.25456 256.0000
[6,]  -0.1250000  -8.00000 -512.0000   0.12500000   8.00000 436.7454
> m6[lower.tri(m6)] <- 0 ## plus one non-zero
> m6[4,2] <- 77
> p <- c(6,4,5,2:1,3); m6 <- m6[p,p]
> expm.t.identity(m6, "Ward") ##  difference; indeed
[1] TRUE
> expm(m6) # is very different from
              [,1]         [,2]         [,3]          [,4]         [,5]
[1,] 4.743903e+189 0.000000e+00 0.000000e+00  0.000000e+00 0.000000e+00
[2,] 1.235301e+189 2.158311e-33 7.038631e-32  1.285266e-31 0.000000e+00
[3,] 2.390628e+189 0.000000e+00 1.133807e-31  0.000000e+00 0.000000e+00
[4,] 2.653829e+188 1.627057e-37 7.493870e-34  1.271561e-33 0.000000e+00
[5,] 1.326905e+188 8.096068e-38 3.743330e-34 -4.016897e-34 2.068543e-33
[6,] 5.312842e+188 5.266261e-37 1.683517e-33  3.098706e-35 0.000000e+00
              [,6]
[1,]  0.000000e+00
[2,] -6.346149e-32
[3,]  0.000000e+00
[4,] -6.406872e-34
[5,] -3.187638e-34
[6,] -1.529642e-35
> expm(m6,"R_Pade")
              [,1]         [,2]         [,3]          [,4]         [,5]
[1,] 4.743903e+189 0.000000e+00 0.000000e+00  0.000000e+00 0.000000e+00
[2,] 1.235301e+189 2.158311e-33 7.038631e-32  1.285266e-31 0.000000e+00
[3,] 2.390628e+189 0.000000e+00 1.133807e-31  0.000000e+00 0.000000e+00
[4,] 2.653829e+188 1.627057e-37 7.493870e-34  1.271561e-33 0.000000e+00
[5,] 1.326905e+188 8.096068e-38 3.743330e-34 -4.016897e-34 2.068543e-33
[6,] 5.312842e+188 5.266261e-37 1.683517e-33  3.098706e-35 0.000000e+00
              [,6]
[1,]  0.000000e+00
[2,] -6.346149e-32
[3,]  0.000000e+00
[4,] -6.406872e-34
[5,] -3.187638e-34
[6,] -1.529642e-35
> 
> str(dm6 <- balanceTst(m6))
List of 4
 $ P       :List of 4
  ..$ z    : num [1:6, 1:6] -75.3 0 0 0 0 ...
  ..$ scale: num [1:6] 3 1 1 1 3 1
  ..$ i1   : int 2
  ..$ i2   : int 4
 $ S       :List of 4
  ..$ z    : num [1:6, 1:6] -75.3 0 0 0 0 ...
  ..$ scale: num [1:6] 1 1 1 1 2 1
  ..$ i1   : int 1
  ..$ i2   : int 6
 $ B       :List of 4
  ..$ z    : num [1:6, 1:6] -75.3 0 0 0 0 ...
  ..$ scale: num [1:6] 3 1 1 1 3 1
  ..$ i1   : int 2
  ..$ i2   : int 4
 $ Sz.eq.Bz: logi FALSE
> ## Now, that's interesting:
> ##
> ## 1.  'S' scales *more* (2 .. 5) than just (2:4 == i1:i2) !
> ##
> ## 2.  'B' has quite different scaling and it does (must!) obey rule
> ##     scale i1:i2 only
> ##
> ## 3. 'B'(oth) is better than  "P" and "S" separately:
> ##
> kappa(eigen(m6)$vectors)#      597.5588
[1] 597.5588
> kappa(eigen(dm6$P$z)$vectors)# 597.5588
[1] 597.5588
> kappa(eigen(dm6$S$z)$vectors)#  42.58396
[1] 621.5746
> kappa(eigen(dm6$B$z)$vectors)#  22.20266
[1] 597.5588
> 
> 
> ## An n=17 example where octave's expm() is wrong too
> m17 <- matrix(c(10,0, 0, 2, 3,-1, 0, 0, 0, 0, 0, 4, 0, 5, 0, 0,-2,
+                 0, 0, 0, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 7, 0,
+                 0, 0,10, 0, 0,-4, 9, 0, 0, 0,-5, 0,-6, 0, 0, 0, 0,
+                 0, 0,-7, 0, 0, 0, 0, 0, 0,10, 0, 0, 0, 0, 0,11, 0,
+                 0, 0, 0, 0, 0, 0,12, 0, 0, 0, 0, 0,-8, 0, 0, 0, 0,
+                 0, 0,-9, 0, 0, 0, 0, 0, 0,-10,0,13,14,-11,-12,-13, 0,
+                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                 0, 0, 0, 0, 0, 0, 0,10, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                 0, 0, 0, 0, 0,-14,16,0,-10,0,17, 0, 0, 0, 0, 0, 0,
+                 0, 0,-16,0, 0,18,19, 0, 0, 0, 0, 0, 0, 0,20, 0, 21,
+                 22,0, 0, 0, 0, 0,-17,0, 0, 0,-10,-19,-20,0,0,0, 0,
+                 0,-21,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                 0, 0, 0,23,24, 0,25,26, 0, 0,27,-22,0,28,-23,0,-24,
+                 0,-25,0,29, 0, 0, 0, 0, 0, 0, 0,30,31, 0, 0, 0, 0,
+                 0, 0,-26,32,0, 0, 0, 0, 0,-27,0,33,34, 0, 0, 0, 0,
+                 0,-28,-29,0,0, 0,35, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                 0, 0, 0, 0, 0, 0,36,37, 0, 0, 0, 0, 0, 0, 0, 0,-10),
+               17, 17)
> str(dm17 <- balanceTst(m17))
List of 4
 $ P       :List of 4
  ..$ z    : num [1:17, 1:17] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ scale: num [1:17] 7 8 9 1 1 1 1 1 1 1 ...
  ..$ i1   : int 4
  ..$ i2   : int 16
 $ S       :List of 4
  ..$ z    : num [1:17, 1:17] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ scale: num [1:17] 1 4 1 2 1 1 1 2 2 1 ...
  ..$ i1   : int 1
  ..$ i2   : int 17
 $ B       :List of 4
  ..$ z    : num [1:17, 1:17] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ scale: num [1:17] 7 8 9 2 2 1 1 2 2 1 ...
  ..$ i1   : int 4
  ..$ i2   : int 16
 $ Sz.eq.Bz: logi FALSE
> sapply(dm17[1:3], `[[`, "scale")
      P   S B
 [1,] 7 1.0 7
 [2,] 8 4.0 8
 [3,] 9 1.0 9
 [4,] 1 2.0 2
 [5,] 1 1.0 2
 [6,] 1 1.0 1
 [7,] 1 1.0 1
 [8,] 1 2.0 2
 [9,] 1 2.0 2
[10,] 1 1.0 1
[11,] 1 1.0 1
[12,] 1 2.0 2
[13,] 1 1.0 1
[14,] 1 1.0 1
[15,] 1 1.0 1
[16,] 1 1.0 1
[17,] 9 0.5 9
> 
> ## The balancing was really rather harmful -- cond(V) *not* improved:
> condX <- function(x) kappa(x, exact=TRUE)
> condX(eigen(m17)$vectors)#      8.9e16
[1] 1.087882e+17
> condX(eigen(dm17$P$z)$vectors)# 1.37e17
[1] 1.36865e+17
> condX(eigen(dm17$S$z)$vectors)# 1.44e17
[1] 6.451363e+16
> condX(eigen(dm17$B$z)$vectors)# 1.43e17 (very slightly smaller)
[1] 1.437207e+17
> 
> proc.time()
   user  system elapsed 
   1.70    0.23    1.95