R Under development (unstable) (2023-11-29 r85646 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 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.

> #### Matrix Factorizations  --- of all kinds
> 
> ## for R_DEFAULT_PACKAGES=NULL :
> library(stats)
> library(utils)
> 
> library(Matrix)
> source(system.file("test-tools.R", package = "Matrix"))# identical3() etc
Loading required package: tools
> options(warn = 0)
> is64bit <- .Machine$sizeof.pointer == 8
> cat("doExtras:", doExtras,";  is64bit:", is64bit, "\n")
doExtras: FALSE ;  is64bit: TRUE 
> 
> ### "sparseQR" : Check consistency of methods
> ##   --------
> data(KNex, package = "Matrix")
> mm <- KNex$mm
> y  <- KNex$y
> stopifnot(is((Y <- Matrix(y)), "dgeMatrix"))
> md <- as(mm, "matrix")                  # dense
> 
> (cS <- system.time(Sq <- qr(mm))) # 0.009
   user  system elapsed 
   0.02    0.00    0.02 
> (cD <- system.time(Dq <- qr(md))) # 0.499 (lynne, 2014 f); 1.04 lynne 2019 ?????
   user  system elapsed 
   0.48    0.02    0.50 
> cD[1] / cS[1] # dense is  much ( ~ 100--170 times) slower
user.self 
       24 
> 
> ## chkQR() in ../inst/test-tools-1.R ;
> 
> if(doExtras) { ## ~ 20 sec {"large" example}   + 2x qr.R() warnings
+     cat("chkQR( <KNex> ) .. takes time .. ")
+     system.time(chkQR(mm, y=y, a.qr = Sq, verbose=TRUE))
+     system.time(chkQR(md, y=y, a.qr = Dq, verbose=TRUE))
+     cat(" done: [Ok]\n")
+ }
> 
> ## consistency of results dense and sparse
> ##	chk.qr.D.S() and  checkQR.DS.both() >>> ../inst/test-tools-Matrix.R
> chk.qr.D.S(Dq, Sq, y, Y)
> 
> ## Another small example with pivoting (and column name "mess"):
> suppressWarnings(RNGversion("3.5.0")); set.seed(1)
> X <- rsparsematrix(9,5, 1/4, dimnames=list(paste0("r", 1:9), LETTERS[1:5]))
> qX <- qr(X); qd <- qr(as(X, "matrix"))
> ## are the same (now, *including* names):
> assert.EQ(print(qr.coef(qX, 1:9)), qr.coef(qd, 1:9), tol=1e-14)
          A           B           C           D           E 
  1.2056600  -6.1033750 -20.4434698 -10.1046230  -0.9090909 
> chk.qr.D.S(d. = qd, s. = qX, y = 1:9)
> 
> 
> 
> ## rank deficient QR cases: ---------------
> 
> ## From Simon (15 Jul 2009) + dimnames (11 May 2015)
> set.seed(10)
> a <- matrix(round(10 * runif(90)), 10,9, dimnames =
+             list(LETTERS[1:10], paste0("c", 1:9)))
> a[a < 7.5] <- 0
> (A <- Matrix(a))# first column = all zeros
10 x 9 sparse Matrix of class "dgCMatrix"
  c1 c2 c3 c4 c5 c6 c7 c8 c9
A  .  .  9  .  .  .  .  .  9
B  .  .  .  .  .  9  .  .  .
C  .  .  8  .  .  .  .  .  .
D  .  .  .  9  .  .  .  .  .
E  .  .  .  .  .  .  8  .  .
F  .  .  .  .  .  .  8  .  .
G  .  .  8  8  .  .  .  .  8
H  .  .  . 10  .  .  . 10  .
I  .  .  8  .  .  .  .  .  .
J  .  8  .  .  8  .  .  .  .
> qD <- chkQR(a, giveRE=TRUE) ## using base qr
Mean relative difference: 4.60537e-16
Mean relative difference: 3.302202e-16
Mean relative difference: 6.056655e-16
Mean relative difference: 7.322761e-16
Mean relative difference: 6.056655e-16
Mean relative difference: 8.208748e-16
> qS <- chkQR(A, giveRE=TRUE) ## using Matrix "sparse qr" -- "structurally rank deficient!
is sparse *structurally* rank deficient:  Qinv.chk=FALSE, QtQ.chk=FALSE
Mean relative difference: 6.18072e-16
Mean relative difference: 6.18072e-16
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
4: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
5: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
> validObject(qS)# with the validity now (2012-11-18) -- ok, also for "bad" case
[1] TRUE
> ## Here, have illegal access  Up[-1] in ../src/cs.c
> try( ## After patch (2016-10-04 - *NOT* committed), this fails
+     ## definitely "fails" (with good singularity message) after c3194 (cs.c):
+ chk.qr.D.S(qD, qS, y = 10 + 1:nrow(A), force=TRUE)# 6 warnings: "structurally rank deficient"
+ )
Error in chk.qr.D.S(qD, qS, y = 10 + 1:nrow(A), force = TRUE) : 
  is.all.equal3(cc, qr.coef(s., y), drop(qr.coef(s., Y)), tol = tol) is not TRUE
In addition: Warning message:
In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
> try( ## NOTE: *Both* checks  currently fail here:
+     chkQR(A, Qinv.chk=TRUE, QtQ.chk=TRUE)
+ )
Error in assert.EQ(drop(qr.qy(a.qr, qr.qty(a.qr, y))), y, giveRE = giveRE,  : 
  all.equal() |-> names for target but not for current
                  Mean relative difference: 0.2131411
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
> 
> 
> ## Larger Scale random testing
> oo <- options(Matrix.quiet.qr.R = TRUE, Matrix.verbose = TRUE, nwarnings = 1e4)
> set.seed(101)
> 
> quiet <- doExtras
> for(N in 1:(if(doExtras) 1008 else 24)) {
+     A <- rsparsematrix(8,5, nnz = rpois(1, lambda=16))
+     cat(sprintf(if(quiet) "%d " else "%4d -", N)); if(quiet && N %% 50 == 0) cat("\n")
+     checkQR.DS.both(A, Qinv.chk= NA, QtQ.chk=NA, quiet=quiet,
+     ##                          --- => FALSE if struct. rank deficient
+ 		    giveRE = FALSE, tol = 1e-12)
+     ## with doExtras = TRUE, 64bit (F34, R 4.3.0-dev. 2022-05): seen 8.188e-13
+ }
   1 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   2 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   3 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   4 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   5 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   6 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   7 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   8 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   9 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  10 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  11 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  12 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  13 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  14 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  15 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  16 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  17 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  18 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  19 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  20 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  21 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  22 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  23 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  24 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
There were 35 warnings (use warnings() to see them)
> summary(warnings())
Summary of (a total of 35) warning messages:
 6x : In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
24x : In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
 1x : In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
 4x : In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 2 row(s) of zeros
> 
> ## Look at single "hard" cases: --------------------------------------
> 
> ## This is *REALLY* nice and small :
> A0 <- new("dgCMatrix", Dim = 4:3, i = c(0:3, 3L), p = c(0L, 3:5), x = rep(1,5))
> A0
4 x 3 sparse Matrix of class "dgCMatrix"
          
[1,] 1 . .
[2,] 1 . .
[3,] 1 . .
[4,] . 1 1
> checkQR.DS.both(A0, Qinv.chk = FALSE, QtQ.chk=FALSE)
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 2.220446e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.110223e-16
Mean relative difference: 2.775558e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 1.110223e-16
Mean relative difference: 3.191891e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.110223e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 1.110223e-16
[Ok]
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
4: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
5: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> ##                                           ----- *both* still needed :
> try( checkQR.DS.both(A0,  TRUE, FALSE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 2.220446e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.110223e-16
Mean relative difference: 2.775558e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 1.110223e-16
Mean relative difference: 3.191891e-16
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Error in assert.EQ(drop(qr.qy(a.qr, qr.qty(a.qr, y))), y, giveRE = giveRE,  : 
  all.equal() |-> Mean relative difference: 0.09622504
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> try( checkQR.DS.both(A0, FALSE,  TRUE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 2.220446e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.110223e-16
Mean relative difference: 2.775558e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 1.110223e-16
Mean relative difference: 3.191891e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.110223e-16
Error in assert.EQ(MM, m, tol = tol, showOnly = showOnly, giveRE = giveRE) : 
  all.equal() |-> Mean relative difference: 1
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> 
> ## and the same when dropping the first row  { --> 3 x 3 }:
> A1 <- A0[-1 ,]
> checkQR.DS.both(A1, Qinv.chk = FALSE, QtQ.chk=FALSE)
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 4.440892e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 4.440892e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 4.440892e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
4: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
5: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> ##                                           ----- *both* still needed :
> try( checkQR.DS.both(A1,  TRUE, FALSE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 4.440892e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 4.440892e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 4.440892e-16
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Error in assert.EQ(drop(qr.qy(a.qr, qr.qty(a.qr, y))), y, giveRE = giveRE,  : 
  all.equal() |-> Mean relative difference: 0.3333333
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> try( checkQR.DS.both(A1, FALSE,  TRUE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 4.440892e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 4.440892e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 4.440892e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Error in assert.EQ(MM, m, tol = tol, showOnly = showOnly, giveRE = giveRE) : 
  all.equal() |-> Mean relative difference: 1
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> 
> 
> qa <- qr(as(A0,"matrix"))
> qA <- qr(A0) # -> message: ".. Matrix structurally rank deficient"
Warning message:
In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> 
> drop0(crossprod( Qd <- qr.Q(qa) ), 1e-15) # perfect = diag( 3 )
3 x 3 sparse Matrix of class "dsCMatrix"
          
[1,] 1 . .
[2,] . 1 .
[3,] . . 1
> drop0(crossprod( Qs <- qr.Q(qA) ), 1e-15) # R[3,3] == 0 -- OOPS!
3 x 3 sparse Matrix of class "dsCMatrix"
          
[1,] 1 . .
[2,] . 1 .
[3,] . . .
Warning message:
In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> ## OTOH, qr.R() is fine, as checked in the checkQR.DS.both(A0, *) above
> 
> 
> ## zero-row *and* zero-column :
> (A2 <- new("dgCMatrix", i = c(0L, 1L, 4L, 7L, 5L, 2L, 4L)
+            , p = c(0L, 3L, 4L, 4L, 5L, 7L)
+            , Dim = c(8L, 5L)
+            , x = c(0.92, 1.06, -1.74, 0.74, 0.19, -0.63, 0.68)))
8 x 5 sparse Matrix of class "dgCMatrix"
                            
[1,]  0.92 .    . .     .   
[2,]  1.06 .    . .     .   
[3,]  .    .    . .    -0.63
[4,]  .    .    . .     .   
[5,] -1.74 .    . .     0.68
[6,]  .    .    . 0.19  .   
[7,]  .    .    . .     .   
[8,]  .    0.74 . .     .   
> checkQR.DS.both(A2, Qinv.chk = FALSE, QtQ.chk=FALSE)
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 3.757678e-16
Mean relative difference: 5.551115e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.790682e-16
Mean relative difference: 6.661338e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 1.790682e-16
Mean relative difference: 9.992007e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean relative difference: 6.637203e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 6.637203e-16
[Ok]
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
4: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
5: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> ##                                           ----- *both* still needed :
> try( checkQR.DS.both(A2,  TRUE, FALSE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 3.757678e-16
Mean relative difference: 5.551115e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.790682e-16
Mean relative difference: 6.661338e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 1.790682e-16
Mean relative difference: 9.992007e-16
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Error in assert.EQ(drop(qr.qy(a.qr, qr.qty(a.qr, y))), y, giveRE = giveRE,  : 
  all.equal() |-> Mean relative difference: 0.875
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> try( checkQR.DS.both(A2, FALSE,  TRUE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 3.757678e-16
Mean relative difference: 5.551115e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.790682e-16
Mean relative difference: 6.661338e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 1.790682e-16
Mean relative difference: 9.992007e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean relative difference: 6.637203e-16
Error in assert.EQ(MM, m, tol = tol, showOnly = showOnly, giveRE = giveRE) : 
  all.equal() |-> Mean absolute difference: 0.3333333
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> 
> 
> ## Case of *NO* zero-row or zero-column:
> (A3 <- new("dgCMatrix", Dim = 6:5
+            , i = c(0L, 2L, 4L, 0L, 1L, 5L, 1L, 3L, 0L)
+            , p = c(0L, 1L, 3L, 6L, 8L, 9L)
+            , x = c(40, -54, -157, -28, 75, 166, 134, 3, -152)))
6 x 5 sparse Matrix of class "dgCMatrix"
                         
[1,] 40    . -28   . -152
[2,]  .    .  75 134    .
[3,]  .  -54   .   .    .
[4,]  .    .   .   3    .
[5,]  . -157   .   .    .
[6,]  .    . 166   .    .
> checkQR.DS.both(A3, Qinv.chk = FALSE, QtQ.chk=FALSE)
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 4.440892e-16
Mean relative difference: 4.440892e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 6.136406e-16
Mean relative difference: 3.091163e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 6.136406e-16
Mean relative difference: 4.87797e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.143119e-15
complete = TRUE: checking  X = Q R P*
Mean relative difference: 1.143119e-15
[Ok]
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
4: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
5: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> ##                                           ----- *both* still needed :
> try( checkQR.DS.both(A3,  TRUE, FALSE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 4.440892e-16
Mean relative difference: 4.440892e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 6.136406e-16
Mean relative difference: 3.091163e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 6.136406e-16
Mean relative difference: 4.87797e-16
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Error in assert.EQ(drop(qr.qy(a.qr, qr.qty(a.qr, y))), y, giveRE = giveRE,  : 
  all.equal() |-> Mean relative difference: 0.2579847
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> try( checkQR.DS.both(A3, FALSE,  TRUE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 4.440892e-16
Mean relative difference: 4.440892e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 6.136406e-16
Mean relative difference: 3.091163e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 6.136406e-16
Mean relative difference: 4.87797e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean relative difference: 1.143119e-15
Error in assert.EQ(MM, m, tol = tol, showOnly = showOnly, giveRE = giveRE) : 
  all.equal() |-> Mean relative difference: 0.3333333
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> 
> 
> 
> (A4 <- new("dgCMatrix", Dim = c(7L, 5L)
+            , i = c(1:2, 4L, 6L, 1L, 5L, 0:3, 0L, 2:4)
+            , p = c(0L, 4L, 6L, 10L, 10L, 14L)
+            , x = c(9, -8, 1, -9, 1, 10, -1, -2, 6, 14, 10, 2, 12, -9)))
7 x 5 sparse Matrix of class "dgCMatrix"
                  
[1,]  .  . -1 . 10
[2,]  9  1 -2 .  .
[3,] -8  .  6 .  2
[4,]  .  . 14 . 12
[5,]  1  .  . . -9
[6,]  . 10  . .  .
[7,] -9  .  . .  .
> checkQR.DS.both(A4, Qinv.chk = FALSE, QtQ.chk=FALSE)
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 3.089316e-16
Mean relative difference: 3.505967e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 5.298712e-16
Mean relative difference: 6.328271e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 5.298712e-16
Mean relative difference: 7.811212e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean relative difference: 4.234339e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 4.234339e-16
[Ok]
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
4: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
5: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> ##                                           ----- *both* still needed :
> try( checkQR.DS.both(A4,  TRUE, FALSE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 3.089316e-16
Mean relative difference: 3.505967e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 5.298712e-16
Mean relative difference: 6.328271e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 5.298712e-16
Mean relative difference: 7.811212e-16
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Error in assert.EQ(drop(qr.qy(a.qr, qr.qty(a.qr, y))), y, giveRE = giveRE,  : 
  all.equal() |-> Mean relative difference: 0.3289631
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> try( checkQR.DS.both(A4, FALSE,  TRUE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 3.089316e-16
Mean relative difference: 3.505967e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 5.298712e-16
Mean relative difference: 6.328271e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 5.298712e-16
Mean relative difference: 7.811212e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean relative difference: 4.234339e-16
Error in assert.EQ(MM, m, tol = tol, showOnly = showOnly, giveRE = giveRE) : 
  all.equal() |-> Mean relative difference: 0.3333333
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> 
> (A5 <- new("dgCMatrix", Dim = c(4L, 4L)
+            , i = c(2L, 2L, 0:1, 0L, 2:3), p = c(0:2, 4L, 7L)
+            , x = c(48, 242, 88, 18, -167, -179, 18)))
4 x 4 sparse Matrix of class "dgCMatrix"
                   
[1,]  .   . 88 -167
[2,]  .   . 18    .
[3,] 48 242  . -179
[4,]  .   .  .   18
> checkQR.DS.both(A5, Qinv.chk = FALSE, QtQ.chk=FALSE)
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 5.921189e-16
Mean relative difference: 3.552714e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 2.732857e-16
Mean relative difference: 3.700743e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 2.732857e-16
Mean relative difference: 3.700743e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean absolute difference: 7.105427e-15
complete = TRUE: checking  X = Q R P*
Mean absolute difference: 7.105427e-15
[Ok]
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
4: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
5: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> ##                                           ----- *both* still needed :
> try( checkQR.DS.both(A5,  TRUE, FALSE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 5.921189e-16
Mean relative difference: 3.552714e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 2.732857e-16
Mean relative difference: 3.700743e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 2.732857e-16
Mean relative difference: 3.700743e-16
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Error in assert.EQ(drop(qr.qy(a.qr, qr.qty(a.qr, y))), y, giveRE = giveRE,  : 
  all.equal() |-> Mean relative difference: 0.6884769
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> try( checkQR.DS.both(A5, FALSE,  TRUE) )
classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
Mean relative difference: 5.921189e-16
Mean relative difference: 3.552714e-16
complete = FALSE: checking  X = Q R P*
Mean relative difference: 2.732857e-16
Mean relative difference: 3.700743e-16
complete = TRUE: checking  X = Q R P*
Mean relative difference: 2.732857e-16
Mean relative difference: 3.700743e-16
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
Mean absolute difference: 7.105427e-15
Error in assert.EQ(MM, m, tol = tol, showOnly = showOnly, giveRE = giveRE) : 
  all.equal() |-> Mean relative difference: 1
In addition: Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> 
> 
> quiet <- doExtras
> for(N in 1:(if(doExtras) 2^12 else 128)) {
+     A <- round(100*rsparsematrix(5,3, nnz = min(15,rpois(1, lambda=10))))
+     if(any(apply(A, 2, function(x) all(x == 0)))) ## "column of all 0"
+         next
+     cat(sprintf(if(quiet) "%d " else "%4d -", N)); if(quiet && N %% 50 == 0) cat("\n")
+     checkQR.DS.both(A, Qinv.chk=NA, giveRE=FALSE, tol = 1e-12, quiet = quiet)
+     ##                         --- => FALSE if struct. rank deficient
+ }
   1 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   2 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   3 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   4 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   5 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   6 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   7 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   8 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
   9 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  10 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  11 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  12 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  13 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  14 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  15 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  16 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  17 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  18 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  19 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  21 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  22 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  24 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  25 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  26 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  28 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  29 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  30 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  31 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  32 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  34 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  35 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  37 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  38 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  39 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  40 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  41 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  42 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  43 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  44 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  45 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  46 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  47 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  48 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  49 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  50 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  51 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  52 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  53 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  54 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  55 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  56 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  57 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  58 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  59 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  60 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  61 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  62 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  63 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  64 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  65 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  66 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  67 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  68 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  69 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  70 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  71 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  72 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  73 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  74 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  75 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  76 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  77 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  78 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  79 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  80 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  81 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  82 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  83 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  84 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  85 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  86 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  87 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  88 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  89 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  90 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  91 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  92 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  93 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  94 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  95 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  96 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  97 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  98 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
  99 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 100 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 101 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 102 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 103 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 104 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 105 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 106 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 107 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 108 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 109 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 110 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 111 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 112 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 113 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 114 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 115 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 116 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 117 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 118 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 119 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 120 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 121 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 122 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 123 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 124 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 126 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 127 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
 128 -classical: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok] ---  sparse: Qinv.chk=TRUE: checking   Q'Q y = y = QQ' y :
complete = FALSE: checking  X = Q R P*
complete = TRUE: checking  X = Q R P*
[Ok]
Warning messages:
1: In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
2: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
3: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
4: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
5: In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> 
> summary(warnings())
Summary of (a total of 5) warning messages:
1x : In .qr.rank.def.warn(r) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
4x : In .qr.rank.def.warn(qr) :
  matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
> 
> 
> options(oo)
> 
> 
> 
> ### "denseLU"
> 
> ## Testing expansions of factorizations {was ./expand.R, then in simple.R }
> ## new: [m x n]  where m and n  may differ
> x. <- c(2^(0:5),9:1,-3:8, round(sqrt(0:16)))
> set.seed(1)
> for(nnn in 1:100) {
+     y <- sample(x., replace=TRUE)
+     m <- sample(2:6, 1)
+     n <- sample(2:7, 1)
+     x <- matrix(seq_len(m*n), m,n)
+     lux <- lu(x)# occasionally a warning about exact singularity
+     xx <- with(expand(lux), (P %*% L %*% U))
+     print(dim(xx))
+     assert.EQ.mat(xx, x, tol = 16*.Machine$double.eps)
+ }
[1] 4 6
[1] 3 2
[1] 4 5
[1] 6 5
[1] 4 7
[1] 3 3
[1] 5 3
[1] 2 6
[1] 5 6
[1] 5 2
[1] 4 7
[1] 6 5
[1] 4 4
[1] 6 3
[1] 2 7
[1] 5 3
[1] 6 6
[1] 3 4
[1] 6 3
[1] 5 7
[1] 3 3
[1] 6 5
[1] 5 4
[1] 4 2
[1] 6 2
[1] 3 7
[1] 3 5
[1] 3 6
[1] 2 4
[1] 2 5
[1] 4 2
[1] 3 6
[1] 5 4
[1] 6 2
[1] 3 2
[1] 6 6
[1] 5 6
[1] 3 7
[1] 2 5
[1] 6 2
[1] 6 3
[1] 2 4
[1] 5 2
[1] 4 2
[1] 6 2
[1] 3 3
[1] 5 7
[1] 4 4
[1] 4 3
[1] 3 7
[1] 5 7
[1] 3 5
[1] 5 5
[1] 6 2
[1] 3 7
[1] 6 7
[1] 3 3
[1] 3 6
[1] 2 5
[1] 4 2
[1] 3 7
[1] 4 7
[1] 4 7
[1] 2 5
[1] 6 6
[1] 6 4
[1] 5 4
[1] 5 7
[1] 2 3
[1] 5 6
[1] 6 2
[1] 5 7
[1] 6 7
[1] 3 6
[1] 3 5
[1] 6 5
[1] 2 2
[1] 6 4
[1] 4 4
[1] 3 2
[1] 6 3
[1] 5 5
[1] 2 2
[1] 4 7
[1] 2 3
[1] 2 2
[1] 2 5
[1] 5 6
[1] 3 4
[1] 5 5
[1] 2 3
[1] 3 6
[1] 2 4
[1] 3 4
[1] 5 6
[1] 3 5
[1] 5 5
[1] 3 3
[1] 2 4
[1] 4 4
There were 50 or more warnings (use warnings() to see the first 50)
> 
> ### "sparseLU"
> por1 <- readMM(system.file("external/pores_1.mtx", package = "Matrix"))
> lu1 <- lu(por1)
> pm <- as(por1, "CsparseMatrix")
> (pmLU <- lu(pm)) # -> show(<MatrixFactorization>)
LU factorization of Formal class 'sparseLU' [package "Matrix"] with 6 slots
  ..@ L       :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
  .. .. ..@ i       : int [1:141] 0 1 18 22 25 26 1 18 22 25 ...
  .. .. ..@ p       : int [1:31] 0 6 11 13 16 19 22 26 31 35 ...
  .. .. ..@ Dim     : int [1:2] 30 30
  .. .. ..@ Dimnames:List of 2
  .. .. .. ..$ : NULL
  .. .. .. ..$ : NULL
  .. .. ..@ x       : num [1:141] 1 -0.993819 -0.004979 -0.000132 0.000132 ...
  .. .. ..@ uplo    : chr "L"
  .. .. ..@ diag    : chr "N"
  ..@ U       :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
  .. .. ..@ i       : int [1:206] 0 0 1 2 2 3 2 4 2 3 ...
  .. .. ..@ p       : int [1:31] 0 1 3 4 6 8 12 15 20 25 ...
  .. .. ..@ Dim     : int [1:2] 30 30
  .. .. ..@ Dimnames:List of 2
  .. .. .. ..$ : NULL
  .. .. .. ..$ : NULL
  .. .. ..@ x       : num [1:206] -7178502 -24613411 -18311731 -6399179 715 ...
  .. .. ..@ uplo    : chr "U"
  .. .. ..@ diag    : chr "N"
  ..@ p       : int [1:30] 1 11 29 27 19 28 9 18 17 7 ...
  ..@ q       : int [1:30] 0 1 29 27 19 28 9 18 17 8 ...
  ..@ Dim     : int [1:2] 30 30
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
> xp <- expand(pmLU)
> ## permute rows and columns of original matrix
> ppm <- pm[pmLU@p + 1:1, pmLU@q + 1:1]
> Ppm <- pmLU@L %*% pmLU@U
> ## identical only as long as we don't keep the original class info:
> stopifnot(identical3(lu1, pmLU, pm@factors$sparseLU),# TODO === por1@factors$LU
+ 	  identical(ppm, with(xp, P %*% pm %*% t(Q))),
+ 	  sapply(xp, is, class2="Matrix"))
> 
> Ipm <- solve(pm, sparse=FALSE)
> Spm <- solve(pm, sparse=TRUE)  # is not sparse at all, here
> assert.EQ.Mat(Ipm, Spm, giveRE=TRUE, tol = 1e-13)# seen 7.36e-15 only on 32-bit
> stopifnot(abs(as.vector(solve(Diagonal(30, x=10) %*% pm) / Ipm) - 1/10) < 1e-7,
+ 	  abs(as.vector(solve(rep.int(4, 30)	  *  pm) / Ipm) - 1/ 4) < 1e-7)
> 
> ## these two should be the same, and `are' in some ways:
> assert.EQ.mat(ppm, as(Ppm, "matrix"), tol = 1e-14, giveRE=TRUE)
Mean relative difference: 3.426772e-16
> ## *however*
> length(ppm@x)# 180
[1] 180
> length(Ppm@x)# 317 !
[1] 317
> table(Ppm@x == 0)# (194, 123) - has 123 "zero" and 14 ``almost zero" entries

FALSE  TRUE 
  194   123 
> 
> ##-- determinant() and det() --- working via LU ---
> m <- matrix(c(0, NA, 0, NA, NA, 0, 0, 0, 1), 3,3)
> m0 <- rbind(0,cbind(0,m))
> M <- as(m,"Matrix"); M ## "dsCMatrix" ...
3 x 3 sparse Matrix of class "dsCMatrix"
            
[1,]  . NA .
[2,] NA NA .
[3,]  .  . 1
> M0 <- rbind(0, cbind(0, M))
> dM  <- as(M, "denseMatrix")
> dM0 <- as(M0,"denseMatrix")
> try( lum  <- lu(M) )# Err: "near-singular A"
Error in .local(x, ...) : 
  LU factorization of .gCMatrix failed: out of memory or near-singular
> (lum  <- lu(M,  errSing=FALSE))# NA --- *BUT* it is not stored in @factors
[1] NA
> (lum0 <- lu(M0, errSing=FALSE))# NA --- and it is stored in M0@factors[["LU"]]
[1] NA
> ## "FIXME" - TODO: Consider
> replNA <- function(x, value) { x[is.na(x)] <- value ; x }
> (EL.1 <- expand(lu.1 <- lu(M.1 <- replNA(M, -10))))
$P
3 x 3 sparse Matrix of class "pMatrix"
          
[1,] . | .
[2,] | . .
[3,] . . |

$L
3 x 3 sparse Matrix of class "dtCMatrix"
          
[1,] 1 . .
[2,] . 1 .
[3,] . . 1

$U
3 x 3 sparse Matrix of class "dtCMatrix"
              
[1,] -10 -10 .
[2,]   . -10 .
[3,]   .   . 1

$Q
3 x 3 sparse Matrix of class "pMatrix"
          
[1,] | . .
[2,] . | .
[3,] . . |

> ## so it's quite clear how  lu() of the *singular* matrix  M	should work
> ## but it's not supported by the C code in ../src/cs.c which errors out
> stopifnot(all.equal(M.1,  with(EL.1, t(P) %*% L %*% U %*% Q)),
+ 	  is.na(det(M)), is.na(det(dM)),
+ 	  is.na(det(M0)), is.na(det(dM0)) )
Warning message:
In .local(A, ...) :
  CHOLMOD warning 'not positive definite' at file '../Cholesky/t_cholmod_rowfac.c', line 430
> 
> ###________ Cholesky() ________
> 
> ##--------  LDL' ---- small exact examples
> 
> set.seed(1)
> for(n in c(5:12)) {
+     cat("\nn = ",n,"\n-------\n")
+     rr <- mkLDL(n)
+     ##    -------- from 'test-tools.R'
+     stopifnot(all(with(rr, A ==
+ 		       as(L %*% D %*% t(L), "symmetricMatrix"))),
+ 	      all(with(rr, A == tcrossprod(L %*% sqrt(D)))))
+     d <- rr$d.half
+     A <- rr$A
+     .A <- as(A, "TsparseMatrix") # 'factors' slot is retained => do chol() _after_ coercion
+     R <- chol(A)
+     assert.EQ.Mat(R, chol(.A)) # gave infinite recursion
+     print(d. <- diag(R))
+     D. <- Diagonal(x= d.^2)
+     L. <- t(R) %*% Diagonal(x = 1/d.)
+     stopifnot(all.equal(as.matrix(D.), as.matrix(rr$ D)),
+               all.equal(as.matrix(L.), as.matrix(rr$ L)))
+     ##
+     CAp <- Cholesky(A)# perm=TRUE --> Permutation:
+     validObject(CAp)
+     p <- CAp@perm + 1L
+     P <- as(p, "pMatrix")
+     ## the inverse permutation:
+     invP <- solve(P)@perm
+     lDet <- sum(2* log(d))# the "true" value
+     ldetp  <-         .diag.dsC(Chx = CAp, res.kind = "sumLog")
+     ldetp. <- sum(log(.diag.dsC(Chx = CAp, res.kind = "diag") ))
+     ##
+     CA	<- Cholesky(A,perm=FALSE)
+     validObject(CA)
+     ldet <- .diag.dsC(Chx = CA, res.kind = "sumLog")
+     ## not printing CAp : ends up non-integer for n >= 11
+     mCAp <- as(CAp, "CsparseMatrix")
+     print(mCA  <- drop0(as(CA, "CsparseMatrix")))
+     stopifnot(identical(A[p,p], as(P %*% A %*% t(P),
+ 				   "symmetricMatrix")),
+ 	      relErr(d.^2, .diag.dsC(Chx= CA, res.kind="diag")) < 1e-14,
+ 	      relErr(A[p,p], tcrossprod(mCAp)) < 1e-14)
+     if(FALSE)
+         rbind(lDet,ldet, ldetp, ldetp.)
+     ##  ==> Empirically, I see  lDet = ldet  !=  ldetp == ldetp.
+     ## if(rr$rcond.A < ...) warning("condition number of A ..." ## <- TODO
+     cat(1,""); assert.EQ.(lDet,  ldet,   tol = 1e-14)
+     cat(2,""); assert.EQ.(ldetp, ldetp., tol = 1e-14)
+     cat(3,""); assert.EQ.(lDet,  ldetp, tol = n^2* 1e-7)# extreme: have seen 0.0011045 !!
+ }## for()

n =  5 
-------
[1] 20 50 40 30 10
5 x 5 sparse Matrix of class "dtCMatrix"
                     
[1,]  20  .   .  .  .
[2,] 100 50   .  .  .
[3,]   .  .  40  .  .
[4,] 140  . 160 30  .
[5,] 120  .   .  . 10
1 2 Mean relative difference: 2.179523e-16
3 Mean relative difference: 2.179523e-16

n =  6 
-------
[1] 30 40 20 60 50 10
6 x 6 sparse Matrix of class "dtCMatrix"
                         
[1,]  30   .  .   .  .  .
[2,]   .  40  .   .  .  .
[3,]   . 240 20   .  .  .
[4,]   .   .  .  60  .  .
[5,] 150 360  .   . 50  .
[6,]   . 280  . 600  . 10
1 Mean relative difference: 1.741974e-16
2 3 
n =  7 
-------
[1] 50 30 10 40 60 20 70
7 x 7 sparse Matrix of class "dtCMatrix"
                            
[1,] 50  .   .   .   .  .  .
[2,]  . 30   .   .   .  .  .
[3,]  .  .  10   .   .  .  .
[4,]  .  . 130  40   .  .  .
[5,]  .  .   .   .  60  .  .
[6,]  .  .  20 480 480 20  .
[7,]  .  .   .   .   .  . 70
1 Mean relative difference: 1.441658e-16
2 3 
n =  8 
-------
[1] 10 80 20 30 50 60 40 70
8 x 8 sparse Matrix of class "dtCMatrix"
                                  
[1,] 10   .   .   .   .   .   .  .
[2,]  .  80   .   .   .   .   .  .
[3,]  .   .  20   .   .   .   .  .
[4,]  .   . 320  30   .   .   .  .
[5,] 70   .  40 150  50   .   .  .
[6,]  .   .   .  30 850  60   .  .
[7,]  .   .   .   .   . 600  40  .
[8,]  . 320   .   .   .   . 440 70
1 2 Mean relative difference: 1.224007e-16
3 Mean relative difference: 1.980492e-11

n =  9 
-------
[1] 30 70 20 50 10 60 40 90 80
9 x 9 sparse Matrix of class "dtCMatrix"
                                        
 [1,]  30   .   .   .   .  .   .    .  .
 [2,]   .  70   .   .   .  .   .    .  .
 [3,]   .   .  20   .   .  .   .    .  .
 [4,]   .   . 420  50   .  .   .    .  .
 [5,]   .   .   . 300  10  .   .    .  .
 [6,] 780 980   .   .   . 60   .    .  .
 [7,]   .   . 540   .   .  .  40    .  .
 [8,]   .   . 200   . 130 60 120   90  .
 [9,] 450   .   . 250  80  . 160 1980 80
1 2 3 Mean relative difference: 1.857049e-12

n =  10 
-------
 [1]  70  90  40  80  20 100  60  50  10  30
10 x 10 sparse Matrix of class "dtCMatrix"
                                               
 [1,]   70    .   .    .   .   .    .  .   .  .
 [2,]    .   90   .    .   .   .    .  .   .  .
 [3,]    .    .  40    .   .   .    .  .   .  .
 [4,]  280 2430   .   80   .   .    .  .   .  .
 [5,]    . 1620   . 2000  20   .    .  .   .  .
 [6,] 2030    . 960    . 220 100    .  .   .  .
 [7,]    .    .   .    .   . 200   60  .   .  .
 [8,]  210    .   . 2080   .   .    . 50   .  .
 [9,] 1400    .   . 1360 180   .    .  .  10  .
[10,]    . 2700   .    .   . 100 1260  . 100 30
1 Mean relative difference: 1.863461e-16
2 3 Mean relative difference: 3.64493e-13

n =  11 
-------
 [1]  50  80 100  40  10  30  90  20 110  60  70
11 x 11 sparse Matrix of class "dtCMatrix"
                                                    
 [1,]   50    .    .   .   .    .    .   .   .  .  .
 [2,]    .   80    .   .   .    .    .   .   .  .  .
 [3,]    .    .  100   .   .    .    .   .   .  .  .
 [4,]    .    . 2800  40   .    .    .   .   .  .  .
 [5,]    .    .    .   .  10    .    .   .   .  .  .
 [6,]    . 2560    .   . 340   30    .   .   .  .  .
 [7,]    .    . 2200   .   . 1110   90   .   .  .  .
 [8,]    .    .    . 760   . 1200    .  20   .  .  .
 [9,] 1950    . 2600   .   .    . 1620   . 110  .  .
[10,]    . 1040    .   .   .    .    . 100   . 60  .
[11,] 1000  720    .   .   .    .  360   .   .  . 70
1 Mean relative difference: 1.658955e-16
2 3 Mean relative difference: 2.296093e-10

n =  12 
-------
 [1] 120 100  40  60  80  70  30  20  10 110  50  90
12 x 12 sparse Matrix of class "dtCMatrix"
                                                           
 [1,]  120    .    .    .    .    .    .   .  .   .    .  .
 [2,]    .  100    .    .    .    .    .   .  .   .    .  .
 [3,]    .    .   40    .    .    .    .   .  .   .    .  .
 [4,]    .    . 1280   60    .    .    .   .  .   .    .  .
 [5,]    . 2300    .    .   80    .    .   .  .   .    .  .
 [6,]    . 1900 1520 2040    .   70    .   .  .   .    .  .
 [7,]  360  700    .    .    .    .   30   .  .   .    .  .
 [8,]    .    .    .    . 1920    .    .  20  .   .    .  .
 [9,]    .    .    .    . 1120    .    . 220 10   .    .  .
[10,]    .    .    .    . 2880    .    .   .  . 110    .  .
[11,] 5640    . 1480    .  320 2030    .   .  .   .   50  .
[12,]    .    .    .    .    .    . 1170   .  . 550 2000 90
1 Mean relative difference: 1.492165e-16
2 Mean relative difference: 2.98433e-16
3 Mean relative difference: 2.671423e-12
> 
> mkCholhash <- function(r.all) {
+     ## r.all %*% (2^(2:0)), but only those that do not have NA / "?" :
+     stopifnot(is.character(rn <- rownames(r.all)),
+               is.matrix(r.all), is.logical(r.all))
+     c.rn <- vapply(rn, function(ch) strsplit(ch, " ")[[1]], character(3))
+     ## Now
+     h1 <- function(i) {
+         ok <- rep.int(TRUE, 3L)
+         if(c.rn[3L, i] == "?")
+             ok[2:3] <- FALSE # no supernodal LDL' factorization !!
+         r.all[i, ok] %*% 2^((2:0)[ok])
+     }
+     vapply(seq_len(nrow(r.all)), h1, numeric(1))
+ }
> 
> set.seed(17)
> (rr <- mkLDL(4))
$A
4 x 4 sparse Matrix of class "dsCMatrix"
                     
[1,] 100   .  300   .
[2,]   . 900  900   .
[3,] 300 900 3400   .
[4,]   .   .    . 400

$L
4 x 4 sparse Matrix of class "dtCMatrix"
            
[1,] 1 . . .
[2,] . 1 . .
[3,] 3 1 1 .
[4,] . . . 1

$d.half
[1] 10 30 40 20

$D
4 x 4 diagonal matrix of class "ddiMatrix"
     [,1] [,2] [,3] [,4]
[1,]  100    .    .    .
[2,]    .  900    .    .
[3,]    .    . 1600    .
[4,]    .    .    .  400

$rcond.A
[1] 0.0112202

$cond.A
NULL

> (CA <- Cholesky(rr$A))
Cholesky factorization of Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
  ..@ x       : num [1:6] 900 1 100 3 1600 400
  ..@ p       : int [1:5] 0 2 4 5 6
  ..@ i       : int [1:6] 0 2 1 2 2 3
  ..@ nz      : int [1:4] 2 2 1 1
  ..@ nxt     : int [1:6] 1 2 3 4 -1 0
  ..@ prv     : int [1:6] 5 0 1 2 3 -1
  ..@ type    : int [1:6] 2 0 0 1 0 0
  ..@ colcount: int [1:4] 2 2 1 1
  ..@ perm    : int [1:4] 1 0 2 3
  ..@ Dim     : int [1:2] 4 4
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
> validObject(CA)
[1] TRUE
> stopifnot(all.equal(determinant(rr$A) -> detA,
+                     determinant(as(rr$A, "matrix"))),
+           is.all.equal3(c(detA$modulus), log(det(rr$D)), sum(log(rr$D@x))))
> A12 <- mkLDL(12, 1/10)
> (r12 <- allCholesky(A12$A))[-1]
Warning in allCholesky(A12$A) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
$dup.r.all
 [1] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

$r.all
       
  p L s  perm   LDL super
  . . . FALSE FALSE FALSE
  | . .  TRUE FALSE FALSE
  . | . FALSE  TRUE FALSE
  | | .  TRUE  TRUE FALSE
  . . | FALSE FALSE  TRUE
  | . |  TRUE FALSE  TRUE
  . | | FALSE FALSE  TRUE
  | | |  TRUE FALSE  TRUE
  . . ? FALSE FALSE FALSE
  | . ?  TRUE FALSE FALSE
  . | ? FALSE  TRUE FALSE
  | | ?  TRUE  TRUE FALSE

$r.uniq
       
  p L s  perm   LDL super
  . . . FALSE FALSE FALSE
  | . .  TRUE FALSE FALSE
  . | . FALSE  TRUE FALSE
  | | .  TRUE  TRUE FALSE
  . . | FALSE FALSE  TRUE
  | . |  TRUE FALSE  TRUE

> aCh.hash <- mkCholhash(r12$r.all)
> if(requireNamespace("sfsmisc"))
+ split(rownames(r12$r.all), sfsmisc::Duplicated(aCh.hash))
Loading required namespace: sfsmisc
$`1`
[1] ". . |" ". | |"

$`2`
[1] "| . |" "| | |"

$`3`
[1] ". . ." ". . ?" ". | ?"

$`4`
[1] "| . ." "| . ?" "| | ?"

> 
> ## TODO: find cases for both choices when we leave it to CHOLMOD to choose
> for(n in 1:50) { ## used to seg.fault at n = 10 !
+     mkA <- mkLDL(1+rpois(1, 30), 1/10, rcond = FALSE, condest = FALSE)
+     cat(sprintf("n = %3d, LDL-dim = %d x %d ", n, nrow(mkA$A), ncol(mkA$A)))
+     r <- allCholesky(mkA$A, silentTry=TRUE)
+     ## Compare .. apart from the NAs that happen from (perm=FALSE, super=TRUE)
+     iNA <- apply(is.na(r$r.all), 1, any)
+     cat(sprintf(" -> %3s NAs\n", if(any(iNA)) format(sum(iNA)) else "no"))
+     stopifnot(aCh.hash[!iNA] == mkCholhash(r$r.all[!iNA,]))
+ ##     cat("--------\n")
+ }
n =   1, LDL-dim = 27 x 27 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =   2, LDL-dim = 27 x 27 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =   3, LDL-dim = 23 x 23 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =   4, LDL-dim = 25 x 25 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =   5, LDL-dim = 32 x 32 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =   6, LDL-dim = 43 x 43 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   6 NAs
n =   7, LDL-dim = 34 x 34 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   3 NAs
n =   8, LDL-dim = 28 x 28 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =   9, LDL-dim = 33 x 33 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  10, LDL-dim = 32 x 32 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  11, LDL-dim = 32 x 32 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  12, LDL-dim = 33 x 33 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   4 NAs
n =  13, LDL-dim = 35 x 35 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  14, LDL-dim = 35 x 35 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   4 NAs
n =  15, LDL-dim = 30 x 30 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   1 NAs
n =  16, LDL-dim = 36 x 36 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   1 NAs
n =  17, LDL-dim = 29 x 29 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  18, LDL-dim = 35 x 35 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  19, LDL-dim = 32 x 32 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   2 NAs
n =  20, LDL-dim = 25 x 25 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  21, LDL-dim = 46 x 46 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   1 NAs
n =  22, LDL-dim = 28 x 28 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  23, LDL-dim = 28 x 28 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  24, LDL-dim = 41 x 41 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   4 NAs
n =  25, LDL-dim = 31 x 31 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   2 NAs
n =  26, LDL-dim = 34 x 34 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  27, LDL-dim = 31 x 31 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  28, LDL-dim = 33 x 33 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  29, LDL-dim = 44 x 44 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   3 NAs
n =  30, LDL-dim = 33 x 33 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  31, LDL-dim = 31 x 31 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   2 NAs
n =  32, LDL-dim = 33 x 33 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  33, LDL-dim = 34 x 34 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  34, LDL-dim = 32 x 32 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  35, LDL-dim = 28 x 28 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  36, LDL-dim = 27 x 27 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  37, LDL-dim = 26 x 26 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  38, LDL-dim = 32 x 32 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  39, LDL-dim = 28 x 28 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  40, LDL-dim = 27 x 27 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  41, LDL-dim = 23 x 23 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  42, LDL-dim = 39 x 39 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   6 NAs
n =  43, LDL-dim = 23 x 23 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  44, LDL-dim = 48 x 48 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   6 NAs
n =  45, LDL-dim = 31 x 31 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  46, LDL-dim = 37 x 37 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  47, LDL-dim = 34 x 34 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   1 NAs
n =  48, LDL-dim = 24 x 24 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
n =  49, LDL-dim = 39 x 39 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->   2 NAs
n =  50, LDL-dim = 24 x 24 Warning in allCholesky(mkA$A, silentTry = TRUE) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
 ->  no NAs
There were 48 warnings (use warnings() to see them)
> 
> 
> ## This is a relatively small "critical example" :
> A. <-
+     new("dsCMatrix", Dim = c(25L, 25L), uplo = "U"
+ 	, i = as.integer(
+           c(0, 1, 2, 3, 4, 2, 5, 6, 0, 8, 8, 9, 3, 4, 10, 11, 6, 12, 13, 4,
+             10, 14, 15, 1, 2, 5, 16, 17, 0, 7, 8, 18, 9, 19, 10, 11, 16, 20,
+             0, 6, 7, 16, 17, 18, 20, 21, 6, 9, 12, 14, 19, 21, 22, 9, 11, 19,
+             20, 22, 23, 1, 16, 24))
+ 	##
+ 	, p = c(0:6, 8:10, 12L, 15:16, 18:19, 22:23, 27:28, 32L, 34L, 38L, 46L, 53L, 59L, 62L)
+ 	##
+ 	, x = c(1, 1, 1, 1, 2, 100, 2, 40, 1, 2, 100, 6700, 100, 100, 13200,
+ 	  1, 50, 4100, 1, 5, 400, 20, 1, 40, 100, 5600, 9100, 5000, 5,
+ 	  100, 100, 5900, 100, 6200, 30, 20, 9, 2800, 1, 100, 8, 10, 8000,
+ 	  100, 600, 23900, 30, 100, 2800, 50, 5000, 3100, 15100, 100, 10,
+ 	  5600, 800, 4500, 5500, 7, 600, 18200))
> validObject(A.)
[1] TRUE
> ## A1: the same pattern as  A.   just simply filled with '1's :
> A1 <- A.; A1@x[] <- 1; A1@factors <- list()
> A1.8 <- A1; diag(A1.8) <- 8
> ##
> nT. <- as(AT <- as(A., "TsparseMatrix"),"nMatrix")
> stopifnot(all(nT.@i <= nT.@j),
+ 	  identical(qr(A1.8), qr(as(A1.8, "generalMatrix"))))
> 
> CA <- Cholesky(A. + Diagonal(x = rowSums(abs(A.)) + 1))
> validObject(CA)
[1] TRUE
> stopifnotValid(CAinv <- solve(CA), "dsCMatrix")
> MA <- as(CA, "CsparseMatrix") # with a confusing warning -- FIXME!
> stopifnotValid(MAinv <- solve(MA), "dtCMatrix")
> ## comparing MAinv with some solve(CA, system="...") .. *not* trivial? - TODO
> ##
> CAinv2 <- solve(CA, Diagonal(nrow(A.)))
> CAinv2 <- as(CAinv2, "symmetricMatrix")
> stopifnot(identical(CAinv, CAinv2))
> 
> ## FINALLY fix "TODO": (not implemented *symbolic* factorization of nMatrix)
> try(    tc <- Cholesky(nT.)  )
Error in .local(A, ...) : leading principal minor of order 8 is zero
In addition: Warning message:
In .local(A, ...) :
  CHOLMOD warning 'not positive definite' at file '../Cholesky/t_cholmod_rowfac.c', line 430
> 
> for(p in c(FALSE,TRUE))
+     for(L in c(FALSE,TRUE))
+         for(s in c(FALSE,TRUE, NA)) {
+             cat(sprintf("p,L,S = (%2d,%2d,%2d): ", p,L,s))
+             r <- tryCatch(Cholesky(A., perm=p, LDL=L, super=s),
+                           error = function(e)e)
+             cat(if(inherits(r, "error")) " *** E ***" else
+                 sprintf("%3d", r@type),"\n", sep="")
+         }
p,L,S = ( 0, 0, 0):  *** E ***
p,L,S = ( 0, 0, 1):  *** E ***
p,L,S = ( 0, 0,NA):  *** E ***
p,L,S = ( 0, 1, 0):   0  0  0  1  0  0
p,L,S = ( 0, 1, 1):  *** E ***
p,L,S = ( 0, 1,NA):   0  0  0  1  0  0
p,L,S = ( 1, 0, 0):  *** E ***
p,L,S = ( 1, 0, 1):  *** E ***
p,L,S = ( 1, 0,NA):  *** E ***
p,L,S = ( 1, 1, 0):   2  0  0  1  0  0
p,L,S = ( 1, 1, 1):  *** E ***
p,L,S = ( 1, 1,NA):   2  0  0  1  0  0
Warning messages:
1: In .local(A, ...) :
  CHOLMOD warning 'not positive definite' at file '../Cholesky/t_cholmod_rowfac.c', line 430
2: In .local(A, ...) :
  CHOLMOD warning 'matrix not positive definite' at file '../Supernodal/t_cholmod_super_numeric.c', line 911
3: In .local(A, ...) :
  CHOLMOD warning 'not positive definite' at file '../Cholesky/t_cholmod_rowfac.c', line 430
4: In .local(A, ...) :
  CHOLMOD warning 'matrix not positive definite' at file '../Supernodal/t_cholmod_super_numeric.c', line 911
5: In .local(A, ...) :
  CHOLMOD warning 'not positive definite' at file '../Cholesky/t_cholmod_rowfac.c', line 430
6: In .local(A, ...) :
  CHOLMOD warning 'matrix not positive definite' at file '../Supernodal/t_cholmod_super_numeric.c', line 911
7: In .local(A, ...) :
  CHOLMOD warning 'not positive definite' at file '../Cholesky/t_cholmod_rowfac.c', line 430
8: In .local(A, ...) :
  CHOLMOD warning 'matrix not positive definite' at file '../Supernodal/t_cholmod_super_numeric.c', line 911
> str(A., max.level=3) ## look at the 'factors'
Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
  ..@ i       : int [1:62] 0 1 2 3 4 2 5 6 0 8 ...
  ..@ p       : int [1:26] 0 1 2 3 4 5 6 8 9 10 ...
  ..@ Dim     : int [1:2] 25 25
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:62] 1 1 1 1 2 100 2 40 1 2 ...
  ..@ uplo    : chr "U"
  ..@ factors :List of 2
  .. ..$ spDCholesky:Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
  .. ..$ sPDCholesky:Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
> 
> facs <- A.@factors
> names(facs) <- sub("Cholesky$", "", names(facs))
> facs <- facs[order(names(facs))]
> 
> sapply(facs, class)
        sPD         spD 
"dCHMsimpl" "dCHMsimpl" 
> str(lapply(facs, slot, "type"))
List of 2
 $ sPD: int [1:6] 2 0 0 1 0 0
 $ spD: int [1:6] 0 0 0 1 0 0
> ## super = TRUE  currently always entails  LDL=FALSE :
> ## hence isLDL is TRUE for ("D" and not "S"):
> sapply(facs, isLDL)
 sPD  spD 
TRUE TRUE 
> 
> chkCholesky <- function(chmf, A) {
+     stopifnot(is(chmf, "CHMfactor"),
+               validObject(chmf),
+               is(A, "Matrix"), isSymmetric(A))
+     if(!is(A, "dsCMatrix"))
+         A <- as(as(as(A, "CsparseMatrix"), "symmetricMatrix", "dMatrix"))
+     L <- drop0(zapsmall(L. <- as(chmf, "CsparseMatrix")))
+     cat("no. nonzeros in L {before / after drop0(zapsmall(.))}: ",
+         c(nnzero(L.), nnzero(L)), "\n") ## 112, 95
+     ecc <- expand(chmf)
+     A... <- with(ecc, crossprod(crossprod(L,P)))
+     stopifnot(all.equal(L., ecc$L, tolerance = 1e-14),
+               all.equal(A,  A...,  tolerance = 1e-14))
+     invisible(ecc)
+ }
> 
> c1.8 <- try(Cholesky(A1.8, super = TRUE))# works "always", interestingly ...
> chkCholesky(c1.8, A1.8)
no. nonzeros in L {before / after drop0(zapsmall(.))}:  71 71 
> 
> 
> 
> ## --- now a "large" (712 x 712) real data example ---------------------------
> 
> data(KNex, package = "Matrix")
> mtm <- with(KNex, crossprod(mm))
> ld.3 <- determinant(Cholesky(mtm, perm = TRUE), sqrt = FALSE)
> stopifnot(identical(names(mtm@factors),
+                     "sPDCholesky"))
> ld.4 <- determinant(Cholesky(mtm, perm = FALSE), sqrt = FALSE)
> stopifnot(identical(names(mtm@factors),
+                     c("sPDCholesky", "spDCholesky")))
> c2 <- Cholesky(mtm, super = TRUE)
> validObject(c2)
[1] TRUE
> stopifnot(identical(names(mtm@factors),
+                     c("sPDCholesky", "spDCholesky", "SPdCholesky")))
> 
> r <- allCholesky(mtm)
Warning in allCholesky(mtm) :
  duplicated( <pLs-matrix> ) differs from duplicated( <CHM-list> )
> r[-1]
$dup.r.all
 [1] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

$r.all
       
  p L s  perm   LDL super
  . . . FALSE FALSE FALSE
  | . .  TRUE FALSE FALSE
  . | . FALSE  TRUE FALSE
  | | .  TRUE  TRUE FALSE
  . . | FALSE FALSE  TRUE
  | . |  TRUE FALSE  TRUE
  . | | FALSE FALSE  TRUE
  | | |  TRUE FALSE  TRUE
  . . ? FALSE FALSE FALSE
  | . ?  TRUE FALSE FALSE
  . | ? FALSE  TRUE FALSE
  | | ?  TRUE  TRUE FALSE

$r.uniq
       
  p L s  perm   LDL super
  . . . FALSE FALSE FALSE
  | . .  TRUE FALSE FALSE
  . | . FALSE  TRUE FALSE
  | | .  TRUE  TRUE FALSE
  . . | FALSE FALSE  TRUE
  | . |  TRUE FALSE  TRUE

> 
> ## is now taken from cache
> c1 <- Cholesky(mtm)
> 
> bv <- 1:nrow(mtm) # even integer
> b <- matrix(bv)
> ## solve(c2, b) by default solves  Ax = b, where A = c2'c2 !
> x <- solve(c2,b)
> stopifnot(identical3(drop(x), solve(c2, bv), drop(solve(c2, b, system = "A"))),
+           all.equal(x, solve(mtm, b)))
> for(sys in c("A", "LDLt", "LD", "DLt", "L", "Lt", "D", "P", "Pt")) {
+     x <- solve(c2, b,  system = sys)
+     cat(sys,":\n"); print(head(x))
+     stopifnot(dim(x) == c(712, 1),
+               identical(drop(x), solve(c2, bv, system = sys)))
+ }
A :
6 x 1 Matrix of class "dgeMatrix"
          [,1]
[1,] -324576.9
[2,] -602422.0
[3,] -475225.1
[4,] -656840.4
[5,] -823533.6
[6,] -853756.8
LDLt :
6 x 1 Matrix of class "dgeMatrix"
          [,1]
[1,] 37696.471
[2,]  4713.432
[3,]  6643.932
[4,]  1866.030
[5,]  4340.785
[6,]  2086.379
LD :
6 x 1 Matrix of class "dgeMatrix"
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
[5,]    5
[6,]    6
DLt :
6 x 1 Matrix of class "dgeMatrix"
          [,1]
[1,]  778.7719
[2,] -422.7658
[3,] -351.6316
[4,] -522.6810
[5,] -544.7115
[6,]  227.2868
L :
6 x 1 Matrix of class "dgeMatrix"
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
[5,]    5
[6,]    6
Lt :
6 x 1 Matrix of class "dgeMatrix"
          [,1]
[1,]  778.7719
[2,] -422.7658
[3,] -351.6316
[4,] -522.6810
[5,] -544.7115
[6,]  227.2868
D :
6 x 1 Matrix of class "dgeMatrix"
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
[5,]    5
[6,]    6
P :
6 x 1 Matrix of class "dgeMatrix"
     [,1]
[1,]  257
[2,]  244
[3,]  243
[4,]  242
[5,]  214
[6,]  694
Pt :
6 x 1 Matrix of class "dgeMatrix"
     [,1]
[1,]  617
[2,]  469
[3,]  473
[4,]  619
[5,]  453
[6,]  458
> 
> ## log(|LL'|) - check if super = TRUE and simplicial give same determinant
> (ld.1 <- determinant(mtm))
$modulus
[1] -343.1384
attr(,"logarithm")
[1] TRUE

$sign
[1] 1

attr(,"class")
[1] "det"
> if(FALSE) {
+ ## MJ: CHMfactor_ldetL2 is unused outside of these tests, so we no longer
+ ##     have it in the namespace { old definition is in ../src/CHMfactor.c }
+ ld1 <- .Call("CHMfactor_ldetL2", c1)
+ ld2 <- .Call("CHMfactor_ldetL2", c2)
+ stopifnot(all.equal(ld1, ld2),
+ 	  all.equal(ld1, as.vector(ld.1$modulus), tolerance = 1e-14),
+           all.equal(ld1, as.vector(ld.3$modulus), tolerance = 1e-14),
+           all.equal(ld1, as.vector(ld.4$modulus), tolerance = 1e-14))
+ } else {
+ stopifnot(all.equal(as.vector(ld.1$modulus), as.vector(ld.3$modulus),
+                     tolerance = 1e-14),
+           all.equal(as.vector(ld.1$modulus), as.vector(ld.4$modulus),
+                     tolerance = 1e-14))
+ }
> 
> ## MJ: ldet[123].dsC() are unused outside of these tests, so we no longer
> ##     have them in the namespace { old definitions are in ../R/determinant.R }
> if(FALSE) {
+ ## Some timing measurements
+ mtm <- with(KNex, crossprod(mm))
+ I <- .symDiagonal(n=nrow(mtm))
+ set.seed(101); r <- runif(100)
+ 
+ system.time(D1 <- sapply(r, function(rho) Matrix:::ldet1.dsC(mtm + (1/rho) * I)))
+ ## 0.842 on fast cmath-5
+ system.time(D2 <- sapply(r, function(rho) Matrix:::ldet2.dsC(mtm + (1/rho) * I)))
+ ## 0.819
+ system.time(D3 <- sapply(r, function(rho) Matrix:::ldet3.dsC(mtm + (1/rho) * I)))
+ ## 0.810
+ stopifnot(is.all.equal3(D1,D2,D3, tol = 1e-13))
+ }
> 
> ## Updating LL'  should remain LL' and not become  LDL' :
> cholCheck <- function(Ut, tol = 1e-12, super = FALSE, LDL = !super) {
+     L <- Cholesky(UtU <- tcrossprod(Ut), super=super, LDL=LDL, Imult = 1)
+     L1 <- update(L, UtU, mult = 1)
+     L2 <- update(L, Ut,  mult = 1)
+     stopifnot(is.all.equal3(L, L1, L2, tol = tol),
+               all.equal(update(L, UtU, mult = pi),
+                         update(L, Ut,  mult = pi), tolerance = tol)
+               )
+ }
> 
> ## Inspired by
> ## data(Dyestuff, package = "lme4")
> ## Zt <- as(Dyestuff$Batch, "sparseMatrix")
> Zt <- new("dgCMatrix", Dim = c(6L, 30L), x = 2*1:30,
+           i = rep(0:5, each=5),
+           p = 0:30, Dimnames = list(LETTERS[1:6], NULL))
> cholCheck(0.78 * Zt, tol=1e-14)
> 
> oo <- options(Matrix.quiet.qr.R = TRUE, warn = 2)# no warnings allowed
> qrZ <- qr(t(Zt))
> Rz <- qr.R(qrZ)
> stopifnot(exprs = {
+     inherits(qrZ, "sparseQR")
+     inherits(Rz, "sparseMatrix")
+     isTriangular(Rz)
+     isDiagonal(Rz) # even though formally a "dtCMatrix"
+     qr2rankMatrix(qrZ, do.warn=FALSE) == 6
+ })
> options(oo)
> 
> ## problematic rank deficient rankMatrix() case -- only seen in large cases ??
> ## MJ: NA in diag(<sparseQR>@R) not seen with Apple Clang 14.0.3
> Z. <- readRDS(system.file("external", "Z_NA_rnk.rds", package="Matrix"))
> (rnkZ. <- rankMatrix(Z., method = "qr")) # gave errors; now warns typically, but not on aarm64 (M1)
[1] NA
attr(,"method")
[1] "qr.R"
attr(,"useGrad")
[1] FALSE
attr(,"tol")
[1] 3.126388e-13
Warning message:
In qr2rankMatrix(q.r, tol = tol, isBqr = x.dense, do.warn = warn.qr) :
  qr2rankMatrix(.): QR with only 684 out of 822 finite diag(R) entries
> qrZ. <- qr(Z.)
> options(warn=1)
> rnk2 <- qr2rankMatrix(qrZ.) # warning ".. only 684 out of 822 finite diag(R) entries"
Warning in qr2rankMatrix(qrZ.) :
  qr2rankMatrix(.): QR with only 684 out of 822 finite diag(R) entries
> oo <- options(warn=2)# no warnings allowed from here
> di.NA <- anyNA(diag(qrZ.@R))
> stopifnot(is(qrZ, "sparseQR"),
+           identical(is.na(rnkZ.), di.NA),
+           identical(is.na(rnk2), di.NA))
> 
> ## The above bug fix was partly wrongly extended to  dense matrices for "qr.R":
> x <- cbind(1, rep(0:9, 18))
> qr.R(qr(x))              # one negative diagonal
          [,1]      [,2]
[1,] -13.41641 -60.37384
[2,]   0.00000  38.53570
> qr.R(qr(x, LAPACK=TRUE)) # two negative diagonals
          [,1]       [,2]
[1,] -71.62402 -11.309056
[2,]   0.00000  -7.218398
> chkRnk <- function(x, rnk) {
+     stopifnot(exprs = {
+         rankMatrix(x) == rnk
+         rankMatrix(x, method="maybeGrad") == rnk ## but "useGrad" is not !
+         rankMatrix(x, method="qrLINPACK") == rnk
+         rankMatrix(x, method="qr.R"     ) == rnk
+     })# the last gave '0' and a warning in Matrix 1.3-0
+ }
> chkRnk(   x,    2)
> chkRnk(diag(1), 1) # had "empty stopifnot" (-> Error in MM's experimental setup) +  warning 'min(<empty>)'
> (m3 <- cbind(2, rbind(diag(pi, 2), 8)))
     [,1]     [,2]     [,3]
[1,]    2 3.141593 0.000000
[2,]    2 0.000000 3.141593
[3,]    2 8.000000 8.000000
> chkRnk(m3, 3)
> chkRnk(matrix(0, 4,3), 0)
> chkRnk(matrix(1, 5,5), 1) # had failed for "maybeGrad"
> chkRnk(matrix(1, 5,2), 1)
> 
> 
> showSys.time(
+ for(i in 1:120) {
+     set.seed(i)
+     M <- rspMat(n=rpois(1,50), m=rpois(1,20), density = 1/(4*rpois(1, 4)))
+     cat(sprintf("%3d: dim(M) = %2dx%2d, rank=%2d, k=%9.4g; ",
+ 		i, nrow(M), ncol(M), rankMatrix(M), kappa(M)))
+     for(super in c(FALSE,TRUE)) {
+         cat("super=",super,"M: ")
+         ## 2018-01-04, Avi Adler: needed 1.2e-12 in Windows 64 (for i=55, l.1):
+         cholCheck( M  , tol=2e-12, super=super); cat(" M': ")
+         cholCheck(t(M), tol=2e-12, super=super)
+     }
+     cat(" [Ok]\n")
+ })
  1: dim(M) = 45x25, rank=21, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
  2: dim(M) = 43x15, rank=15, k=    5.083; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
  3: dim(M) = 43x18, rank=17, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
  4: dim(M) = 51x17, rank=17, k=    26.07; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
  5: dim(M) = 44x17, rank=16, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
  6: dim(M) = 51x17, rank=12, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
  7: dim(M) = 66x14, rank=14, k=    5.556; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
  8: dim(M) = 49x23, rank=23, k=    13.49; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
  9: dim(M) = 44x16, rank=16, k=    5.585; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 10: dim(M) = 50x19, rank=19, k=    8.837; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 11: dim(M) = 45x10, rank=10, k=    7.917; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 12: dim(M) = 39x17, rank=17, k=    11.81; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 13: dim(M) = 53x18, rank=13, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 14: dim(M) = 45x20, rank=18, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 15: dim(M) = 51x28, rank=27, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 16: dim(M) = 53x19, rank=18, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 17: dim(M) = 42x23, rank=21, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 18: dim(M) = 56x28, rank=28, k=    14.55; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 19: dim(M) = 41x13, rank=12, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 20: dim(M) = 58x17, rank=12, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 21: dim(M) = 55x22, rank=18, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 22: dim(M) = 46x20, rank=17, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 23: dim(M) = 51x18, rank=16, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 24: dim(M) = 46x20, rank=17, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 25: dim(M) = 48x15, rank=11, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 26: dim(M) = 34x23, rank=21, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 27: dim(M) = 63x25, rank=25, k=    18.79; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 28: dim(M) = 36x25, rank=21, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 29: dim(M) = 40x17, rank=11, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 30: dim(M) = 40x19, rank=19, k=    8.338; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 31: dim(M) = 50x19, rank=14, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 32: dim(M) = 50x23, rank=23, k=    11.41; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 33: dim(M) = 49x19, rank=13, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 34: dim(M) = 49x25, rank=25, k=    14.28; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 35: dim(M) = 57x20, rank=20, k=     13.3; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 36: dim(M) = 52x23, rank=21, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 37: dim(M) = 50x21, rank=19, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 38: dim(M) = 48x15, rank=13, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 39: dim(M) = 48x14, rank=13, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 40: dim(M) = 53x22, rank=22, k=    6.927; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 41: dim(M) = 44x15, rank=12, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 42: dim(M) = 59x17, rank=17, k=      4.9; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 43: dim(M) = 49x12, rank=11, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 44: dim(M) = 54x20, rank=20, k=    10.87; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 45: dim(M) = 52x16, rank=15, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 46: dim(M) = 43x18, rank=17, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 47: dim(M) = 64x23, rank=23, k=    20.77; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 48: dim(M) = 51x28, rank=22, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 49: dim(M) = 47x21, rank=21, k=    10.29; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 50: dim(M) = 53x16, rank=16, k=    4.912; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 51: dim(M) = 55x17, rank=15, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 52: dim(M) = 43x10, rank= 8, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 53: dim(M) = 51x14, rank=14, k=    19.61; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 54: dim(M) = 63x22, rank=22, k=    8.823; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 55: dim(M) = 50x11, rank=11, k=     6.33; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 56: dim(M) = 48x17, rank=17, k=    11.92; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 57: dim(M) = 45x15, rank=11, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 58: dim(M) = 46x27, rank=26, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 59: dim(M) = 36x22, rank=18, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 60: dim(M) = 55x22, rank=22, k=    30.21; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 61: dim(M) = 47x24, rank=20, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 62: dim(M) = 55x21, rank=21, k=    7.926; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 63: dim(M) = 59x11, rank=11, k=    28.41; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 64: dim(M) = 63x20, rank=20, k=    10.79; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 65: dim(M) = 41x12, rank=12, k=    4.169; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 66: dim(M) = 66x20, rank=20, k=    13.37; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 67: dim(M) = 58x19, rank=19, k=    7.797; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 68: dim(M) = 60x18, rank=18, k=    5.237; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 69: dim(M) = 50x21, rank=21, k=    14.43; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 70: dim(M) = 39x14, rank=14, k=    6.119; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 71: dim(M) = 46x16, rank=14, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 72: dim(M) = 59x25, rank=25, k=    16.74; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 73: dim(M) = 48x21, rank=20, k=3.136e+18; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 74: dim(M) = 53x16, rank=16, k=    8.899; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 75: dim(M) = 44x18, rank=18, k=    8.835; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 76: dim(M) = 52x19, rank=18, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 77: dim(M) = 46x27, rank=25, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 78: dim(M) = 55x21, rank=17, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 79: dim(M) = 57x22, rank=21, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 80: dim(M) = 48x23, rank=22, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 81: dim(M) = 42x 9, rank= 8, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 82: dim(M) = 41x19, rank=17, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 83: dim(M) = 33x27, rank=23, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 84: dim(M) = 55x25, rank=25, k=    11.44; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 85: dim(M) = 49x17, rank=13, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 86: dim(M) = 55x25, rank=25, k=    237.1; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 87: dim(M) = 62x17, rank=17, k=    4.925; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 88: dim(M) = 48x22, rank=15, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 89: dim(M) = 39x24, rank=20, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 90: dim(M) = 50x19, rank=19, k=    7.178; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 91: dim(M) = 50x20, rank=14, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 92: dim(M) = 39x21, rank=17, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 93: dim(M) = 47x18, rank=18, k=    6.972; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 94: dim(M) = 58x21, rank=21, k=    10.51; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 95: dim(M) = 42x13, rank=13, k=    22.02; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 96: dim(M) = 50x26, rank=22, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 97: dim(M) = 37x15, rank=10, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 98: dim(M) = 49x18, rank=18, k=    16.69; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
 99: dim(M) = 51x22, rank=21, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
100: dim(M) = 46x12, rank= 9, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
101: dim(M) = 47x21, rank=21, k=     12.6; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
102: dim(M) = 51x23, rank=23, k=    12.84; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
103: dim(M) = 44x20, rank=20, k=    7.502; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
104: dim(M) = 47x28, rank=28, k=    16.27; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
105: dim(M) = 40x18, rank=16, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
106: dim(M) = 43x16, rank=13, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
107: dim(M) = 50x20, rank=20, k=    12.23; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
108: dim(M) = 49x18, rank=17, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
109: dim(M) = 36x29, rank=27, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
110: dim(M) = 52x26, rank=25, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
111: dim(M) = 51x18, rank=18, k=    8.149; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
112: dim(M) = 47x28, rank=28, k=    14.75; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
113: dim(M) = 50x26, rank=26, k=    27.23; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
114: dim(M) = 51x25, rank=25, k=    14.32; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
115: dim(M) = 54x22, rank=21, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
116: dim(M) = 54x16, rank=15, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
117: dim(M) = 54x20, rank=18, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
118: dim(M) = 38x21, rank=17, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
119: dim(M) = 34x29, rank=25, k=      Inf; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
120: dim(M) = 48x18, rank=18, k=    7.804; super= FALSE M:  M': super= TRUE M:  M':  [Ok]
Time    user  system elapsed 
Time    1.07    0.00    1.06 
> 
> .updateCHMfactor
function (object, parent, mult = 0) 
.Call(CHMfactor_update, object, parent, mult)
<bytecode: 0x000002a23cd1eba8>
<environment: namespace:Matrix>
> ## TODO: (--> ../TODO "Cholesky"):
> ## ----
> ## allow Cholesky(A,..) when A is not symmetric *AND*
> ## we really want to factorize  AA' ( + beta * I)
> 
> 
> ## Schur() ----------------------
> checkSchur <- function(A, SchurA = Schur(A), tol = 1e-14) {
+     stopifnot(is(SchurA, "Schur"),
+               isOrthogonal(Q <- SchurA@Q),
+               all.equal(as.mat(A),
+                         as.mat(Q %*% SchurA@T %*% t(Q)), tolerance = tol))
+ }
> 
> SH <- Schur(H5 <- Hilbert(5))
> checkSchur(H5, SH)
> checkSchur(Diagonal(x = 9:3))
> 
> p <- 4L
> uTp <- new("dtpMatrix", x=c(2, 3, -1, 4:6, -2:1), Dim = c(p,p))
> (uT <- as(uTp, "unpackedMatrix"))
4 x 4 Matrix of class "dtrMatrix"
     [,1] [,2] [,3] [,4]
[1,]    2    3    4   -2
[2,]    .   -1    5   -1
[3,]    .    .    6    0
[4,]    .    .    .    1
> ## Schur ( <general> )  <--> Schur( <triangular> )
> Su <- Schur(uT) ;   checkSchur(uT, Su)
> gT <- as(uT,"generalMatrix")
> Sg <- Schur(gT) ;   checkSchur(gT, Sg)
> Stg <- Schur(t(gT));checkSchur(t(gT), Stg)
> Stu <- Schur(t(uT));checkSchur(t(uT), Stu)
> 
> stopifnot(identical3(Sg@T, uT, Su@T),
+           identical(Sg@Q, as(diag(p), "generalMatrix")),
+           identical(Stg@T, as(t(gT[,p:1])[,p:1], "triangularMatrix")),
+           identical(Stg@Q, as(diag(p)[,p:1], "generalMatrix")),
+           identical(Stu@T, Stg@T))
> assert.EQ.mat(Stu@Q, as(Stg@Q,"matrix"), tol=0)
> 
> ## the pedigreemm example where solve(.) failed:
> p <- new("dtCMatrix", i = c(2L, 3L, 2L, 5L, 4L, 4:5), p = c(0L, 2L, 4:7, 7L),
+ 	 Dim = c(6L, 6L), Dimnames = list(as.character(1:6), NULL),
+ 	 x = rep.int(-0.5, 7), uplo = "L", diag = "U")
> Sp <- Schur(p)
> Sp. <- Schur(as(p,"generalMatrix"))
> Sp.p <- Schur(crossprod(p))
> ## the last two failed
> ip <- solve(p)
> assert.EQ.mat(solve(ip), as(p,"matrix"))
> 
> 
> ## chol2inv() for a traditional matrix
> assert.EQ.mat(     crossprod(chol2inv(chol(Diagonal(x = 5:1)))),
+               C <- crossprod(chol2inv(chol(    diag(x = 5:1)))))
> stopifnot(all.equal(C, diag((5:1)^-2)))
> ## failed in some versions because of a "wrong" implicit generic
> 
> U <- cbind(1:0, 2*(1:2))
> (sU <- as(U, "CsparseMatrix"))
2 x 2 sparse Matrix of class "dtCMatrix"
        
[1,] 1 2
[2,] . 4
> validObject(sS <- crossprod(sU))
[1] TRUE
> C. <- chol(sS)
> stopifnot(all.equal(C., sU, tolerance=1e-15))
> ## chol(<triangular sparse which is diagonal>)
> tC7 <- .trDiagonal(7, 7:1)
> stopifnotValid(tC7, "dtCMatrix")
> ch7  <- chol(tC7) ## this (and the next 2) failed: 'no slot .. "factors" ..."dtCMatrix"'
> chT7 <- chol(tT7 <- as(tC7, "TsparseMatrix"))
> chR7 <- chol(tR7 <- as(tC7, "RsparseMatrix"))
> stopifnot(expr = {
+     isDiagonal(ch7)
+     identical(chT7, ch7) # "ddiMatrix" all of them
+     identical(chR7, ch7) # "ddiMatrix" all of them
+     all.equal(sqrt(7:1), diag(ch7 ))
+ })
> 
> 
> 
> ## From [Bug 14834] New: chol2inv *** caught segfault ***
> n <- 1e6 # was 595362
> A <- chol( D <- Diagonal(n) )
> stopifnot(identical(A,D)) # A remains (unit)diagonal
> is(tA <- as(A,"triangularMatrix"))# currently a dtTMatrix
[1] "dtCMatrix"        "CsparseMatrix"    "dsparseMatrix"    "triangularMatrix"
[5] "dMatrix"          "sparseMatrix"     "Matrix"           "replValueSp"     
> stopifnotValid(tA, "dsparseMatrix")
> CA <- as(tA, "CsparseMatrix")
> 
> selectMethod(solve, c("dtCMatrix","missing"))
Method Definition:

function (a, b, ...) 
{
    .local <- function (a, b, sparse = TRUE, ...) 
    {
        if (a@diag != "N") 
            a <- ..diagU2N(a)
        .Call(dtCMatrix_solve, a, NULL, sparse)
    }
    .local(a, b = b, ...)
}
<bytecode: 0x000002a235d175a0>
<environment: namespace:Matrix>

Signatures:
        a           b        
target  "dtCMatrix" "missing"
defined "dtCMatrix" "missing"
> ##--> .Call(dtCMatrix_sparse_solve, a, .trDiagonal(n))  in ../src/dtCMatrix.c
> sA  <- solve(CA)## -- R_CheckStack() segfault in Matrix <= 1.0-4
> nca <- diagU2N(CA)
> stopifnot(identical(sA, nca))
> ## same check with non-unit-diagonal D :
> A <- chol(D <- Diagonal(n, x = 0.5))
> ia <- chol2inv(A)
> stopifnot(is(ia, "diagonalMatrix"),
+ 	  all.equal(ia@x, rep(2,n), tolerance = 1e-15))
> 
> ##------- Factor caches must be cleaned - even after scalar-Ops such as "2 *"
> set.seed(7)
> d <- 5
> S <- 10*Diagonal(d) + rsparsematrix(d,d, 1/4)
> class(M <- as(S, "denseMatrix")) # dgeMatrix
[1] "dgeMatrix"
attr(,"package")
[1] "Matrix"
> m <- as.matrix(M)
> (dS <- determinant(S))
$modulus
[1] 11.47196
attr(,"logarithm")
[1] TRUE

$sign
[1] 1

attr(,"class")
[1] "det"
> stopifnot(exprs = {
+     all.equal(determinant(m), dS, tolerance=1e-15)
+     all.equal(dS, determinant(M), tolerance=1e-15)
+     ## These had failed, as the "LU" factor cache was kept unchanged in 2*M :
+     all.equal(determinant(2*S), determinant(2*M) -> d2M)
+     all.equal(determinant(S^2), determinant(M^2) -> dM2)
+     all.equal(determinant(m^2), dM2)
+     all.equal(d*log(2), c(d2M$modulus - dS$modulus))
+ })
> 
> ## misc. bugs found in Matrix 1.4-1
> L. <- new("dtCMatrix", Dim = c(1L, 1L), uplo = "L",
+           p = c(0L, 1L), i = 0L, x = 1)
> S. <- forceSymmetric(L.)
> lu(S.)
LU factorization of Formal class 'sparseLU' [package "Matrix"] with 6 slots
  ..@ L       :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
  .. .. ..@ i       : int 0
  .. .. ..@ p       : int [1:2] 0 1
  .. .. ..@ Dim     : int [1:2] 1 1
  .. .. ..@ Dimnames:List of 2
  .. .. .. ..$ : NULL
  .. .. .. ..$ : NULL
  .. .. ..@ x       : num 1
  .. .. ..@ uplo    : chr "L"
  .. .. ..@ diag    : chr "N"
  ..@ U       :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
  .. .. ..@ i       : int 0
  .. .. ..@ p       : int [1:2] 0 1
  .. .. ..@ Dim     : int [1:2] 1 1
  .. .. ..@ Dimnames:List of 2
  .. .. .. ..$ : NULL
  .. .. .. ..$ : NULL
  .. .. ..@ x       : num 1
  .. .. ..@ uplo    : chr "U"
  .. .. ..@ diag    : chr "N"
  ..@ p       : int 0
  ..@ q       : int 0
  ..@ Dim     : int [1:2] 1 1
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
> stopifnot(validObject(lu(L.)), # was invalid
+           identical(names(S.@factors), "sparseLU")) # was "lu"
> 
> ## chol() should give matrix with 'Dimnames',
> ## even if 'Dimnames' are not cached
> D. <- as(diag(3), "CsparseMatrix")
> D.@Dimnames <- dn <- list(zzz = letters[1:3], ZZZ = LETTERS[1:3])
> cd1 <- chol(D.) # "fresh"
> stopifnot(identical(cd1@Dimnames, rep(dn[2L], 2L)))
> cd2 <- chol(D.) # from cache
> stopifnot(identical(cd1, cd2))
> 
> ## lu(<m-by-0>), lu(<0-by-n>), BunchKaufman(<0-by-0>), chol(<0-by-0>)
> stopifnot(identical(lu(new("dgeMatrix", Dim = c(2L, 0L))),
+                     new("denseLU", Dim = c(2L, 0L))),
+           identical(lu(new("dgeMatrix", Dim = c(0L, 2L))),
+                     new("denseLU", Dim = c(0L, 2L))),
+           identical(BunchKaufman(new("dsyMatrix", uplo = "U")),
+                     new("BunchKaufman", uplo = "U")),
+           identical(BunchKaufman(new("dspMatrix", uplo = "L")),
+                     new("pBunchKaufman", uplo = "L")),
+           identical(Cholesky(new("dpoMatrix", uplo = "U")),
+                     new("Cholesky", uplo = "U")),
+           identical(Cholesky(new("dppMatrix", uplo = "L")),
+                     new("pCholesky", uplo = "L")))
> 
> ## determinant(<ds[yp]Matrix>) going via Bunch-Kaufman
> set.seed(15742)
> n <- 10L
> syU <- syL <- new("dsyMatrix", Dim = c(n, n), x = rnorm(n * n))
> spU <- spL <- new("dspMatrix", Dim = c(n, n), x = rnorm((n * (n + 1L)) %/% 2L))
> syL@uplo <- spL@uplo <- "L"
> for(m in list(syU, syL, spU, spL))
+     for(givelog in c(FALSE, TRUE))
+         stopifnot(all.equal(determinant(   m,            givelog),
+                             determinant(as(m, "matrix"), givelog)))
> 
> ## was an error at least in Matrix 1.5-4 ...
> BunchKaufman(as.matrix(1))
Bunch-Kaufman factorization of Formal class 'BunchKaufman' [package "Matrix"] with 5 slots
  ..@ uplo    : chr "U"
  ..@ x       : num 1
  ..@ perm    : int 1
  ..@ Dim     : int [1:2] 1 1
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
> 
> 
> ## 'expand2': product of listed factors should reproduce factorized matrix
> ## FIXME: many of our %*% methods still mangle dimnames or names(dimnames) ...
> ##        hence for now we coerce the factors to matrix before multiplying
> chkMF <- function(X, Y, FUN, ...) {
+     ## t(x)@factors may preserve factorizations with x@uplo
+     X@factors <- list()
+ 
+     mf <- FUN(X, ...)
+     e2.mf <- expand2(mf)
+     e1.mf <- sapply(names(e2.mf), expand1, x = mf, simplify = FALSE)
+ 
+     m.e2.mf <- lapply(e2.mf, as, "matrix")
+     m.e1.mf <- lapply(e1.mf, as, "matrix")
+ 
+     identical(m.e1.mf, lapply(m.e2.mf, unname)) &&
+         isTRUE(all.equal(Reduce(`%*%`, m.e2.mf), Y))
+ }
> set.seed(24831)
> n <- 16L
> mS <- tcrossprod(matrix(rnorm(n * n), n, n,
+                         dimnames = list(A = paste0("s", seq_len(n)), NULL)))
> sS <- as(pS <- as(S <- as(mS, "dpoMatrix"), "packedMatrix"), "CsparseMatrix")
> stopifnot(exprs = {
+     chkMF(   S , mS,    Schur)
+     chkMF(  pS , mS,    Schur)
+     chkMF(   S , mS,       lu)
+     chkMF(  pS , mS,       lu)
+     chkMF(  sS , mS,       lu)
+     chkMF(  sS , mS,       qr)
+     chkMF(   S , mS, BunchKaufman)
+     chkMF(  pS , mS, BunchKaufman)
+     chkMF(t( S), mS, BunchKaufman)
+     chkMF(t(pS), mS, BunchKaufman)
+     chkMF(   S , mS, Cholesky)
+     chkMF(  pS , mS, Cholesky)
+     chkMF(t( S), mS, Cholesky)
+     chkMF(t(pS), mS, Cholesky)
+     chkMF(  sS , mS, Cholesky, super = FALSE, LDL =  TRUE)
+     chkMF(  sS , mS, Cholesky, super = FALSE, LDL = FALSE)
+     chkMF(  sS , mS, Cholesky, super =  TRUE, LDL = FALSE)
+ })
> 
> 
> cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
Time elapsed:  5.84 0.18 6.01 NA NA 
> 
> proc.time()
   user  system elapsed 
   5.84    0.18    6.01