Running 'Class+Meth.R' [3s]
  Running 'Simple.R' [10s]
  Running 'abIndex-tsts.R' [4s]
  Running 'base-matrix-fun.R' [3s]
  Running 'bind.R' [4s]
  Comparing 'bind.Rout' to 'bind.Rout.save' ... OK
  Running 'dg_Matrix.R' [3s]
  Running 'dpo-test.R' [5s]
  Running 'dtpMatrix.R' [3s]
  Running 'factorizing.R' [18s]
  Running 'group-methods.R' [5s]
  Running 'indexing.R' [12s]
  Running 'matprod.R' [6s]
  Running 'matr-exp.R' [3s]
  Running 'other-pkgs.R' [2s]
  Running 'spModel.matrix.R' [4s]
  Running 'validObj.R' [6s]
  Running 'write-read.R' [3s]
Running the tests in 'tests/dpo-test.R' failed.
Complete output:
  > ### Testing positive definite matrices
  > library(Matrix)
  > source(system.file("test-tools.R", package = "Matrix"))# identical3() etc
  Loading required package: tools
  > cat("doExtras:",doExtras,"\n")
  doExtras: FALSE 
  > 
  > h9 <- Hilbert(9)
  > stopifnot(c(0,0) == dim(Hilbert(0)),
  +           c(9,9) == dim(h9),
  + 	  identical(h9@factors, list()))
  > str(h9)# no 'factors'		32b:	-96.73694669	2.08e-8
  Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
    ..@ x       : num [1:81] 1 0.5 0.333 0.25 0.2 ...
    ..@ Dim     : int [1:2] 9 9
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ uplo    : chr "U"
    ..@ factors : list()
  > assert.EQ.(c(determinant(h9)$modulus),	-96.7369487, tol = 8e-8)
  Mean relative difference: 2.08569e-08
  > ##				64b:	-96.73695078	2.15e-8 then 6.469e-8
  > 
  > ## determinant() now working via chol(): ==> h9 now has factorization
  > stopifnot(names(h9@factors) == "Cholesky",
  +           identical(ch9 <- chol(h9), h9@factors$Cholesky))
  > round(ch9, 3) ## round() preserves 'triangular' !
  9 x 9 Matrix of class "Cholesky"
        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] 
   [1,] 1.000 0.500 0.333 0.250 0.200 0.167 0.143 0.125 0.111
   [2,]     . 0.289 0.289 0.260 0.231 0.206 0.186 0.168 0.154
   [3,]     .     . 0.075 0.112 0.128 0.133 0.133 0.130 0.126
   [4,]     .     .     . 0.019 0.038 0.052 0.063 0.070 0.075
   [5,]     .     .     .     . 0.005 0.012 0.019 0.027 0.033
   [6,]     .     .     .     .     . 0.001 0.004 0.007 0.010
   [7,]     .     .     .     .     .     . 0.000 0.001 0.002
   [8,]     .     .     .     .     .     .     . 0.000 0.000
   [9,]     .     .     .     .     .     .     .     . 0.000
  > str(f9 <- as(ch9, "dtrMatrix"))
  Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
    ..@ x       : num [1:81] 1 0 0 0 0 0 0 0 0 0.5 ...
    ..@ Dim     : int [1:2] 9 9
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ uplo    : chr "U"
    ..@ diag    : chr "N"
  > stopifnot(all.equal(rcond(h9), 9.0938e-13),
  +           all.equal(rcond(f9), 9.1272e-7, tolerance = 1e-6))# more precision fails
  > options(digits=4)
  > (cf9 <- crossprod(f9))# looks the same as  h9 :
  9 x 9 Matrix of class "dpoMatrix"
          [,1]   [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
   [1,] 1.0000 0.5000 0.33333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111
   [2,] 0.5000 0.3333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000
   [3,] 0.3333 0.2500 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091
   [4,] 0.2500 0.2000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333
   [5,] 0.2000 0.1667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692
   [6,] 0.1667 0.1429 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143
   [7,] 0.1429 0.1250 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667
   [8,] 0.1250 0.1111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250
   [9,] 0.1111 0.1000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250 0.05882
  > assert.EQ.mat(h9, as(cf9,"matrix"), tol=1e-15)
  > 
  > h9. <- round(h9, 2)# actually loses pos.def. "slightly"
  >                    # ==> the above may be invalid in the future
  > h9p  <- as(h9,  "dppMatrix")
  > h9.p <- as(h9., "dppMatrix")
  > ch9p <- chol(h9p)
  > stopifnot(identical(ch9p, h9p@factors$pCholesky),
  + 	  identical(names(h9p@factors), c("Cholesky", "pCholesky")))
  > h4  <- h9.[1:4, 1:4] # this and the next
  > h9.[1,1] <- 10       # had failed in 0.995-14
  > h9p[1,1] <- 10
  > stopifnotValid(h9., "symmetricMatrix")
  > stopifnotValid(h9p, "symmetricMatrix")
  > stopifnotValid(h4,  "symmetricMatrix")
  > 
  > h9p[1,2] <- 99
  > stopifnot(class(h9p) == "dgeMatrix", h9p[1,1:2] == c(10,99))
  > 
  > str(h9p <- as(h9, "dppMatrix"))# {again}
  Formal class 'dppMatrix' [package "Matrix"] with 5 slots
    ..@ x       : num [1:45] 1 0.5 0.333 0.333 0.25 ...
    ..@ Dim     : int [1:2] 9 9
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ uplo    : chr "U"
    ..@ factors :List of 1
    .. ..$ Cholesky:Formal class 'Cholesky' [package "Matrix"] with 5 slots
    .. .. .. ..@ x       : num [1:81] 1 0 0 0 0 0 0 0 0 0.5 ...
    .. .. .. ..@ Dim     : int [1:2] 9 9
    .. .. .. ..@ Dimnames:List of 2
    .. .. .. .. ..$ : NULL
    .. .. .. .. ..$ : NULL
    .. .. .. ..@ uplo    : chr "U"
    .. .. .. ..@ diag    : chr "N"
  > h6 <- h9[1:6,1:6]
  > stopifnot(all(h6 == Hilbert(6)), length(h6@factors) == 0)
  > stopifnotValid(th9p <- t(h9p), "dppMatrix")
  > stopifnotValid(h9p@factors$Cholesky,"Cholesky")
  > H6  <- as(h6, "dspMatrix")
  > pp6 <- as(H6, "dppMatrix")
  > po6 <- as(pp6,"dpoMatrix")
  > hs <- as(h9p, "dspMatrix")
  > stopifnot(names(H6@factors)  == "pCholesky",
  + 	  names(pp6@factors) == "pCholesky",
  + 	  names(hs@factors)  == "Cholesky") # for now
  > chol(hs) # and that is cached in 'hs' too :
  'MatrixFactorization' of Formal class 'pCholesky' [package "Matrix"] with 5 slots
    ..@ x       : num [1:45] 1 0.5 0.289 0.333 0.289 ...
    ..@ Dim     : int [1:2] 9 9
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ uplo    : chr "U"
    ..@ diag    : chr "N"
  > stopifnot(names(hs@factors) %in% c("Cholesky","pCholesky"),
  + 	  all.equal(h9, crossprod(hs@factors$pCholesky), tolerance =1e-13),
  + 	  all.equal(h9, crossprod(hs@factors$ Cholesky), tolerance =1e-13))
  > 
  > hs@x <- 1/h9p@x # is not pos.def. anymore
  > validObject(hs) # "but" this does not check
  [1] TRUE
  > stopifnot(diag(hs) == seq(1, by = 2, length = 9))
  > 
  > s9 <- solve(h9p, seq(nrow(h9p)))
  > signif(t(s9)/10000, 4)# only rounded numbers are platform-independent
  1 x 9 Matrix of class "dgeMatrix"
         [,1]  [,2]  [,3]   [,4] [,5]  [,6] [,7]  [,8] [,9]
  [1,] 0.0729 -5.76 109.5 -864.9 3468 -7668 9459 -6095 1597
  > (I9 <- h9p %*% s9)
  9 x 1 Matrix of class "dgeMatrix"
        [,1]
   [1,]    1
   [2,]    2
   [3,]    3
   [4,]    4
   [5,]    5
   [6,]    6
   [7,]    7
   [8,]    8
   [9,]    9
  > m9 <- matrix(1:9, dimnames = list(NULL,NULL))
  > stopifnot(all.equal(m9, .asmatrix(I9), tolerance = 2e-9))
  > 
  > ### Testing nearPD() --- this is partly in  ../man/nearPD.Rd :
  > pr <- Matrix(c(1,     0.477, 0.644, 0.478, 0.651, 0.826,
  +                0.477, 1,     0.516, 0.233, 0.682, 0.75,
  +                0.644, 0.516, 1,     0.599, 0.581, 0.742,
  +                0.478, 0.233, 0.599, 1,     0.741, 0.8,
  +                0.651, 0.682, 0.581, 0.741, 1,     0.798,
  +                0.826, 0.75,  0.742, 0.8,   0.798, 1),
  +              nrow = 6, ncol = 6)
  > 
  > nL <-
  +     list(r   = nearPD(pr, conv.tol = 1e-7), # default
  + 	 r.1 = nearPD(pr, conv.tol = 1e-7,		corr = TRUE),
  + 	 rs  = nearPD(pr, conv.tol = 1e-7, doDyk=FALSE),
  + 	 rs1 = nearPD(pr, conv.tol = 1e-7, doDyk=FALSE, corr = TRUE),
  + 	 rH  = nearPD(pr, conv.tol = 1e-15),
  +          rH.1= nearPD(pr, conv.tol = 1e-15, corr = TRUE))
  > 
  > sapply(nL, `[`, c("iterations", "normF"))
             r       r.1     rs      rs1     rH      rH.1   
  iterations 2       12      2       11      2       30     
  normF      0.06262 0.07429 0.06262 0.07438 0.06262 0.07429
  > 
  > allnorms <- function(d) vapply(c("1","I","F","M","2"),
  +                                norm, x = d, double(1))
  > 
  > ## "F" and "M" distances are larger for the (corr=TRUE) constrained:
  > 100 * sapply(nL, function(rr) allnorms((pr - rr  $ mat)))
        r   r.1    rs   rs1    rH  rH.1
  1 8.889 8.805 8.889 8.746 8.889 8.805
  I 8.889 8.805 8.889 8.746 8.889 8.805
  F 6.262 7.429 6.262 7.438 6.262 7.429
  M 2.655 2.860 2.655 2.785 2.655 2.860
  2 6.262 6.418 6.262 6.461 6.262 6.418
  > 
  > ## But indeed, the 'corr = TRUE' constraint yield a better solution,
  > ## if you need the constraint :  cov2cor() does not just fix it up :
  > 100 * (nn <- sapply(nL, function(rr) allnorms((pr - cov2cor(rr  $ mat)))))
        r   r.1    rs   rs1    rH  rH.1
  1 9.994 8.805 9.994 8.746 9.994 8.805
  I 9.994 8.805 9.994 8.746 9.994 8.805
  F 8.391 7.429 8.391 7.438 8.391 7.429
  M 3.544 2.860 3.544 2.785 3.544 2.860
  2 6.322 6.418 6.322 6.461 6.322 6.418
  > 
  > stopifnot(
  + all.equal(nn["1",],
  +           c(r =0.0999444286984696, r.1= 0.0880468666522317,
  +             rs=0.0999444286984702, rs1= 0.0874614179943388,
  +             rH=0.0999444286984696, rH.1=0.0880468927726625),
  +           tolerance=1e-9))
  > 
  > nr <- nL $rH.1 $mat
  > stopifnot(
  +     all.equal(nr[lower.tri(nr)],
  + 	      c(0.4877861230299, 0.6429309061748, 0.4904554299278, 0.6447150779852,
  + 		0.8082100656035, 0.514511537243, 0.2503412693503, 0.673249718642,
  + 		0.7252316891977, 0.5972811755863, 0.5818673040157, 0.7444549621769,
  + 		0.7308954865819, 0.7713984381710, 0.8124321235679),
  + 	      tolerance = 1e-9))
  > showProc.time()
  Time elapsed:  1.06 0.03 1.16 
  > 
  > 
  > set.seed(27)
  > m9 <- h9 + rnorm(9^2)/1000 ; m9 <- (m9 + t(m9))/2
  > nm9 <- nearPD(m9)
  > showProc.time()
  Time elapsed:  0.06 0 0.06 
  > 
  > nRep <- if(doExtras) 50 else 4
  > CPU <- 0
  > for(M in c(5, 12))
  +     for(i in 1:nRep) {
  + 	m <- matrix(round(rnorm(M^2),2), M, M)
  + 	m <- m + t(m)
  + 	diag(m) <- pmax(0, diag(m)) + 1
  + 	m <- cov2cor(m)
  + 	CPU <- CPU + system.time(n.m <- nearPD(m))[1]
  + 	X <- as(n.m$mat, "matrix")
  + 	stopifnot(all.equal(X, (X + t(X))/2, tolerance = 8*.Machine$double.eps),
  + 		  all.equal(eigen(n.m$mat, only.values=TRUE)$values,
  + 			    n.m$eigenvalues, tolerance = 4e-8))
  +     }
  > cat('Time elapsed for ',nRep, 'nearPD(): ', CPU,'\n')
  Time elapsed for  4 nearPD():  0.01 
  > showProc.time()
  Time elapsed:  1.48 0 1.48 
  > 
  > ## cov2cor()
  > m <- diag(6:1) %*% as(pr,"matrix") %*% diag(6:1) # so we can "vector-index"
  > m[upper.tri(m)] <- 0
  > ltm <- which(lower.tri(m))
  > ne <- length(ltm)
  > set.seed(17)
  > m[ltm[sample(ne, 3/4*ne)]] <- 0
  > m <- (m + t(m))/2 # now is a covariance matrix with many 0 entries
  > (spr <- Matrix(m))
  6 x 6 sparse Matrix of class "dsCMatrix"
                                     
  [1,] 36.000  .      . 4.302 . 2.478
  [2,]  .     25.000  . .     . 1.875
  [3,]  .      .     16 .     . .    
  [4,]  4.302  .      . 9.000 . 1.200
  [5,]  .      .      . .     4 .    
  [6,]  2.478  1.875  . 1.200 . 1.000
  > cspr <- cov2cor(spr)
  > ev <- eigen(cspr, only.v = TRUE)$values
  > stopifnot(is(spr, "dsCMatrix"),
  +           is(cspr,"dsCMatrix"),
  +           all.equal(ev, c(1.5901626099,  1.1902658504, 1, 1,
  +                           0.80973414959, 0.40983739006), tolerance=1e-10))
  Error: ev and c(1.5901626099, 1.1902658504, 1, 1, 0.80973414959, 0.40983739006) are not equal:
    Mean relative difference: 0.09731
  Execution halted
Running the tests in 'tests/indexing.R' failed.
Complete output:
  > #### For both 'Extract' ("[") and 'Replace' ("[<-") Method testing
  > ####    aka    subsetting     and  subassignment
  > 
  > if(interactive()) {
  +     options(error = recover, warn = 1)
  + } else if(FALSE) { ## MM @ testing *manually* only
  +     options(error = recover, Matrix.verbose = TRUE, warn = 1)
  + } else {
  +     options(                 Matrix.verbose = TRUE, warn = 1)
  + }
  > ## Matrix.verbose = TRUE (*before* loading 'Matrix' pkg)
  > ## ==> will also show method dispath ambiguity messages: getOption("ambiguousMethodSelection")
  > 
  > #### suppressPackageStartupMessages(...)  as we have an *.Rout.save to Rdiff against
  > stopifnot(suppressPackageStartupMessages(require(Matrix)))
  > 
  > source(system.file("test-tools.R", package = "Matrix"), keep.source = FALSE)
  Loading required package: tools
  > ##-> identical3() etc
  > cat("doExtras:",doExtras,"\n")
  doExtras: FALSE 
  > 
  > 
  > ### Dense Matrices
  > 
  > m <- Matrix(1:28 +0, nrow = 7)
  > validObject(m)
  [1] TRUE
  > stopifnot(identical(m, m[]),
  + 	  identical(m[2, 3],  16), # simple number
  + 	  identical(m[2, 3:4], c(16,23)), # simple numeric of length 2
  + 	  identical(m[NA,NA], as(Matrix(NA, 7,4), "dMatrix")))
  > 
  > m[2, 3:4, drop=FALSE] # sub matrix of class 'dgeMatrix'
  1 x 2 Matrix of class "dgeMatrix"
       [,1] [,2]
  [1,]   16   23
  > m[-(4:7), 3:4]        # ditto; the upper right corner of 'm'
  3 x 2 Matrix of class "dgeMatrix"
       [,1] [,2]
  [1,]   15   22
  [2,]   16   23
  [3,]   17   24
  > 
  > ## rows or columns only:
  > m[1,]     # first row, as simple numeric vector
  [1]  1  8 15 22
  > m[,2]     # 2nd column
  [1]  8  9 10 11 12 13 14
  > m[,1:2]   # sub matrix of first two columns
  7 x 2 Matrix of class "dgeMatrix"
       [,1] [,2]
  [1,]    1    8
  [2,]    2    9
  [3,]    3   10
  [4,]    4   11
  [5,]    5   12
  [6,]    6   13
  [7,]    7   14
  > m[-(1:6),, drop=FALSE] # not the first 6 rows, i.e. only the 7th
  1 x 4 Matrix of class "dgeMatrix"
       [,1] [,2] [,3] [,4]
  [1,]    7   14   21   28
  > m[integer(0),] #-> 0 x 4 Matrix
  0 x 4 Matrix of class "dgeMatrix"
       [,1] [,2] [,3] [,4]
  > m[2:4, numeric(0)] #-> 3 x 0 Matrix
  3 x 0 Matrix of class "dgeMatrix"
      
  [1,]
  [2,]
  [3,]
  > 
  > ## logical indexing
  > stopifnot(identical(m[2,3], m[(1:nrow(m)) == 2, (1:ncol(m)) == 3]),
  +           identical(m[2,], m[(1:nrow(m)) == 2, ]),
  +           identical(m[,3:4], m[, (1:4) >= 3]))
  > 
  > ## dimnames indexing:
  > mn <- m
  > dimnames(mn) <- list(paste("r",letters[1:nrow(mn)],sep=""),
  +                      LETTERS[1:ncol(mn)])
  > checkMatrix(mn)
   norm(m [7 x 4]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  > mn["rd", "D"]
  [1] 25
  > msr <- ms <- as(mn,"sparseMatrix")
  > mnr <- mn
  > v <- rev(as(ms, "vector"))
  > mnr[] <- v
  > msr[] <- v # [<- "sparse" -- not very sensical; did fail w/o a message
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  > z <- msr; z[] <- 0
  > zz <- as(array(0, dim(z)), "sparseMatrix")
  > a.m <- as(mnr,"matrix")
  > stopifnot(identical(mn["rc", "D"], mn[3,4]), mn[3,4] == 24,
  + 	  identical(mn[, "A"], mn[,1]), mn[,1] == 1:7,
  + 	  identical(mn[c("re", "rb"), "B"], mn[c(5,2), 2]),
  + 	  identical(ms["rc", "D"], ms[3,4]), ms[3,4] == 24,
  + 	  identical(ms[, "A"], ms[,1]), ms[,1] == 1:7,
  + 	  identical(ms[ci <- c("re", "rb"), "B"], ms[c(5,2), 2]),
  + 	  identical(rownames(mn[ci, ]), ci),
  + 	  identical(rownames(ms[ci, ]), ci),
  + 	  identical(colnames(mn[,cj <- c("B","D")]), cj),
  + 	  identical(colnames(ms[,cj]), cj),
  + 	  identical(a.m, as(msr,"matrix")),
  + 	  identical(unname(z), zz),
  + 	  identical(a.m, array(v, dim=dim(mn), dimnames=dimnames(mn)))
  + 	  )
  > showProc.time()
  Time elapsed:  1.08 0.06 1.17 
  > 
  > ## Bug found thanks to Timothy Mak, Feb 3, 2017:
  > ## sparseMatrix logical indexing with (partial) NA:
  > a.m <- as(mn,"matrix")
  > assert.EQ(as(ms,"matrix"), a.m) # incl. dimnames
  > iN4 <- c(NA, TRUE, FALSE, TRUE)
  > assert.EQ(as(mn[,iN4],"matrix"), a.m[,iN4]) # (incl. dimnames)
  > ##assert.EQ(as.matrix(ms[,iN4]), a.m[,iN4]) # ms[, <with_NA>]  fails still :
  > try(ms[,iN4])
  Error in intI(j, n = x@Dim[2], dn[[2]], give.dn = FALSE) : 
    'NA' indices are not (yet?) supported for sparse Matrices
  > try(ms[,iN4] <- 100) ## <- segfaulted in Matrix <= 1.2-8  (!)
  Error in intI(i, n = di[margin], dn = dn[[margin]], give.dn = FALSE) : 
    'NA' indices are not (yet?) supported for sparse Matrices
  > 
  > ## R-forge Matrix bug #2556: Subsetting a sparse matrix did remove  names(dimnames(.)) :
  > m44 <- matrix(1:16, 4, 4, dimnames=list(row=c('a','b','c','d'), col=c('x','y','z','w')))
  > ## Dense matrix: ------------------------------------------
  > a <- Matrix(m44)
  > identical(
  +     dimnames(m44[,FALSE, drop=FALSE]),
  +     dimnames(  a[,FALSE, drop=FALSE]))
  [1] TRUE
  > chk.ndn <- function(a, a0=m44)
  +     stopifnot(identical(names(dimnames(a)), names(dimnames(a0))))
  > i <- 1:2
  > chk.ndn(a[i,]); chk.ndn(a[i, i])
  > ## Sparse matrix: -----------------------------------------
  > s <- as(a %% 3 == 1, "sparseMatrix")
  > ts <- as(s,"TsparseMatrix")
  > b <- sparseMatrix(i=1:3, j=rep(2,3), dims=c(4,4), dimnames=dimnames(s))
  > tb <- as(b,"TsparseMatrix")
  > stopifnot(identical5(
  +     dimnames(a), dimnames(s), dimnames(ts),
  +     dimnames(b), dimnames(tb)))
  > 
  > chk.ndn(b [i, i]); chk.ndn(b [i, ])
  > chk.ndn(s [i, i]); chk.ndn(s [i, ])
  > chk.ndn(tb[i, i]); chk.ndn(tb[i, ])
  > chk.ndn(ts[i, i]); chk.ndn(ts[i, ])
  > chk.ndn( b[ , 1, drop=FALSE]); chk.ndn( s[i, 2, drop=FALSE])
  > chk.ndn(tb[ , 1, drop=FALSE]); chk.ndn(ts[i, 2, drop=FALSE])
  > 
  > ## Printing sparse colnames:
  > ms[sample(28, 20)] <- 0
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  > ms <- t(rbind2(ms, 3*ms))
  > cnam1 <- capture.output(show(ms))[2] ; op <- options("sparse.colnames" = "abb3")
     [[ suppressing 14 column names 'ra', 'rb', 'rc' ... ]]
  > cnam2 <- capture.output(show(ms))[2] ; options(op) # revert
  > stopifnot(## sparse printing
  + 	  grep("^ +$", cnam1) == 1, # cnam1 is empty
  + 	  identical(cnam2,
  + 		    paste(" ", paste(rep(rownames(mn), 2), collapse=" "))))
  > 
  > mo <- m
  > m[2,3] <- 100
  > m[1:2, 4] <- 200
  > m[, 1] <- -1
  > m[1:3,]
  3 x 4 Matrix of class "dgeMatrix"
       [,1] [,2] [,3] [,4]
  [1,]   -1    8   15  200
  [2,]   -1    9  100  200
  [3,]   -1   10   17   24
  > 
  > m. <- .asmatrix(m)
  > 
  > ## m[ cbind(i,j) ] indexing:
  > iN <- ij <- cbind(1:6, 2:3)
  > iN[2:3,] <- iN[5,2] <- NA
  > stopifnot(identical(m[ij], m.[ij]),
  + 	  identical(m[iN], m.[iN]))
  > 
  > ## testing operations on logical Matrices rather more than indexing:
  > g10 <- m [ m > 10 ]
  > stopifnot(18 == length(g10))
  > stopifnot(10 == length(m[ m <= 10 ]))
  > sel <- (20 <  m) & (m <  150)
  > sel.<- (20 <  m.)& (m.<  150)
  > nsel <-(20 >= m) | (m >= 150)
  > (ssel <- as(sel, "sparseMatrix"))
  7 x 4 sparse Matrix of class "lgCMatrix"
              
  [1,] . . . .
  [2,] . . | .
  [3,] . . . |
  [4,] . . . |
  [5,] . . . |
  [6,] . . . |
  [7,] . . | |
  > stopifnot(is(sel, "lMatrix"), is(ssel, "lsparseMatrix"),
  + 	  identical3(as.mat(sel.), as.mat(sel), as.mat(ssel)),
  + 	  identical3(!sel, !ssel, nsel), # !<sparse> is typically dense
  + 	  identical3(m[ sel],  m[ ssel], .asmatrix(m)[.asmatrix( ssel)]),
  + 	  identical3(m[!sel],  m[!ssel], .asmatrix(m)[.asmatrix(!ssel)])
  + 	  )
  > showProc.time()
  Time elapsed:  0.45 0.01 0.48 
  > 
  > ## more sparse Matrices --------------------------------------
  > 
  > ##' @title Check sparseMatrix sub-assignment   m[i,j] <- v
  > ##' @param ms sparse Matrix
  > ##' @param mm its [traditional matrix]-equivalent
  > ##' @param k (approximate) length of index vectors (i,j)
  > ##' @param n.uniq (approximate) number of unique values in  i,j
  > ##' @param vRNG function(n) for random 'v' generation
  > ##' @param show logical; if TRUE, it will not stop on error
  > ##' @return
  > ##' @author Martin Maechler
  > chkAssign <- function(ms, mm = as(ms, "matrix"),
  +                       k = min(20,dim(mm)), n.uniq = k %/% 3,
  +                       vRNG = { if(is.numeric(mm) || is.complex(mm))
  +                                    function(n) rpois(n,lambda= 0.75)# <- about 47% zeros
  +                       else ## logical
  +                           function(n) runif(n) > 0.8 }, ## 80% zeros
  +                       showOnly=FALSE)
  + {
  +     stopifnot(is(ms,"sparseMatrix"))
  +     d <- dim(ms)
  +     s1 <- function(n) sample(n, pmin(n, pmax(1, rpois(1, n.uniq))))
  +     i <- sample(s1(d[1]), k/2+ rpois(1, k/2), replace = TRUE)
  +     j <- sample(s1(d[2]), k/2+ rpois(1, k/2), replace = TRUE)
  +     assert.EQ.mat(ms[i,j], mm[i,j])
  +     ms2 <- ms. <- ms; mm. <- mm # save
  +     ## now sub*assign* to these repeated indices, and then compare -----
  +     v <- vRNG(length(i) * length(j))
  +     mm[i,j] <- v
  +     ms[i,j] <- v
  +     ## useful to see (ii,ij), but confusing R/ESS when additionally debugging:
  +     ## if(!showOnly && interactive()) { op <- options(error = recover); on.exit(options(op)) }
  +     assert.EQ.mat(ms, mm, show=showOnly)
  +     ## vector indexing m[cbind(i,j)] == m[i + N(j-1)] ,  N = nrow(.)
  +     ii <- seq_len(min(length(i), length(j)))
  +     i <- i[ii]
  +     j <- j[ii]
  +     ij <- cbind(i, j)
  +     ii <- i + nrow(ms)*(j - 1)
  +     ord.i <- order(ii)
  +     iio <- ii[ord.i]
  +     ui <- unique(iio) # compare these with :
  +     neg.ii <- - setdiff(seq_len(prod(d)), ii)
  +     stopifnot(identical(mm[ii], mm[ij]),
  +               identical(ms.[ui], ms.[neg.ii]),
  + 	      ms.[ij] == mm.[ii], ## M[ cbind(i,j) ] was partly broken; now checking
  +               ms.[ii] == mm.[ii])
  +     v <- v[seq_len(length(i))]
  +     if(is(ms,"nMatrix")) v <- as.logical(v)  # !
  +     mm.[ij] <- v
  +     ms.[ii] <- v
  +     nodup <- (length(ui) == length(ii)) ## <==>  ! anyDuplicated(iio)
  +     if(nodup) { cat("[-]") # rare, unfortunately
  + 	ms2[neg.ii] <- v[ord.i]
  + 	stopifnot(identical(ms2, ms.))
  +     }
  +     assert.EQ.mat(ms., mm., show=showOnly)
  + } ##{chkAssign}
  > 
  > ## Get duplicated index {because these are "hard" (and rare)
  > getDuplIndex <- function(n, k) {
  +     repeat {
  +         i <- sample(n, k, replace=TRUE) # 3 4 6 9 2 9 :  9 is twice
  +         if(anyDuplicated(i)) break
  +     }
  +     i
  + }
  > 
  > ## From package 'sfsmisc':
  > repChar <- function (char, no) paste(rep.int(char, no), collapse = "")
  > 
  > m <- 1:800
  > set.seed(101) ; m[sample(800, 600)] <- 0
  > m0 <- Matrix(m, nrow = 40)
  > m1 <- add.simpleDimnames(m0)
  > for(kind in c("n", "l", "d")) {
  +  for(m in list(m0,m1)) { ## -- with and without dimnames -------------------------
  +     kClass <-paste0(kind, "Matrix"  )
  +     Ckind <- paste0(kind, "gCMatrix")
  +     Tkind <- paste0(kind, "gTMatrix")
  +     str(mC <- as(m, Ckind))
  +     str(mT <- as(as(as(m, kClass), "TsparseMatrix"), Tkind))
  +     mm <- as(mC, "matrix") # also logical or double
  +     IDENT <- if(kind == "n") function(x,y) Q.eq2(x,y, tol=0) else identical
  +     stopifnot(identical(mT, as(as(mC, "TsparseMatrix"), Tkind)),
  +               identical(mC, as(mT, Ckind)),
  +               Qidentical(mC[0,0], new(Ckind)),
  +               Qidentical(mT[0,0], new(Tkind)),
  +               identical(unname(mT[0,]), new(Tkind, Dim = c(0L,ncol(m)))),
  +               identical(unname(mT[,0]), new(Tkind, Dim = c(nrow(m),0L))),
  +               IDENT(mC[0,], as(mT[FALSE,], Ckind)),
  +               IDENT(mC[,0], as(mT[,FALSE], Ckind)),
  +               sapply(pmin(min(dim(mC)), c(0:2, 5:10)),
  +                      function(k) {i <- seq_len(k); all(mC[i,i] == mT[i,i])}),
  +               TRUE)
  +     cat("ok\n")
  +     show(mC[,1])
  +     show(mC[1:2,])
  +     show(mC[7,  drop = FALSE])
  +     assert.EQ.mat(mC[1:2,],   mm[1:2,])
  +     assert.EQ.mat(mC[0,],     mm[0,])
  +     assert.EQ.mat(mC[,FALSE], mm[,FALSE])
  +     ##
  +     ## *repeated* (aka 'duplicated') indices - did not work at all ...
  +     i <- pmin(nrow(mC), rep(8:10,2))
  +     j <- c(2:4, 4:3)
  +     assert.EQ.mat(mC[i,], mm[i,])
  +     assert.EQ.mat(mC[,j], mm[,j])
  +     ## FIXME? assert.EQ.mat(mC[,NA], mm[,NA]) -- mC[,NA] is all 0 "instead" of all NA
  +     ## MM currently thinks we should  NOT  allow  <sparse>[ <NA> ]
  +     assert.EQ.mat(mC[i, 2:1], mm[i, 2:1])
  +     assert.EQ.mat(mC[c(4,1,2:1), j], mm[c(4,1,2:1), j])
  +     assert.EQ.mat(mC[i,j], mm[i,j])
  +     ##
  +     ## set.seed(7)
  +     op <- options(Matrix.verbose = FALSE)
  +     cat(" for(): ")
  +     for(n in 1:(if(doExtras) 50 else 5)) { # (as chkAssign() is random)
  + 	chkAssign(mC, mm)
  + 	chkAssign(mC[-3,-2], mm[-3,-2])
  +         cat(".")
  +     }
  +     options(op)
  +     cat(sprintf("\n[Ok]%s\n\n", repChar("-", 64)))
  +  }
  +  cat(sprintf("\nok( %s )\n== ###%s\n\n", kind, repChar("=", 70)))
  + }## end{for}---------------------------------------------------------------
  Formal class 'ngCMatrix' [package "Matrix"] with 5 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ p       : int [1:21] 0 7 20 28 38 49 59 65 73 84 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ factors : list()
  Formal class 'ngTMatrix' [package "Matrix"] with 5 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ j       : int [1:200] 0 0 0 0 0 0 0 1 1 1 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ factors : list()
  Note: method with signature 'nsparseMatrix#sparseMatrix' chosen for function '==',
   target signature 'ngCMatrix#ngTMatrix'.
   "nMatrix#nMatrix", "sparseMatrix#nsparseMatrix" would also be valid
  ok
   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
  [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [25] FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE
  [37] FALSE FALSE FALSE FALSE
  2 x 20 sparse Matrix of class "ngCMatrix"
                                              
  [1,] . | | | . . . . | . . | . . . | . . . |
  [2,] . | . . | . . . | . . . . . | | . . . |
  [1] FALSE
   for(): .....
  [Ok]----------------------------------------------------------------
  
  Formal class 'ngCMatrix' [package "Matrix"] with 5 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ p       : int [1:21] 0 7 20 28 38 49 59 65 73 84 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
    .. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
    ..@ factors : list()
  Formal class 'ngTMatrix' [package "Matrix"] with 5 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ j       : int [1:200] 0 0 0 0 0 0 0 1 1 1 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
    .. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
    ..@ factors : list()
  ok
     r1    r2    r3    r4    r5    r6    r7    r8    r9   r10   r11   r12   r13 
  FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE 
    r14   r15   r16   r17   r18   r19   r20   r21   r22   r23   r24   r25   r26 
  FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
    r27   r28   r29   r30   r31   r32   r33   r34   r35   r36   r37   r38   r39 
   TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE 
    r40 
  FALSE 
  2 x 20 sparse Matrix of class "ngCMatrix"
     [[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
                                            
  r1 . | | | . . . . | . . | . . . | . . . |
  r2 . | . . | . . . | . . . . . | | . . . |
  [1] FALSE
   for(): .....
  [Ok]----------------------------------------------------------------
  
  
  ok( n )
  == ###======================================================================
  
  Formal class 'lgCMatrix' [package "Matrix"] with 6 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ p       : int [1:21] 0 7 20 28 38 49 59 65 73 84 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ x       : logi [1:200] TRUE TRUE TRUE TRUE TRUE TRUE ...
    ..@ factors : list()
  Formal class 'lgTMatrix' [package "Matrix"] with 6 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ j       : int [1:200] 0 0 0 0 0 0 0 1 1 1 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ x       : logi [1:200] TRUE TRUE TRUE TRUE TRUE TRUE ...
    ..@ factors : list()
  ok
   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
  [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [25] FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE
  [37] FALSE FALSE FALSE FALSE
  2 x 20 sparse Matrix of class "lgCMatrix"
                                              
  [1,] . | | | . . . . | . . | . . . | . . . |
  [2,] . | . . | . . . | . . . . . | | . . . |
  [1] FALSE
   for(): .....
  [Ok]----------------------------------------------------------------
  
  Formal class 'lgCMatrix' [package "Matrix"] with 6 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ p       : int [1:21] 0 7 20 28 38 49 59 65 73 84 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
    .. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
    ..@ x       : logi [1:200] TRUE TRUE TRUE TRUE TRUE TRUE ...
    ..@ factors : list()
  Formal class 'lgTMatrix' [package "Matrix"] with 6 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ j       : int [1:200] 0 0 0 0 0 0 0 1 1 1 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
    .. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
    ..@ x       : logi [1:200] TRUE TRUE TRUE TRUE TRUE TRUE ...
    ..@ factors : list()
  ok
     r1    r2    r3    r4    r5    r6    r7    r8    r9   r10   r11   r12   r13 
  FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE 
    r14   r15   r16   r17   r18   r19   r20   r21   r22   r23   r24   r25   r26 
  FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
    r27   r28   r29   r30   r31   r32   r33   r34   r35   r36   r37   r38   r39 
   TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE 
    r40 
  FALSE 
  2 x 20 sparse Matrix of class "lgCMatrix"
     [[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
                                            
  r1 . | | | . . . . | . . | . . . | . . . |
  r2 . | . . | . . . | . . . . . | | . . . |
  [1] FALSE
   for(): .....
  [Ok]----------------------------------------------------------------
  
  
  ok( l )
  == ###======================================================================
  
  Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ p       : int [1:21] 0 7 20 28 38 49 59 65 73 84 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ x       : num [1:200] 8 27 29 30 32 33 35 41 42 44 ...
    ..@ factors : list()
  Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ j       : int [1:200] 0 0 0 0 0 0 0 1 1 1 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ x       : num [1:200] 8 27 29 30 32 33 35 41 42 44 ...
    ..@ factors : list()
  ok
   [1]  0  0  0  0  0  0  0  8  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
  [26]  0 27  0 29 30  0 32 33  0 35  0  0  0  0  0
  2 x 20 sparse Matrix of class "dgCMatrix"
                                                              
  [1,] . 41 81 121   . . . . 321 . . 441 . .   . 601 . . . 761
  [2,] . 42  .   . 162 . . . 322 . .   . . . 562 602 . . . 762
  [1] 0
   for(): .....
  [Ok]----------------------------------------------------------------
  
  Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ p       : int [1:21] 0 7 20 28 38 49 59 65 73 84 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
    .. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
    ..@ x       : num [1:200] 8 27 29 30 32 33 35 41 42 44 ...
    ..@ factors : list()
  Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
    ..@ i       : int [1:200] 7 26 28 29 31 32 34 0 1 3 ...
    ..@ j       : int [1:200] 0 0 0 0 0 0 0 1 1 1 ...
    ..@ Dim     : int [1:2] 40 20
    ..@ Dimnames:List of 2
    .. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
    .. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
    ..@ x       : num [1:200] 8 27 29 30 32 33 35 41 42 44 ...
    ..@ factors : list()
  ok
   r1  r2  r3  r4  r5  r6  r7  r8  r9 r10 r11 r12 r13 r14 r15 r16 r17 r18 r19 r20 
    0   0   0   0   0   0   0   8   0   0   0   0   0   0   0   0   0   0   0   0 
  r21 r22 r23 r24 r25 r26 r27 r28 r29 r30 r31 r32 r33 r34 r35 r36 r37 r38 r39 r40 
    0   0   0   0   0   0  27   0  29  30   0  32  33   0  35   0   0   0   0   0 
  2 x 20 sparse Matrix of class "dgCMatrix"
     [[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
                                                            
  r1 . 41 81 121   . . . . 321 . . 441 . .   . 601 . . . 761
  r2 . 42  .   . 162 . . . 322 . .   . . . 562 602 . . . 762
  [1] 0
   for(): .....
  [Ok]----------------------------------------------------------------
  
  
  ok( d )
  == ###======================================================================
  
  > showProc.time()
  Time elapsed:  2 0 2 
  > 
  > if(doExtras) {### {was ./AAA_index.R, MM-only}
  +     ## an nsparse-example
  +     A <- Matrix(c(rep(c(1,0,0),2), rep(c(2,0),7), c(0,0,2), rep(0,4)), 3,9)
  +     i <- c(3,1:2)
  +     j <- c(3, 5, 9, 5, 9)
  +     vv <- logical(length(i)*length(j)); vv[6:9] <- TRUE
  + 
  +     print(An <- as(A,"nMatrix")); an <- as(An, "matrix")
  +     assert.EQ.mat(An, an)
  +     An[i, j] <- vv
  +     an[i, j] <- vv
  +     assert.EQ.mat(An, an)## error
  +     if(!all(An == an)) show(drop0(An - an))
  +     ## all are +1
  + 
  +     options("Matrix.subassign.verbose" = TRUE)# output from C
  +     An <- as(A,"nMatrix"); An[i, j] <- vv
  +     ## and compare with this:
  +     Al <- as(A,"lMatrix"); Al[i, j] <- vv
  +     options("Matrix.subassign.verbose" = FALSE)
  + 
  +     ##--- An interesting not small not large example  for  M[i,j] <- v ------------
  +     ##
  +     M <- Matrix(c(1, rep(0,7), 1:4), 3,4)
  +     N0 <- kronecker(M,M)
  +     mkN1 <- function(M) {
  +         stopifnot(length(d <- dim(M)) == 2)
  +         isC <- is(M,"CsparseMatrix")
  +         M[,d[2]] <- c(0,2,0)
  +         N <- kronecker(diag(x=1:2), M)## remains sparse if 'M' is
  +         if(isC) N <- as(N, "CsparseMatrix")
  +         diag(N[-1,]) <- -2
  +         N[9,]  <- 1:4   # is recycled
  +         N[,12] <- -7:-9 # ditto
  +         N
  +     }
  + 
  +     show(N1 <- t(N <- mkN1(N0)))    # transpose {for display reasons}
  +     C1 <- t(C <- mkN1(as(N0,"CsparseMatrix")))
  +     stopifnot(all(C == N))
  +     assert.EQ.mat(C, mkN1(.asmatrix(N0)))
  + 
  +     C. <- C1
  +     show(N <- N1) ; n <- .asmatrix(N); str(N)
  +     sort(i <- c(6,8,19,11,21,20,10,7,12,9,5,18,17,22,13))## == c(5:13, 17:22))
  +     sort(j <- c(3,8,6,15,10,4,14,13,16,2,11,17,7,5))## == c(2:8, 10:11, 13:17)
  +     val <- v.l <- 5*c(0,6,0,7,0,0,8:9, 0,0)
  +     show(spv <- as(val, "sparseVector")); str(spv)
  + 
  +     n [i,j] <- v.l
  +     N [i,j] <- val# is recycled, too
  +     C.[i,j] <- val
  +     assert.EQ.mat(N,n) ; stopifnot(all(C. == N))
  +     ## and the same *again*:
  +     n [i,j] <- v.l
  +     N [i,j] <- val
  +     C.[i,j] <- val
  +     assert.EQ.mat(N,n)
  +     stopifnot(all(C. == N))
  + 
  +     print(load(system.file("external", "symA.rda", package="Matrix"))) # "As"
  +     stopifnotValid(As, "dsCMatrix"); stopifnot(identical(As@factors, list()))
  +     R. <- drop0(chol(As))
  +     stopifnot(1:32 == sort(diag(R.)), ## !
  +               R.@x == as.integer(R.@x),## so it is an integer-valued chol-decomp !
  +               ## shows that (1) As is *not* singular  (2) the matrix is not random
  +               all.equal(crossprod(R.), As, tolerance =1e-15))
  +     print(summary(evA <- eigen(As, only.values=TRUE)$val))
  +     print(tail(evA)) ## largest three ~= 10^7,  smallest two *negative*
  +     print(rcond(As)) # 1.722 e-21 == very bad !
  +     ##-> this *is* a border line case, i.e. very close to singular !
  +     ## and also determinant(.) is rather random here!
  +     cc0 <- Cholesky(As)# no problem
  +     try({
  +         cc <- Cholesky(As, super=TRUE)
  +         ## gives --on 32-bit only--
  +         ## Cholmod error 'matrix not positive definite' at file:../Supernodal/t_cholmod_super_numeric.c, line 613
  +         ecc <- expand(cc)
  +         L.P <- with(ecc, crossprod(L,P))  ## == L'P
  +         ## crossprod(L.P) == (L'P)' L'P == P'LL'P
  +         stopifnot( all.equal(crossprod(L.P), As) )
  +     })
  +     ##---- end{ eigen( As ) -----------
  + 
  + } ## only if(doExtras)
  > 
  > 
  > ##---- Symmetric indexing of symmetric Matrix ----------
  > m. <- mC
  > m.[, c(2, 7:12)] <- 0
  > stopifnotValid(S <- crossprod(add.simpleDimnames(m.) %% 100), "dsCMatrix")
  > ss <- as(S, "matrix")
  > ds <- as(S, "denseMatrix")
  > ## NA-indexing of *dense* Matrices: should work as traditionally
  > assert.EQ.mat(ds[NA,NA], ss[NA,NA])
  > assert.EQ.mat(ds[NA,  ], ss[NA,])
  > assert.EQ.mat(ds[  ,NA], ss[,NA])
  > T <- as(S, "TsparseMatrix")
  > stopifnot(identical(ds[2 ,NA], ss[2,NA]),
  + 	  identical(ds[NA, 1], ss[NA, 1]),
  + 	  identical(S, as(T, "CsparseMatrix")) )
  > 
  > ## non-repeated indices:
  > i <- c(7:5, 2:4);assert.EQ.mat(T[i,i], ss[i,i])
  > ## NA in indices  -- check that we get a helpful error message:
  > i[2] <- NA
  > er <- tryCatch(T[i,i], error = function(e)e)
  > stopifnot(as.logical(grep("indices.*sparse Matrices", er$message)))
  > 
  > N <- nrow(T)
  > set.seed(11)
  > for(n in 1:(if(doExtras) 50 else 3)) {
  +     i <- sample(N, max(2, sample(N,1)), replace = FALSE)
  +     validObject(Tii <- T[i,i]) ; tTi <- t(T)[i,i]
  +     stopifnot(is(Tii, "dsTMatrix"), # remained symmetric Tsparse
  + 	      is(tTi, "dsTMatrix"), # may not be identical when *sorted* differently
  + 	      identical(as(t(Tii),"CsparseMatrix"), as(tTi,"CsparseMatrix")))
  +     assert.EQ.mat(Tii, ss[i,i])
  + }
  > 
  > b <- diag(1:2)[,c(1,1,2,2)]
  > cb <- crossprod(b)
  > cB <- crossprod(Matrix(b, sparse=TRUE))
  > a <- matrix(0, 6, 6)
  > a[1:4, 1:4] <- cb
  > A1 <- A2 <- Matrix(0, 6, 6)#-> sparse
  > A1[1:4, 1:4] <- cb
  > A2[1:4, 1:4] <- cB
  > assert.EQ.mat(A1, a)# indeed
  > stopifnot(identical(A1, A2), is(A1, "dsCMatrix"))
  > 
  > ## repeated ones ``the challenge'' (to do smartly):
  > j <- c(4, 4, 9, 12, 9, 4, 17, 3, 18, 4, 12, 18, 4, 9)
  > assert.EQ.mat(T[j,j], ss[j,j])
  > ## and another two sets  (a, A) &  (a., A.) :
  > a <- matrix(0, 6,6)
  > a[upper.tri(a)] <- (utr <- c(2, 0,-1, 0,0,5, 7,0,0,0, 0,0,-2,0,8))
  > ta <- t(a); ta[upper.tri(a)] <- utr; a <- t(ta)
  > diag(a) <- c(0,3,0,4,6,0)
  > A <- as(Matrix(a), "TsparseMatrix")
  > A. <- A
  > diag(A.) <- 10 * (1:6)
  > a. <- as(A., "matrix")
  > ## More testing {this was not working for a long time..}
  > set.seed(1)
  > for(n in 1:(if(doExtras) 100 else 6)) {
  +     i <- sample(1:nrow(A), 3+2*rpois(1, lam=3), replace=TRUE)
  +     Aii  <- A[i,i]
  +     A.ii <- A.[i,i]
  +     stopifnot(class(Aii) == class(A),
  +               class(A.ii) == class(A.))
  +     assert.EQ.mat(Aii , a [i,i])
  +     assert.EQ.mat(A.ii, a.[i,i])
  +     assert.EQ.mat(T[i,i], ss[i,i])
  + }
  > showProc.time()
  Time elapsed:  0.56 0 0.56 
  > 
  > stopifnot(all.equal(mC[,3], mm[,3]),
  + 	  identical(mC[ij], mC[ij + 0.4]),
  + 	  identical(mC[ij], mm[ij]),
  + 	  identical(mC[iN], mm[iN]))
  > ## out of bound indexing must be detected:
  > assertError(mC[cbind(ij[,1] - 5, ij[,2])])
  > assertError(mC[cbind(ij[,1],     ij[,2] + ncol(mC))])
  > 
  > assert.EQ.mat(mC[7, , drop=FALSE], mm[7, , drop=FALSE])
  > identical    (mC[7,   drop=FALSE], mm[7,   drop=FALSE]) # *vector* indexing
  [1] TRUE
  > 
  > stopifnot(dim(mC[numeric(0), ]) == c(0,20), # used to give warnings
  +           dim(mC[, integer(0)]) == c(40,0),
  +           identical(mC[, integer(0)], mC[, FALSE]))
  > validObject(print(mT[,c(2,4)]))
  40 x 2 sparse Matrix of class "dgTMatrix"
      c2  c4
  r1  41 121
  r2  42   .
  r3   .   .
  r4  44   .
  r5  45   .
  r6   .   .
  r7   .   .
  r8  48 128
  r9   .   .
  r10  .   .
  r11  .   .
  r12  . 132
  r13  .   .
  r14  . 134
  r15  . 135
  r16  . 136
  r17  .   .
  r18  .   .
  r19 59   .
  r20  .   .
  r21 61   .
  r22  .   .
  r23 63   .
  r24  .   .
  r25  .   .
  r26 66 146
  r27  .   .
  r28 68   .
  r29 69   .
  r30  .   .
  r31 71   .
  r32  .   .
  r33  . 153
  r34  . 154
  r35  .   .
  r36  .   .
  r37 77   .
  r38  .   .
  r39  .   .
  r40  . 160
  [1] TRUE
  > stopifnot(all.equal(mT[2,], mm[2,]),
  + 	  ## row or column indexing in combination with t() :
  + 	  Q.C.identical(mT[2,], t(mT)[,2]),
  + 	  Q.C.identical(mT[-2,], t(t(mT)[,-2])),
  + 	  Q.C.identical(mT[c(2,5),], t(t(mT)[,c(2,5)])) )
  > assert.EQ.mat(mT[4,, drop = FALSE], mm[4,, drop = FALSE])
  > stopifnot(identical3(mm[,1], mC[,1], mT[,1]),
  + 	  identical3(mm[3,], mC[3,], mT[3,]),
  + 	  identical3(mT[2,3], mC[2,3], 0),
  + 	  identical(mT[], mT),
  +           identical4(       mm[c(3,7), 2:4],  as.mat( m[c(3,7), 2:4]),
  +                      as.mat(mT[c(3,7), 2:4]), as.mat(mC[c(3,7), 2:4]))
  +           )
  > 
  > x.x <- crossprod(mC)
  > stopifnot(class(x.x) == "dsCMatrix",
  +           class(x.x. <- round(x.x / 10000)) == "dsCMatrix",
  +           identical(x.x[cbind(2:6, 2:6)],
  +                     diag(x.x [2:6, 2:6])))
  > head(x.x.) # Note the *non*-structural 0's printed as "0"
  6 x 20 sparse Matrix of class "dgCMatrix"
     [[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
                                                             
  c1 1 0 .  1  1  1  .  .  1  5  .  6  5  .  .  2  4  3  7  6
  c2 0 5 2  2  1  3  4  .  5  4  6  5  9 11  7 13  5 20 11 23
  c3 . 2 8  4  6  5  3  5  7  .  9 19  5  7 12  5  . 13 22 20
  c4 1 2 4 20  2  9 12 17 19  5 31 11 28  . 25 15 10 49 52 52
  c5 1 1 6  2 35 13  9  5 24 13  8 36  . 11 30 42 23 59 26 51
  c6 1 3 5  9 13 51 12 21 30 17 19 44 45 24 40 53 16 15 49 35
  > tail(x.x., -3) # all but the first three lines
  17 x 20 sparse Matrix of class "dgCMatrix"
     [[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
                                                                             
  c4  1  2  4 20  2  9 12  17  19   5  31  11  28   .  25  15  10  49  52  52
  c5  1  1  6  2 35 13  9   5  24  13   8  36   .  11  30  42  23  59  26  51
  c6  1  3  5  9 13 51 12  21  30  17  19  44  45  24  40  53  16  15  49  35
  c7  .  4  3 12  9 12 42   9  27   .  35   .  28   .  16  16  33  54  41  40
  c8  .  .  5 17  5 21  9  71  11  13  38  55  45   .  53  17  41  61 110  45
  c9  1  5  7 19 24 30 27  11 127  25  30  47  34  39  59  79  68  47  50 157
  c10 5  4  .  5 13 17  .  13  25 160  15 108  58  41  43  71  77 106 141 117
  c11 .  6  9 31  8 19 35  38  30  15 144  19  64  22  50   .  55 148 126  32
  c12 6  5 19 11 36 44  .  55  47 108  19 302  71  51  81  56 124  96 173  68
  c13 5  9  5 28  . 45 28  45  34  58  64  71 253  27  59  93   .  70 183 117
  c14 . 11  7  . 11 24  .   .  39  41  22  51  27 179  33 102  35  75   .  86
  c15 .  7 12 25 30 40 16  53  59  43  50  81  59  33 302  71   .   . 125  86
  c16 2 13  5 15 42 53 16  17  79  71   .  56  93 102  71 380   . 128  44 238
  c17 4  5  . 10 23 16 33  41  68  77  55 124   .  35   .   . 351 230 102 105
  c18 3 20 13 49 59 15 54  61  47 106 148  96  70  75   . 128 230 682 103 268
  c19 7 11 22 52 26 49 41 110  50 141 126 173 183   . 125  44 102 103 660 228
  c20 6 23 20 52 51 35 40  45 157 117  32  68 117  86  86 238 105 268 228 841
  > 
  > lx.x <- as(x.x, "lsCMatrix") # FALSE only for "structural" 0
  > (l10 <- lx.x[1:10, 1:10])# "lsC"
  10 x 10 sparse Matrix of class "lsCMatrix"
     [[ suppressing 10 column names 'c1', 'c2', 'c3' ... ]]
                         
  c1  | | . | | | . . | |
  c2  | | | | | | | . | |
  c3  . | | | | | | | | .
  c4  | | | | | | | | | |
  c5  | | | | | | | | | |
  c6  | | | | | | | | | |
  c7  . | | | | | | | | .
  c8  . . | | | | | | | |
  c9  | | | | | | | | | |
  c10 | | . | | | . | | |
  > (l3 <-  lx.x[1:3, ])
  3 x 20 sparse Matrix of class "lgCMatrix"
     [[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
                                            
  c1 | | . | | | . . | | . | | . . | | | | |
  c2 | | | | | | | . | | | | | | | | | | | |
  c3 . | | | | | | | | . | | | | | | . | | |
  > m.x <- as.mat(x.x) # as.mat() *drops* (NULL,NULL) dimnames
  > stopifnot(class(l10) == "lsCMatrix", # symmetric indexing -> symmetric !
  +           identical(as.mat(lx.x), m.x != 0),
  +           identical(as.logical(lx.x), as.logical(m.x)),
  +           identical(as.mat(l10), m.x[1:10, 1:10] != 0),
  +           identical(as.mat(l3 ), m.x[1:3, ] != 0)
  +           )
  > 
  > ##-- Sub*assignment* with repeated / duplicated index:
  > A <- Matrix(0,4,3) ; A[c(1,2,1), 2] <- 1 ; A
  4 x 3 sparse Matrix of class "dgCMatrix"
            
  [1,] . 1 .
  [2,] . 1 .
  [3,] . . .
  [4,] . . .
  > B <- A;              B[c(1,2,1), 2] <- 1:3; B; B. <- B
  4 x 3 sparse Matrix of class "dgCMatrix"
            
  [1,] . 3 .
  [2,] . 2 .
  [3,] . . .
  [4,] . . .
  > B.[3,] <- rbind(4:2)
  > ## change the diagonal and the upper and lower subdiagonal :
  > diag(B.) <- 10 * diag(B.)
  > diag(B.[,-1]) <- 5* diag(B.[,-1])
  > diag(B.[-1,]) <- 4* diag(B.[-1,]) ; B.
  4 x 3 sparse Matrix of class "dgCMatrix"
              
  [1,] . 15  .
  [2,] . 20  .
  [3,] 4 12 20
  [4,] .  .  .
  > C <- B.; C[,2] <- C[,2];  C[1,] <- C[1,]; C[2:3,2:1] <- C[2:3,2:1]
  > stopifnot(identical(unname(.asmatrix(A)),
  + 		    local({a <- matrix(0,4,3); a[c(1,2,1), 2] <-  1 ; a})),
  + 	  identical(unname(.asmatrix(B)),
  + 		    local({a <- matrix(0,4,3); a[c(1,2,1), 2] <- 1:3; a})),
  + 	  identical(C, drop0(B.)))
  > ## <sparse>[<logicalSparse>] <- v  failed in the past
  > T <- as(C,"TsparseMatrix"); C. <- C
  > T[T>0] <- 21
  Note: method with signature 'TsparseMatrix#Matrix#missing#replValue' chosen for function '[<-',
   target signature 'dgTMatrix#lgTMatrix#missing#numeric'.
   "Matrix#lsparseMatrix#missing#replValue" would also be valid
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  > C[C>0] <- 21
  Note: method with signature 'CsparseMatrix#Matrix#missing#replValue' chosen for function '[<-',
   target signature 'dgCMatrix#lgCMatrix#missing#numeric'.
   "Matrix#lsparseMatrix#missing#replValue" would also be valid
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  > a. <- local({a <- .asmatrix(C.); a[a>0] <- 21; a})
  > assert.EQ.mat(C, a.)
  > stopifnot(identical(C, as(T, "CsparseMatrix")))
  > 
  > ## used to fail
  > n <- 5 ## or much larger
  > sm <- new("dsTMatrix", i=1L, j=1L, Dim=as.integer(c(n,n)), x = 1)
  > (cm <- as(sm, "CsparseMatrix"))
  5 x 5 sparse Matrix of class "dsCMatrix"
                
  [1,] . . . . .
  [2,] . 1 . . .
  [3,] . . . . .
  [4,] . . . . .
  [5,] . . . . .
  > sm[2,]
  [1] 0 1 0 0 0
  > stopifnot(sm[2,] == c(0:1, rep.int(0,ncol(sm)-2)),
  + 	  sm[2,] == cm[2,],
  + 	  sm[,3] == sm[3,],
  + 	  all(sm[,-(1:3)] == t(sm[-(1:3),])), # all(<lge.>)
  + 	  all(sm[,-(1:3)] == 0)
  + 	  )
  > showProc.time()
  Time elapsed:  0.5 0 0.53 
  > 
  > ##---  "nsparse*" sub-assignment :----------
  > M <- Matrix(c(1, rep(0,7), 1:4), 3,4)
  > N0 <- kronecker(M,M)
  > Nn <- as(N0, "nMatrix"); nn <- as(Nn,"matrix")
  > (Nn00 <- Nn0 <- Nn); nn00 <- nn0 <- nn
  9 x 16 sparse Matrix of class "ngTMatrix"
                                       
   [1,] | . . | . . . . . . . . | . . |
   [2,] . . . | . . . . . . . . . . . |
   [3,] . . | | . . . . . . . . . . | |
   [4,] . . . . . . . . . . . . | . . |
   [5,] . . . . . . . . . . . . . . . |
   [6,] . . . . . . . . . . . . . . | |
   [7,] . . . . . . . . | . . | | . . |
   [8,] . . . . . . . . . . . | . . . |
   [9,] . . . . . . . . . . | | . . | |
  > 
  > set.seed(1)
  > Nn0 <- Nn00; nn0 <- nn00
  > for(i in 1:(if(doExtras) 200 else 25)) {
  +     Nn <- Nn0
  +     nn <- nn0
  +     i. <- getDuplIndex(nrow(N0), 6)
  +     j. <- getDuplIndex(ncol(N0), 4)
  +     vv <- sample(c(FALSE,TRUE),
  +                  length(i.)*length(j.), replace=TRUE)
  +     cat(",")
  +     Nn[i., j.] <- vv
  +     nn[i., j.] <- vv
  +     assert.EQ.mat(Nn, nn)
  +     if(!all(Nn == nn)) {
  +         cat("i=",i,":\n i. <- "); dput(i.)
  +         cat("j. <- "); dput(j.)
  +         cat("which(vv): "); dput(which(vv))
  +         cat("Difference matrix:\n")
  +         show(drop0(Nn - nn))
  +     }
  +     cat("k")
  +     ## sub-assign double precision to logical sparseMatrices: now *with* warning:
  +     ##  {earlier: gave *no* warning}:
  +     assertWarning(Nn[1:2,] <- -pi)
  +     assertWarning(Nn[, 5] <- -pi)
  +     assertWarning(Nn[2:4, 5:8] <- -pi)
  +     stopifnotValid(Nn,"nsparseMatrix")
  +     ##
  +     cat(".")
  +     if(i %% 10 == 0) cat("\n")
  +     if(i == 100) {
  +         Nn0 <- as(Nn0, "CsparseMatrix")
  +         cat("Now: class", class(Nn0)," :\n~~~~~~~~~~~~~~~~~\n")
  +     }
  + }
  ,k.,k.,k.,k.,k.,k.,k.,k.,k.,k.
  ,k.,k.,k.,k.,k.,k.,k.,k.,k.,k.
  ,k.,k.,k.,k.,k.> showProc.time()
  Time elapsed:  0.47 0 0.47 
  > Nn <- Nn0
  > ## Check that  NA is interpreted as TRUE (with a warning), for "nsparseMatrix":
  > assertWarning(Nn[ii <-     3  ] <- NA); stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii])
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  > assertWarning(Nn[ii <-   22:24] <- NA); stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii])
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  > assertWarning(Nn[ii <- -(1:99)] <- NA); stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii])
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  > assertWarning(Nn[ii <- 3:4  ] <- c(0,NA))
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  > stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii] == 0:1)
  > assertWarning(Nn[ii <- 25:27] <- c(0,1,NA))
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  > stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii] == c(FALSE,TRUE,TRUE))
  > 
  > m0 <- Diagonal(5)
  > stopifnot(identical(m0[2,], m0[,2]),
  + 	  identical(m0[,1], c(1,0,0,0,0)))
  > ### Diagonal -- Sparse:
  > (m1 <- as(m0, "TsparseMatrix")) # dtTMatrix unitriangular
  5 x 5 sparse Matrix of class "dtTMatrix" (unitriangular)
                
  [1,] I . . . .
  [2,] . I . . .
  [3,] . . I . .
  [4,] . . . I .
  [5,] . . . . I
  > (m2 <- as(m0, "CsparseMatrix")) # dtCMatrix unitriangular
  5 x 5 sparse Matrix of class "dtCMatrix" (unitriangular)
                
  [1,] I . . . .
  [2,] . I . . .
  [3,] . . I . .
  [4,] . . . I .
  [5,] . . . . I
  > m1g <- as(m1, "generalMatrix")
  > tr1 <- as(m1, "denseMatrix") # dtrMatrix unitriangular
  > stopifnotValid(m1g, "dgTMatrix")
  > diag(tr1) <- 100
  > stopifnot(diag(tr1) == 100)# failed when 'diag<-' did not recycle
  > assert.EQ.mat(m2[1:3,],    diag(5)[1:3,])
  > assert.EQ.mat(m2[,c(4,1)], diag(5)[,c(4,1)])
  > stopifnot(identical(m2[1:3,], as(m1[1:3,], "CsparseMatrix")),
  +           identical(uniqTsparse(m1[, c(4,2)]),
  +                     uniqTsparse(as(m2[, c(4,2)], "TsparseMatrix")))
  +           )## failed in 0.9975-11
  > 
  > (uTr <- new("dtTMatrix", Dim = c(3L,3L), diag="U"))
  3 x 3 sparse Matrix of class "dtTMatrix" (unitriangular)
            
  [1,] I . .
  [2,] . I .
  [3,] . . I
  > uTr[1,] <- 0
  > assert.EQ.mat(uTr, cbind(0, rbind(0,diag(2))))
  > 
  > M <- m0; M[1,] <- 0
  > Z <- m0; Z[] <- 0; z <- array(0, dim(M))
  > stopifnot(identical(M, Diagonal(x=c(0, rep(1,4)))),
  +           all(Z == 0), Qidentical(as(Z, "matrix"), z))
  > M <- m0; M[,3] <- 3 ; M ; stopifnot(is(M, "sparseMatrix"), M[,3] == 3)
  5 x 5 sparse Matrix of class "dgCMatrix"
                
  [1,] 1 . 3 . .
  [2,] . 1 3 . .
  [3,] . . 3 . .
  [4,] . . 3 1 .
  [5,] . . 3 . 1
  > checkMatrix(M)
  Note: method with signature 'sparseMatrix#ldiMatrix' chosen for function '==',
   target signature 'nsCMatrix#ldiMatrix'.
   "nsparseMatrix#sparseMatrix", "nMatrix#lMatrix" would also be valid
   norm(m [5 x 5]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > M <- m0; M[1:3, 3] <- 0 ;M
  5 x 5 diagonal matrix of class "ddiMatrix"
       [,1] [,2] [,3] [,4] [,5]
  [1,]    1    .    .    .    .
  [2,]    .    1    .    .    .
  [3,]    .    .    0    .    .
  [4,]    .    .    .    1    .
  [5,]    .    .    .    .    1
  > T <- m0; T[1:3, 3] <- 10
  > stopifnot(identical(M, Diagonal(x=c(1,1, 0, 1,1))),
  +           isValid(T, "triangularMatrix"), identical(T[,3], c(10,10,10,0,0)))
  > 
  > M <- m1; M[1,] <- 0 ; M ; assert.EQ.mat(M, diag(c(0,rep(1,4))), tol=0)
  5 x 5 sparse Matrix of class "dtTMatrix"
                
  [1,] . . . . .
  [2,] . 1 . . .
  [3,] . . 1 . .
  [4,] . . . 1 .
  [5,] . . . . 1
  > M <- m1; M[,3] <- 3 ; stopifnot(is(M,"sparseMatrix"), M[,3] == 3)
  M[i,j] <- v :  coercing symmetric M[] into non-symmetric
  > Z <- m1; Z[] <- 0
  > checkMatrix(M)
   norm(m [5 x 5]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: Note: method with signature 'sparseMatrix#ldiMatrix' chosen for function '&',
   target signature 'nsCMatrix#ldiMatrix'.
   "nsparseMatrix#sparseMatrix", "nMatrix#lMatrix" would also be valid
  ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > M <- m1; M[1:3, 3] <- 0 ;M
  5 x 5 sparse Matrix of class "dtTMatrix"
                
  [1,] 1 . . . .
  [2,] . 1 . . .
  [3,] . . . . .
  [4,] . . . 1 .
  [5,] . . . . 1
  > assert.EQ.mat(M, diag(c(1,1, 0, 1,1)), tol=0)
  > T <- m1; T[1:3, 3] <- 10; checkMatrix(T)
   norm(m [5 x 5]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  as(<triangular (ge)matrix>, dtTMatrix): valid: TRUE 
  > stopifnot(is(T, "triangularMatrix"), identical(T[,3], c(10,10,10,0,0)),
  + 	  Qidentical(as(Z, "matrix"), z))
  > 
  > M <- m2; M[1,] <- 0 ; M ; assert.EQ.mat(M, diag(c(0,rep(1,4))), tol=0)
  5 x 5 sparse Matrix of class "dtCMatrix"
                
  [1,] . . . . .
  [2,] . 1 . . .
  [3,] . . 1 . .
  [4,] . . . 1 .
  [5,] . . . . 1
  > M <- m2; M[,3] <- 3 ; stopifnot(is(M,"sparseMatrix"), M[,3] == 3)
  > checkMatrix(M)
   norm(m [5 x 5]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > Z <- m2; Z[] <- 0
  > M <- m2; M[1:3, 3] <- 0 ;M
  5 x 5 sparse Matrix of class "dtCMatrix"
                
  [1,] 1 . . . .
  [2,] . 1 . . .
  [3,] . . . . .
  [4,] . . . 1 .
  [5,] . . . . 1
  > assert.EQ.mat(M, diag(c(1,1, 0, 1,1)), tol=0)
  > T <- m2; T[1:3, 3] <- 10; checkMatrix(T)
   norm(m [5 x 5]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  as(<triangular (ge)matrix>, dtCMatrix): valid: TRUE 
  > stopifnot(is(T, "dtCMatrix"), identical(T[,3], c(10,10,10,0,0)),
  + 	  Qidentical(as(Z, "matrix"), z))
  > showProc.time()
  Time elapsed:  1.57 0.02 1.61 
  > 
  > 
  > ## "Vector indices" -------------------
  > asLogical <- function(x) {
  +     stopifnot(is.atomic(x))
  +     storage.mode(x) <- "logical"
  +     x
  + }
  > .iniDiag.example <- expression({
  +     D <- Diagonal(6)
  +     M <- as(D,"dgeMatrix")
  +     m <- as(D,"matrix")
  +     s <- as(D,"TsparseMatrix"); N <- as(s,"nMatrix")
  +     S <- as(s,"CsparseMatrix"); C <- as(S,"nMatrix")
  + })
  > eval(.iniDiag.example)
  > i <- c(3,1,6); v <- c(10,15,20)
  > ## (logical,value) which both are recycled:
  > L <- c(TRUE, rep(FALSE,8)) ; z <- c(50,99)
  > 
  > ## vector subassignment, both with integer & logical
  > ## these now work correctly {though not very efficiently; hence warnings}
  > m[i] <- v # the role model: only first column is affected
  > M[i] <- v; assert.EQ.mat(M,m) # dge
  > D[i] <- v; assert.EQ.mat(D,m) # ddi -> dtT -> dgT
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  > s[i] <- v; assert.EQ.mat(s,m) # dtT -> dgT
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  > S[i] <- v; assert.EQ.mat(S,m); S # dtC -> dtT -> dgT -> dgC
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  6 x 6 sparse Matrix of class "dgCMatrix"
                   
  [1,] 15 . . . . .
  [2,]  . 1 . . . .
  [3,] 10 . 1 . . .
  [4,]  . . . 1 . .
  [5,]  . . . . 1 .
  [6,] 20 . . . . 1
  > m.L <- asLogical(m)
  > C[i] <- v # - with a warning: C is nMatrix, v not T/F
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
  Warning in `[<-`(`*tmp*`, i, value = c(10, 15, 20)) :
    x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
  > assert.EQ.mat(C,m.L); validObject(C)
  [1] TRUE
  > N[i] <- v # - with a warning
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
  Warning in `[<-`(`*tmp*`, i, value = c(10, 15, 20)) :
    x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
  > assert.EQ.mat(N,m.L); validObject(N)
  [1] TRUE
  > stopifnot(Q.C.identical(D,s, checkClass=FALSE))
  > ## logical *vector* indexing
  > eval(.iniDiag.example)
  > m[L] <- z; m.L <- asLogical(m)
  > M[L] <- z; assert.EQ.mat(M,m)
  > D[L] <- z; assert.EQ.mat(D,m)
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  > s[L] <- z; assert.EQ.mat(s,m)
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  > S[L] <- z; assert.EQ.mat(S,m) ; S
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  6 x 6 sparse Matrix of class "dgCMatrix"
                      
  [1,] 50  . . 50  . .
  [2,]  .  1 .  .  . .
  [3,]  .  . 1  .  . .
  [4,]  . 99 .  1 99 .
  [5,]  .  . .  .  1 .
  [6,]  .  . .  .  . 1
  > C[L] <- z; assert.EQ.mat(C,m.L) # with a good warning
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
  Warning in `[<-`(`*tmp*`, i, value = c(50, 99)) :
    x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
  > N[L] <- z; assert.EQ.mat(N,m.L)
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
  Warning in `[<-`(`*tmp*`, L, value = c(50, 99)) :
    x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
  > 
  > 
  > ## indexing [i]  vs  [i,] --- now ok
  > eval(.iniDiag.example)
  > stopifnot(identical5(m[i], M[i], D[i], s[i], S[i]), identical3(as.logical(m[i]), C[i], N[i]),
  +           identical5(m[L], M[L], D[L], s[L], S[L]), identical3(as.logical(m[L]), C[L], N[L]))
  <sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
  <sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
  <sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
  <sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
  <sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
  > ## bordercase ' drop = .' *vector* indexing {failed till 2009-04-..)
  > stopifnot(identical5(m[i,drop=FALSE], M[i,drop=FALSE], D[i,drop=FALSE],
  + 		     s[i,drop=FALSE], S[i,drop=FALSE]),
  + 	  identical3(as.logical(m[i,drop=FALSE]),
  + 		     C[i,drop=FALSE], N[i,drop=FALSE]))
  > stopifnot(identical5(m[L,drop=FALSE], M[L,drop=FALSE], D[L,drop=FALSE],
  + 		     s[L,drop=FALSE], S[L,drop=FALSE]),
  + 	  identical3(as.logical(m[L,drop=FALSE]),
  + 		     C[L,drop=FALSE], N[L,drop=FALSE]))
  > ## using L for row-indexing should give an error
  > assertError(m[L,]); assertError(m[L,, drop=FALSE])
  > ## these did not signal an error, upto (including) 0.999375-30:
  > assertError(s[L,]); assertError(s[L,, drop=FALSE])
  > assertError(S[L,]); assertError(S[L,, drop=FALSE])
  > assertError(N[L,]); assertError(N[L,, drop=FALSE])
  > 
  > ## row indexing:
  > assert.EQ.mat(D[i,], m[i,])
  > assert.EQ.mat(M[i,], m[i,])
  > assert.EQ.mat(s[i,], m[i,])
  > assert.EQ.mat(S[i,], m[i,])
  > assert.EQ.mat(C[i,], asLogical(m[i,]))
  > assert.EQ.mat(N[i,], asLogical(m[i,]))
  > ## column indexing:
  > assert.EQ.mat(D[,i], m[,i])
  > assert.EQ.mat(M[,i], m[,i])
  > assert.EQ.mat(s[,i], m[,i])
  > assert.EQ.mat(S[,i], m[,i])
  > assert.EQ.mat(C[,i], asLogical(m[,i]))
  > assert.EQ.mat(N[,i], asLogical(m[,i]))
  > 
  > 
  > ### --- negative indices ----------
  > 
  > ## 1) negative *vector* indexing
  > eval(.iniDiag.example)
  > i <- -(2:30)
  > stopifnot(identical5(m[i], M[i], D[i], s[i], S[i]),
  +           identical3(as.logical(m[i]), C[i], N[i]))
  > ##  negative vector subassignment :
  > v <- seq_along(m[i])
  > m[i] <- v; m.L <- asLogical(m)
  > M[i] <- v; assert.EQ.mat(M,m) # dge
  > D[i] <- v; assert.EQ.mat(D,m) # ddi -> dtT -> dgT
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  > s[i] <- v; assert.EQ.mat(s,m) # dtT -> dgT
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  > S[i] <- v; assert.EQ.mat(S,m); S # dtC -> dtT -> dgT -> dgC
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  6 x 6 sparse Matrix of class "dgCMatrix"
                  
  [1,] 1 . . . . 2
  [2,] . 1 . . . 3
  [3,] . . 1 . . 4
  [4,] . . . 1 . 5
  [5,] . . . . 1 6
  [6,] . . . . . 7
  > N[i] <- v # with a good warning
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
  Warning in `[<-`(`*tmp*`, i, value = 1:7) :
    x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
  > assert.EQ.mat(N,m.L); N
  6 x 6 sparse Matrix of class "ngTMatrix"
                  
  [1,] | . . . . |
  [2,] . | . . . |
  [3,] . . | . . |
  [4,] . . . | . |
  [5,] . . . . | |
  [6,] . . . . . |
  > C[i] <- v #  ..  .  ..  warning
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
  Warning in `[<-`(`*tmp*`, i, value = 1:7) :
    x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
  > assert.EQ.mat(C,m.L); C #
  6 x 6 sparse Matrix of class "ngCMatrix"
                  
  [1,] | . . . . |
  [2,] . | . . . |
  [3,] . . | . . |
  [4,] . . . | . |
  [5,] . . . . | |
  [6,] . . . . . |
  > 
  > options(warn = 2) #---------------------# NO WARNINGS from here -----------------
  > 					# =====================
  > ## 2) negative [i,j] indices
  > mc <- mC[1:5, 1:7]
  > mt <- mT[1:5, 1:7]
  > ## sub matrix
  > assert.EQ.mat(mC[1:2, 0:3], mm[1:2, 0:3]) # test 0-index
  > stopifnot(identical(mc[-(3:5), 0:2], mC[1:2, 0:2]),
  +           identical(mt[-(3:5), 0:2], mT[1:2, 0:2]),
  +           identical(mC[2:3, 4],      mm[2:3, 4]))
  > assert.EQ.mat(mC[1:2,], mm[1:2,])
  > ## sub vector
  > stopifnot(identical4(mc[-(1:4), ], mC[5, 1:7],
  +                      mt[-(1:4), ], mT[5, 1:7]))
  > stopifnot(identical4(mc[-(1:4), -(2:4)], mC[5, c(1,5:7)],
  +                      mt[-(1:4), -(2:4)], mT[5, c(1,5:7)]))
  > 
  > ## mixing of negative and positive must give error
  > assertError(mT[-1:1,])
  > showProc.time()
  Time elapsed:  0.64 0 0.64 
  > 
  > ## Sub *Assignment* ---- now works (partially):
  > mt0 <- mt
  > nt <- as(mt, "nMatrix")
  > mt[1, 4] <- -99
  > mt[2:3, 1:6] <- 0
  > mt
  5 x 7 sparse Matrix of class "dgTMatrix"
     c1 c2 c3  c4  c5 c6 c7
  r1  . 41 81 -99   .  .  .
  r2  .  .  .   .   .  .  .
  r3  .  .  .   .   .  .  .
  r4  . 44 84   .   .  .  .
  r5  . 45 85   . 165  .  .
  > m2 <- mt+mt
  > m2[1,4] <- -200
  > m2[c(1,3), c(5:6,2)] <- 1:6
  > stopifnot(m2[1,4] == -200,
  +           as.vector(m2[c(1,3), c(5:6,2)]) == 1:6)
  > mt[,3] <- 30
  > mt[2:3,] <- 250
  > mt[1:5 %% 2 == 1, 3] <- 0
  > mt[3:1, 1:7 > 5] <- 0
  > mt
  5 x 7 sparse Matrix of class "dgTMatrix"
      c1  c2  c3  c4  c5 c6 c7
  r1   .  41   . -99   .  .  .
  r2 250 250 250 250 250  .  .
  r3 250 250   . 250 250  .  .
  r4   .  44  30   .   .  .  .
  r5   .  45   .   . 165  .  .
  > 
  > tt <- as(mt,"matrix")
  > ii <- c(0,2,5)
  > jj <- c(2:3,5)
  > tt[ii, jj] <- 1:6 # 0 is just "dropped"
  > mt[ii, jj] <- 1:6
  > assert.EQ.mat(mt, tt)
  > 
  > mt[1:5, 2:6]
  5 x 5 sparse Matrix of class "dgTMatrix"
      c2 c3  c4  c5 c6
  r1  41  . -99   .  .
  r2   1  3 250   5  .
  r3 250  . 250 250  .
  r4  44 30   .   .  .
  r5   2  4   .   6  .
  > as((mt0 - mt)[1:5,], "dsparseMatrix")# [1,5] and lines 2:3
  5 x 7 sparse Matrix of class "dgCMatrix"
       c1   c2 c3   c4  c5 c6 c7
  r1    .    0 81  220   .  .  .
  r2 -250   41 -3 -250 157  .  .
  r3 -250 -250  . -250 -87  .  .
  r4    .    0 54    .   .  .  .
  r5    .   43 81    . 159  .  .
  > 
  > mt[c(2,4), ] <- 0; stopifnot(as(mt[c(2,4), ],"matrix") == 0)
  > mt[2:3, 4:7] <- 33
  > checkMatrix(mt)
   norm(m [5 x 7]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  > mt
  5 x 7 sparse Matrix of class "dgTMatrix"
      c1  c2 c3  c4 c5 c6 c7
  r1   .  41  . -99  .  .  .
  r2   .   .  .  33 33 33 33
  r3 250 250  .  33 33 33 33
  r4   .   .  .   .  .  .  .
  r5   .   2  4   .  6  .  .
  > 
  > mc[1,4] <- -99 ; stopifnot(mc[1,4] == -99)
  > mc[1,4] <-  00 ; stopifnot(mc[1,4] ==  00)
  > mc[1,4] <- -99 ; stopifnot(mc[1,4] == -99)
  > mc[1:2,4:3] <- 4:1; stopifnot(.asmatrix(mc[1:2,4:3]) == 4:1)
  > 
  > mc[-1, 3] <- -2:1 # 0 should not be entered; 'value' recycled
  > mt[-1, 3] <- -2:1
  > stopifnot(mc@x != 0, mt@x != 0,
  + 	  mc[-1,3] == -2:1, mt[-1,3] == -2:1) ## failed earlier
  > 
  > mc0 <- mc
  > mt0 <- as(mc0, "TsparseMatrix")
  > m0  <- as(mc0, "matrix")
  > set.seed(1); options(Matrix.verbose = FALSE)
  > for(i in 1:(if(doExtras) 50 else 4)) {
  +     mc <- mc0; mt <- mt0 ; m <- m0
  +     ev <- 1:5 %% 2 == round(runif(1))# 0 or 1
  +     j <- sample(ncol(mc), 1 + round(runif(1)))
  +     nv <- rpois(sum(ev) * length(j), lambda = 1)
  +     mc[ev, j] <- nv
  +      m[ev, j] <- nv
  +     mt[ev, j] <- nv
  +     if(i %% 10 == 1) print(mc[ev,j, drop = FALSE])
  +     stopifnot(as.vector(mc[ev, j]) == nv, ## failed earlier...
  +               as.vector(mt[ev, j]) == nv)
  +     validObject(mc) ; assert.EQ.mat(mc, m)
  +     validObject(mt) ; assert.EQ.mat(mt, m)
  + }
  2 x 1 sparse Matrix of class "dgCMatrix"
     c7
  r2  2
  r4  .
  > showProc.time()
  Time elapsed:  0.61 0 0.62 
  > options(Matrix.verbose = TRUE)
  > 
  > mc # no longer has non-structural zeros
  5 x 7 sparse Matrix of class "dgCMatrix"
     c1 c2 c3 c4  c5 c6 c7
  r1  . 41  2  4   .  .  .
  r2  . 42 -2  3 162  .  .
  r3  .  . -1  . 163  .  .
  r4  1 44  .  .   .  .  .
  r5  . 45  1  . 165  .  .
  > mc[ii, jj] <- 1:6
  > mc[c(2,5), c(3,5)] <- 3.2
  > checkMatrix(mc)
   norm(m [5 x 7]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  > m. <- mc
  > mc[4,] <- 0
  > mc
  5 x 7 sparse Matrix of class "dgCMatrix"
     c1 c2   c3 c4    c5 c6 c7
  r1  . 41  2.0  4   .    .  .
  r2  .  1  3.2  3   3.2  .  .
  r3  .  . -1.0  . 163.0  .  .
  r4  .  .  .    .   .    .  .
  r5  .  2  3.2  .   3.2  .  .
  > 
  > S <- as(Diagonal(5),"TsparseMatrix")
  > H <- Hilbert(9)
  > Hc <- as(round(H, 3), "dsCMatrix")# a sparse matrix with no 0 ...
  > (trH <- tril(Hc[1:5, 1:5]))
  5 x 5 sparse Matrix of class "dtCMatrix"
                                    
  [1,] 1.000 .     .     .     .    
  [2,] 0.500 0.333 .     .     .    
  [3,] 0.333 0.250 0.200 .     .    
  [4,] 0.250 0.200 0.167 0.143 .    
  [5,] 0.200 0.167 0.143 0.125 0.111
  > stopifnot(is(trH, "triangularMatrix"), trH@uplo == "L",
  +           is(S, "triangularMatrix"))
  > 
  > ## triangular assignment
  > ## the slick (but inefficient in case of sparse!) way to assign sub-diagonals:
  > ## equivalent to tmp <- `diag<-`(S[,-1], -2:1); S[,-1] <- tmp
  > ## which dispatches to (x="TsparseMatrix", i="missing",j="index", value="replValue")
  > diag(S[,-1]) <- -2:1 # used to give a wrong warning
  M[i,j] <- v :  coercing symmetric M[] into non-symmetric
  > S <- as(S,"triangularMatrix")
  > assert.EQ.mat(S, local({s <- diag(5); diag(s[,-1]) <- -2:1; s}))
  > 
  > trH[c(1:2,4), c(2:3,5)] <- 0 # gave an *error* upto Jan.2008
  > trH[ lower.tri(trH) ] <- 0   # ditto, because of callNextMethod()
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
  > 
  > m <- Matrix(0+1:28, nrow = 4)
  > m[-3,c(2,4:5,7)] <- m[ 3, 1:4] <- m[1:3, 6] <- 0
  > mT <- as(m, "dgTMatrix")
  > stopifnot(identical(mT[lower.tri(mT)],
  +                     m [lower.tri(m) ]))
  <sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
  > lM <- upper.tri(mT, diag=TRUE)
  > mT[lM] <- 0
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  >  m[lM] <- 0
  > assert.EQ.mat(mT, as(m,"matrix"))
  > mT[lM] <- -1:0
  diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
  >  m[lM] <- -1:0
  > assert.EQ.mat(mT, as(m,"matrix"))
  > (mT <- drop0(mT))
  4 x 7 sparse Matrix of class "dgCMatrix"
                           
  [1,] -1  .  . -1 -1 -1 -1
  [2,]  2 -1 -1  .  .  .  .
  [3,]  .  .  . -1 -1 -1 -1
  [4,]  4  . 12  .  .  .  .
  > 
  > i <- c(1:2, 4, 6:7); j <- c(2:4,6)
  > H[i,j] <- 0
  > (H. <- round(as(H, "sparseMatrix"), 3)[ , 2:7])
  9 x 6 sparse Matrix of class "dgCMatrix"
                                           
   [1,] .     .     .     0.200 .     0.143
   [2,] .     .     .     0.167 .     0.125
   [3,] 0.250 0.200 0.167 0.143 0.125 0.111
   [4,] .     .     .     0.125 .     0.100
   [5,] 0.167 0.143 0.125 0.111 0.100 0.091
   [6,] .     .     .     0.100 .     0.083
   [7,] .     .     .     0.091 .     0.077
   [8,] 0.111 0.100 0.091 0.083 0.077 0.071
   [9,] 0.100 0.091 0.083 0.077 0.071 0.067
  > Hc. <- Hc
  > Hc.[i,j] <- 0 ## now "works", but setting "non-structural" 0s
  > stopifnot(.asmatrix(Hc.[i,j]) == 0)
  > Hc.[, 1:6]
  9 x 6 sparse Matrix of class "dgCMatrix"
                                           
   [1,] 1.000 .     .     .     0.200 .    
   [2,] 0.500 .     .     .     0.167 .    
   [3,] 0.333 0.250 0.200 0.167 0.143 0.125
   [4,] 0.250 .     .     .     0.125 .    
   [5,] 0.200 0.167 0.143 0.125 0.111 0.100
   [6,] 0.167 .     .     .     0.100 .    
   [7,] 0.143 .     .     .     0.091 .    
   [8,] 0.125 0.111 0.100 0.091 0.083 0.077
   [9,] 0.111 0.100 0.091 0.083 0.077 0.071
  > 
  > ## an example that failed for a long time
  > sy3 <- new("dsyMatrix", Dim = as.integer(c(2, 2)), x = c(14, -1, 2, -7))
  > checkMatrix(dm <- kronecker(Diagonal(2), sy3))# now sparse with new kronecker
  Note: method with signature 'sparseMatrix#ANY' chosen for function 'kronecker',
   target signature 'dtTMatrix#dsyMatrix'.
   "ANY#Matrix" would also be valid
   norm(m [4 x 4]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > dm <- Matrix(.asmatrix(dm))# -> "dsyMatrix"
  > (s2 <- as(dm, "sparseMatrix"))
  4 x 4 sparse Matrix of class "dsCMatrix"
                  
  [1,] 14  2  .  .
  [2,]  2 -7  .  .
  [3,]  .  . 14  2
  [4,]  .  .  2 -7
  > checkMatrix(st <- as(s2, "TsparseMatrix"))
   norm(m [4 x 4]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: suboptimal 'Arith' implementation of  'dsC*  o  dsC*'
  identical
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: suboptimal 'Arith' implementation of  'dsC*  o  dsC*'
  ok;  determinant(): ok
  > stopifnot(is(s2, "symmetricMatrix"),
  + 	  is(st, "symmetricMatrix"))
  > checkMatrix(s.32  <- st[1:3,1:2]) ## 3 x 2 - and *not* dsTMatrix
   norm(m [3 x 2]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  > checkMatrix(s2.32 <- s2[1:3,1:2])
   norm(m [3 x 2]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  > I <- c(1,4:3)
  > stopifnot(is(s2.32, "generalMatrix"),
  +           is(s.32,  "generalMatrix"),
  +           identical(as.mat(s.32), as.mat(s2.32)),
  +           identical3(dm[1:3,-1], asD(s2[1:3,-1]), asD(st[1:3,-1])),
  +           identical4(2, dm[4,3], s2[4,3], st[4,3]),
  +           identical3(diag(dm), diag(s2), diag(st)),
  +           is((cI <- s2[I,I]), "dsCMatrix"),
  +           is((tI <- st[I,I]), "dsTMatrix"),
  +           identical4(as.mat(dm)[I,I], as.mat(dm[I,I]), as.mat(tI), as.mat(cI))
  +           )
  > 
  > ## now sub-assign  and check for consistency
  > ## symmetric subassign should keep symmetry
  > st[I,I] <- 0; checkMatrix(st); stopifnot(is(st,"symmetricMatrix"))
   norm(m [4 x 4]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: suboptimal 'Arith' implementation of  'dsC*  o  dsC*'
  identical
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: suboptimal 'Arith' implementation of  'dsC*  o  dsC*'
  ok;  determinant(): ok
  > s2[I,I] <- 0; checkMatrix(s2); stopifnot(is(s2,"symmetricMatrix"))
   norm(m [4 x 4]) : 1  I  F  M  ok
   Summary:  ok
  as(., "nMatrix") giving full nonzero-pattern: ok
  2*m =?= m+m: suboptimal 'Arith' implementation of  'dsC*  o  dsC*'
  identical
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: suboptimal 'Arith' implementation of  'dsC*  o  dsC*'
  ok;  determinant(): ok
  > ##
  > m <- as.mat(st)
  >  m[2:1,2:1] <- 4:1
  > st[2:1,2:1] <- 4:1
  M[i,j] <- v :  coercing symmetric M[] into non-symmetric
  > s2[2:1,2:1] <- 4:1
  > stopifnot(identical(m, as.mat(st)),
  + 	  1:4 == as.vector(s2[1:2,1:2]),
  + 	  identical(m, as.mat(s2)))
  > 
  > ## now a slightly different situation for 's2' (had bug)
  > s2 <- as(dm, "sparseMatrix")
  > s2[I,I] <- 0; diag(s2)[2:3] <- -(1:2)
  > stopifnot(is(s2,"symmetricMatrix"), diag(s2) == c(0:-2,0))
  > t2 <- as(s2, "TsparseMatrix")
  > m <- as.mat(s2)
  > s2[2:1,2:1] <- 4:1
  > t2[2:1,2:1] <- 4:1
  M[i,j] <- v :  coercing symmetric M[] into non-symmetric
  >  m[2:1,2:1] <- 4:1
  > assert.EQ.mat(t2, m)
  > assert.EQ.mat(s2, m)
  > ## and the same (for a different s2 !)
  > s2[2:1,2:1] <- 4:1
  > t2[2:1,2:1] <- 4:1
  > assert.EQ.mat(t2, m)# ok
  > assert.EQ.mat(s2, m)# failed in 0.9975-8
  > showProc.time()
  Time elapsed:  1.52 0 1.53 
  > 
  > 
  > ## m[cbind(i,j)] <- value: (2-column matrix subassignment):
  > m.[ cbind(3:5, 1:3) ] <- 1:3
  > stopifnot(m.[3,1] == 1, m.[4,2] == 2)
  > nt. <- nt ; nt[rbind(2:3, 3:4, c(3,3))] <- FALSE
  > s. <- m. ; m.[cbind(3,4:6)] <- 0 ## assigning 0 where there *is* 0 ..
  > stopifnot(identical(nt.,nt),       ## should not have changed
  + 	  identical(s., m.))
  Error: identical(s., m.) is not TRUE
  Execution halted
Running the tests in 'tests/matprod.R' failed.
Complete output:
  > library(Matrix)
  > 
  > ### Matrix Products including  cross products
  > 
  > source(system.file("test-tools.R", package = "Matrix")) # is.EQ.mat(), dnIdentical() ..etc
  Loading required package: tools
  > doExtras
  [1] FALSE
  > options(warn=1, # show as they happen
  + 	Matrix.verbose = doExtras)
  > 
  > ##' Check matrix multiplications with (unit) Diagonal matrices
  > chkDiagProd <- function(M) {
  +     stopifnot(is.matrix(M) || is(M,"Matrix"))
  +     I.l  <- Diagonal(nrow(M)) # I_n -- "unit" Diagonal
  +     I.r  <- Diagonal(ncol(M)) # I_d
  +     D2.l <- Diagonal(nrow(M), x = 2) # D_n
  +     D2.r <- Diagonal(ncol(M), x = 2) # I_d
  +     stopifnot(is.EQ.mat(M, M %*% I.r),
  + 	      is.EQ.mat(M, I.l %*% M),
  + 	      is.EQ.mat(2*M, M %*% D2.r),
  + 	      is.EQ.mat(M*2, D2.l %*% M),
  + 	      ## crossprod
  + 	      is.EQ.mat(t(M), crossprod(M, I.l)),
  + 	      is.EQ.mat(  M , crossprod(I.l, M)),
  + 	      is.EQ.mat(t(2*M), crossprod(M, D2.l)),
  + 	      is.EQ.mat(  M*2 , crossprod(D2.l, M)),
  + 	      ## tcrossprod
  + 	      is.EQ.mat(  M , tcrossprod(M, I.r)),
  + 	      is.EQ.mat(t(M), tcrossprod(I.r, M)),
  + 	      is.EQ.mat(  2*M , tcrossprod(M, D2.r)),
  + 	      is.EQ.mat(t(M*2), tcrossprod(D2.r, M)))
  + }
  > 
  > ### dimnames -- notably for matrix products ----------------
  > #
  > ##' Checks that matrix products are the same, including dimnames
  > ##'
  > ##' @param m matrix = traditional-R-matrix version of M
  > ##' @param M optional Matrix = "Matrix class version of m"
  > ##' @param browse
  > chkDnProd <- function(m = as(M, "matrix"), M = Matrix(m), browse=FALSE, warn.ok=FALSE) {
  +     ## TODO:
  +     ## if(browse) stopifnot <- f.unction(...)  such that it enters browser() when it is not fulfilled
  +     if(!warn.ok) { # NO warnings allowd
  +         op <- options(warn = 2)
  +         on.exit(options(op))
  +     }
  +     stopifnot(is.matrix(m), is(M, "Matrix"), identical(dim(m), dim(M)), dnIdentical(m,M))
  +     ## m is  n x d  (say)
  +     is.square <- nrow(m) == ncol(m)
  + 
  +     p1 <- (tm <- t(m)) %*% m ## d x d
  +     p1. <- crossprod(m)
  +     stopifnot(is.EQ.mat3(p1, p1., crossprod(m,m)))
  + 
  +     t1 <- m %*% tm ## n x n
  +     t1. <- tcrossprod(m)
  +     stopifnot(is.EQ.mat3(t1, t1., tcrossprod(m,m)))
  + 
  +     if(is.square) {
  + 	mm <- m %*% m
  + 	stopifnot(is.EQ.mat3(mm, crossprod(tm,m), tcrossprod(m,tm)))
  +     }
  + 
  +     chkDiagProd(m)## was not ok in Matrix 1.2.0
  + 
  +     ## Now the	"Matrix" ones -- should match the "matrix" above
  +     M0 <- M
  +     cat("sparse: ")
  +     for(sparse in c(TRUE, FALSE)) {
  + 	cat(sparse, "; ")
  + 	M <- as(M0, if(sparse) "sparseMatrix" else "denseMatrix")
  + 	P1 <- (tM <- t(M)) %*% M
  + 	P1. <- crossprod(M)
  + 	stopifnot(is.EQ.mat3(P1, P1., p1),
  + 		  is.EQ.mat3(P1., crossprod(M,M), crossprod(M,m)),
  + 		  is.EQ.mat (P1., crossprod(m,M)))
  + 
  + 	## P1. is "symmetricMatrix" -- semantically "must have" symm.dimnames
  + 	PP1 <- P1. %*% P1. ## still  d x d
  + 	R <- triu(PP1); r <- as(R, "matrix") # upper - triangular
  + 	L <- tril(PP1); l <- as(L, "matrix") # lower - triangular
  + 	stopifnot(isSymmetric(P1.), isSymmetric(PP1),
  + 		  isDiagonal(L) || is(L,"triangularMatrix"),
  + 		  isDiagonal(R) || is(R,"triangularMatrix"),
  + 		  isTriangular(L, upper=FALSE), isTriangular(R, upper=TRUE),
  + 		  is.EQ.mat(PP1, (pp1 <- p1 %*% p1)),
  + 		  dnIdentical(PP1, R),
  + 		  dnIdentical(L, R))
  + 
  + 	T1 <- M %*% tM
  + 	T1. <- tcrossprod(M)
  + 	stopifnot(is.EQ.mat3(T1, T1., t1),
  + 		  is.EQ.mat3(T1., tcrossprod(M,M), tcrossprod(M,m)),
  + 		  is.EQ.mat (T1., tcrossprod(m,M)),
  + 		  is.EQ.mat(tcrossprod(T1., tM),
  + 			    tcrossprod(t1., tm)),
  + 		  is.EQ.mat(crossprod(T1., M),
  + 			    crossprod(t1., m)))
  + 
  + 	## Now, *mixing*  Matrix x matrix:
  + 	stopifnot(is.EQ.mat3(tM %*% m, tm %*% M, tm %*% m))
  + 
  + 	if(is.square)
  + 	    stopifnot(is.EQ.mat (mm, M %*% M),
  + 		      is.EQ.mat3(mm, crossprod(tM,M), tcrossprod(M,tM)),
  + 		      ## "mixing":
  + 		      is.EQ.mat3(mm, crossprod(tm,M), tcrossprod(m,tM)),
  + 		      is.EQ.mat3(mm, crossprod(tM,m), tcrossprod(M,tm)))
  + 
  + 	## Symmetric and Triangular
  + 	stopifnot(is.EQ.mat(PP1 %*% tM, pp1 %*% tm),
  + 		  is.EQ.mat(R %*% tM, r %*% tm),
  + 		  is.EQ.mat(L %*% tM, L %*% tm))
  + 
  + 	## Diagonal :
  + 	chkDiagProd(M)
  +     }
  +     cat("\n")
  +     invisible(TRUE)
  + }
  > 
  > ## All these are ok  {now, (2012-06-11) also for dense
  > (m <- matrix(c(0, 0, 2:0), 3, 5))
       [,1] [,2] [,3] [,4] [,5]
  [1,]    0    1    0    0    2
  [2,]    0    0    2    0    1
  [3,]    2    0    1    0    0
  > m00 <- m # *no* dimnames
  > dimnames(m) <- list(LETTERS[1:3], letters[1:5])
  > (m.. <- m) # has *both* dimnames
    a b c d e
  A 0 1 0 0 2
  B 0 0 2 0 1
  C 2 0 1 0 0
  > m0. <- m.0 <- m..
  > dimnames(m0.)[1] <- list(NULL); m0.
       a b c d e
  [1,] 0 1 0 0 2
  [2,] 0 0 2 0 1
  [3,] 2 0 1 0 0
  > dimnames(m.0)[2] <- list(NULL); m.0
    [,1] [,2] [,3] [,4] [,5]
  A    0    1    0    0    2
  B    0    0    2    0    1
  C    2    0    1    0    0
  > d <- diag(3); dimnames(d) <- list(c("u","v","w"), c("X","Y","Z")); d
    X Y Z
  u 1 0 0
  v 0 1 0
  w 0 0 1
  > dU <- diagN2U(Matrix(d)) # unitriangular sparse
  > (T <- new("dtrMatrix", diag = "U", x= c(0,0,5,0), Dim= c(2L,2L),
  +           Dimnames= list(paste0("r",1:2),paste0("C",1:2)))) # unitriangular dense
  2 x 2 Matrix of class "dtrMatrix" (unitriangular)
     C1 C2
  r1  1  5
  r2  .  1
  > ##                                                                   ^^^^^^^^^^^^
  > ## currently many warnings about sub-optimal matrix products :
  > chkDnProd(m..)
  sparse: TRUE ; FALSE ; 
  > chkDnProd(m0.)
  sparse: TRUE ; FALSE ; 
  > chkDnProd(m.0)
  sparse: TRUE ; FALSE ; 
  > chkDnProd(m00)
  sparse: TRUE ; FALSE ; 
  > chkDnProd(M =   T)
  sparse: TRUE ; FALSE ; 
  > chkDnProd(M = t(T))
  sparse: TRUE ; FALSE ; 
  > chkDnProd(M =   dU)
  sparse: TRUE ; FALSE ; 
  > chkDnProd(M = t(dU))
  sparse: TRUE ; FALSE ; 
  > ## all the above failed in 1.2-0 and 1.1-5, 1.1-4 some even earlier
  > chkDnProd(M = Diagonal(4))
  sparse: TRUE ; FALSE ; 
  > chkDnProd(diag(x=3:1))
  sparse: TRUE ; FALSE ; 
  > chkDnProd(d)
  sparse: TRUE ; FALSE ; 
  > chkDnProd(M = as(d, "denseMatrix"))# M: dtrMatrix (diag = "N")
  sparse: TRUE ; FALSE ; 
  > 
  > 
  > m5 <- 1 + as(diag(-1:4)[-5,], "dgeMatrix")
  > ## named dimnames:
  > dimnames(m5) <- list(Rows= LETTERS[1:5], paste("C", 1:6, sep=""))
  > tr5 <- tril(m5[,-6])
  > m. <- as(m5, "matrix")
  > m5.2 <- local({t5 <- as.matrix(tr5); t5 %*% t5})
  > stopifnotValid(tr5, "dtrMatrix")
  > stopifnot(dim(m5) == 5:6,
  + 	  class(cm5 <- crossprod(m5)) == "dpoMatrix")
  > assert.EQ.mat(t(m5) %*% m5, as(cm5, "matrix"))
  > assert.EQ.mat(tr5.2 <- tr5 %*% tr5, m5.2)
  > stopifnotValid(tr5.2, "dtrMatrix")
  > stopifnot(as.vector(rep(1,6) %*% cm5) == colSums(cm5),
  + 	  as.vector(cm5 %*% rep(1,6)) == rowSums(cm5))
  > ## uni-diagonal dtrMatrix with "wrong" entries in diagonal
  > ## {the diagonal can be anything: because of diag = "U" it should never be used}:
  > tru <- Diagonal(3, x=3); tru[i.lt <- lower.tri(tru, diag=FALSE)] <- c(2,-3,4)
  > tru@diag <- "U" ; stopifnot(diag(trm <- as.matrix(tru)) == 1)
  > ## TODO: Also add *upper-triangular*  *packed* case !!
  > stopifnot((tru %*% tru)[i.lt] ==
  + 	  (trm %*% trm)[i.lt])
  > 
  > ## crossprod() with numeric vector RHS and LHS
  > i5 <- rep.int(1, 5)
  > isValid(S5 <- tcrossprod(tr5), "dpoMatrix")# and inherits from "dsy"
  [1] TRUE
  > G5 <- as(S5, "generalMatrix")# "dge"
  > assert.EQ.mat( crossprod(i5, m5), rbind( colSums(m5)))
  > assert.EQ.mat( crossprod(i5, m.), rbind( colSums(m5)))
  > assert.EQ.mat( crossprod(m5, i5), cbind( colSums(m5)))
  > assert.EQ.mat( crossprod(m., i5), cbind( colSums(m5)))
  > assert.EQ.mat( crossprod(i5, S5), rbind( colSums(S5))) # failed in Matrix 1.1.4
  > ## tcrossprod() with numeric vector RHS and LHS :
  > stopifnot(identical(tcrossprod(i5, S5), # <- lost dimnames
  + 		    tcrossprod(i5, G5) -> m51),
  + 	  identical(dimnames(m51), list(NULL, LETTERS[1:5]))
  + 	  )
  > m51 <- m5[, 1, drop=FALSE] # [6 x 1]
  > m.1 <- m.[, 1, drop=FALSE] ; assert.EQ.mat(m51, m.1)
  > ## The only (M . v) case
  > assert.EQ.mat(tcrossprod(m51, 5:1), tcrossprod(m.1, 5:1))
  > ## The two  (v . M) cases:
  > assert.EQ.mat(tcrossprod(rep(1,6), m.), rbind( rowSums(m5)))# |v| = ncol(m)
  > assert.EQ.mat(tcrossprod(rep(1,3), m51),
  +               tcrossprod(rep(1,3), m.1))# ncol(m) = 1
  > 
  > ## classes differ
  > tc.m5 <- m5 %*% t(m5)	 # "dge*", no dimnames (FIXME)
  > (tcm5 <- tcrossprod(m5)) # "dpo*"  w/ dimnames
  5 x 5 Matrix of class "dpoMatrix"
    A  B  C  D  E
  A 5  5  6  7  9
  B 5  6  7  8 10
  C 6  7  9  9 11
  D 7  8  9 14 12
  E 9 10 11 12 30
  > assert.EQ.mat(tc.m5, mm5 <- as(tcm5, "matrix"))
  > ## tcrossprod(x,y) :
  > assert.EQ.mat(tcrossprod(m5, m5), mm5)
  > assert.EQ.mat(tcrossprod(m5, m.), mm5)
  > assert.EQ.mat(tcrossprod(m., m5), mm5)
  > 
  > M50 <- m5[,FALSE, drop=FALSE]
  > M05 <- t(M50)
  > s05 <- as(M05, "sparseMatrix")
  > s50 <- t(s05)
  > assert.EQ.mat(M05, matrix(1, 0,5))
  > assert.EQ.mat(M50, matrix(1, 5,0))
  > assert.EQ.mat(tcrossprod(M50), tcrossprod(as(M50, "matrix")))
  > assert.EQ.mat(tcrossprod(s50), tcrossprod(as(s50, "matrix")))
  > assert.EQ.mat( crossprod(s50),	crossprod(as(s50, "matrix")))
  > stopifnot(identical( crossprod(s50), tcrossprod(s05)),
  + 	  identical( crossprod(s05), tcrossprod(s50)))
  > (M00 <- crossprod(M50))## used to fail -> .Call(dgeMatrix_crossprod, x, FALSE)
  0 x 0 Matrix of class "dpoMatrix"
  <0 x 0 matrix>
  > stopifnot(identical(M00, tcrossprod(M05)),
  + 	  all(M00 == t(M50) %*% M50), dim(M00) == 0)
  > 
  > ## simple cases with 'scalars' treated as 1x1 matrices:
  > d <- Matrix(1:5)
  > d %*% 2
  5 x 1 Matrix of class "dgeMatrix"
       [,1]
  [1,]    2
  [2,]    4
  [3,]    6
  [4,]    8
  [5,]   10
  > 10 %*% t(d)
  1 x 5 Matrix of class "dgeMatrix"
       [,1] [,2] [,3] [,4] [,5]
  [1,]   10   20   30   40   50
  > assertError(3 %*% d)		 # must give an error , similar to
  > assertError(5 %*% as.matrix(d))	 # -> error
  > 
  > ## right and left "numeric" and "matrix" multiplication:
  > (p1 <- m5 %*% c(10, 2:6))
  5 x 1 Matrix of class "dgeMatrix"
    [,1]
  A   20
  B   30
  C   33
  D   38
  E   54
  > (p2 <- c(10, 2:5) %*% m5)
  1 x 6 Matrix of class "dgeMatrix"
       C1 C2 C3 C4 C5 C6
  [1,] 14 24 27 32 24 44
  > (pd1 <- m5 %*% diag(1:6))
  5 x 6 Matrix of class "dgeMatrix"
    [,1] [,2] [,3] [,4] [,5] [,6]
  A    0    2    3    4    5    6
  B    1    2    3    4    5    6
  C    1    2    6    4    5    6
  D    1    2    3   12    5    6
  E    1    2    3    4    5   30
  > (pd. <- m5 %*% Diagonal(x = 1:6))
  5 x 6 Matrix of class "dgeMatrix"
      
  Rows C1 C2 C3 C4 C5 C6
     A  0  2  3  4  5  6
     B  1  2  3  4  5  6
     C  1  2  6  4  5  6
     D  1  2  3 12  5  6
     E  1  2  3  4  5 30
  > (pd2 <- diag (10:6)	   %*% m5)
  5 x 6 Matrix of class "dgeMatrix"
       C1 C2 C3 C4 C5 C6
  [1,]  0 10 10 10 10 10
  [2,]  9  9  9  9  9  9
  [3,]  8  8 16  8  8  8
  [4,]  7  7  7 21  7  7
  [5,]  6  6  6  6  6 30
  > (pd..<- Diagonal(x = 10:6) %*% m5)
  5 x 6 Matrix of class "dgeMatrix"
      
  Rows C1 C2 C3 C4 C5 C6
     A  0 10 10 10 10 10
     B  9  9  9  9  9  9
     C  8  8 16  8  8  8
     D  7  7  7 21  7  7
     E  6  6  6  6  6 30
  > stopifnot(dim(crossprod(t(m5))) == c(5,5),
  + 	  c(class(p1),class(p2),class(pd1),class(pd2),
  + 	    class(pd.),class(pd..)) == "dgeMatrix")
  > assert.EQ.mat(p1, cbind(c(20,30,33,38,54)))
  > assert.EQ.mat(pd1, m. %*% diag(1:6))
  > assert.EQ.mat(pd2, diag(10:6) %*% m.)
  > assert.EQ.mat(pd., as(pd1,"matrix"))
  > assert.EQ.mat(pd..,as(pd2,"matrix"))
  > 
  > ## check that 'solve' and '%*%' are inverses
  > set.seed(1)
  > A <- Matrix(rnorm(25), nc = 5)
  > y <- rnorm(5)
  > all.equal((A %*% solve(A, y))@x, y)
  [1] TRUE
  > Atr <- new("dtrMatrix", Dim = A@Dim, x = A@x, uplo = "U")
  > all.equal((Atr %*% solve(Atr, y))@x, y)
  [1] TRUE
  > 
  > ## R-forge bug 5933 by Sebastian Herbert,
  > ## https://r-forge.r-project.org/tracker/index.php?func=detail&aid=5933&group_id=61&atid=294
  > mLeft  <- matrix(data = double(0), nrow = 3, ncol = 0)
  > mRight <- matrix(data = double(0), nrow = 0, ncol = 4)
  > MLeft   <- Matrix(data = double(0), nrow = 3, ncol = 0)
  > MRight  <- Matrix(data = double(0), nrow = 0, ncol = 4)
  > stopifnot(identical3(class(mLeft), class(mRight), "matrix"))
  > stopifnot(class(MLeft) ==  class(MRight),
  +           class(MLeft) ==  "dgeMatrix")
  > 
  > Qidentical3 <- function(a,b,c) Q.eq(a,b) && Q.eq(b,c)
  > Qidentical4 <- function(a,b,c,d) Q.eq(a,b) && Q.eq(b,c) && Q.eq(c,d)
  > chkP <- function(mLeft, mRight, MLeft, MRight, cl = class(MLeft)) {
  +     ident4 <- if(extends(cl, "generalMatrix"))
  + 		  function(a,b,c,d) identical4(as(a, cl), b, c, d)
  + 	      else function(a,b,c,d) {
  + 		  assert.EQ.mat(M=b, m=a, tol=0)
  + 		  Qidentical3(b,c,d)
  + 	      }
  +     mm <- mLeft %*% mRight # ok
  +     m.m <- crossprod(mRight)
  +     mm. <- tcrossprod(mLeft, mLeft)
  +     stopifnot(mm == 0,
  +               ident4(mm,
  +                      mLeft %*% MRight,
  +                      MLeft %*% mRight,
  +                      MLeft %*% MRight),# now ok
  + 	      m.m == 0, identical(m.m, crossprod(mRight, mRight)),
  + 	      mm. == 0, identical(mm., tcrossprod(mLeft, mLeft)))
  +     stopifnot(ident4(m.m,
  + 		     crossprod(MRight, MRight),
  + 		     crossprod(MRight, mRight),
  + 		     crossprod(mRight, MRight)))
  +     stopifnot(ident4(mm.,
  + 		     tcrossprod(mLeft, MLeft),
  + 		     tcrossprod(MLeft, MLeft),
  + 		     tcrossprod(MLeft, mLeft)))
  + }
  > 
  > chkP(mLeft, mRight, MLeft, MRight, "dgeMatrix")
  > m0 <- mLeft[FALSE,] # 0 x 0
  > for(cls in c("triangularMatrix", "symmetricMatrix")) {
  +     cat(cls,": "); stopifnotValid(ML0 <- as(MLeft[FALSE,], cls), cls)
  +     chkP(m0, mRight, ML0, MRight, class(ML0))
  +     chkP(mLeft, m0 , MLeft, ML0 , class(ML0))
  +     chkP(m0,    m0 , ML0,   ML0 , class(ML0)); cat("\n")
  + }
  triangularMatrix : 
  symmetricMatrix : 
  > 
  > 
  > ## New in R 3.2.0 -- for traditional matrix m and vector v
  > for(spV in c(FALSE,TRUE)) {
  +     cat("sparseV:", spV, "\n")
  +     v <- if(spV) as(1:3, "sparseVector") else 1:3
  +     stopifnot(identical(class(v2 <- v[1:2]), class(v)))
  +     assertError(crossprod(v, v2)) ; assertError(v %*% v2)
  +     assertError(crossprod(v, 1:2)); assertError(v %*% 1:2)
  +     assertError(crossprod(v, 2))  ; assertError(v %*% 2)
  +     assertError(crossprod(1:2, v)); assertError(1:2 %*% v)
  +     if(spV || getRversion() >= "3.2.0") {
  + 	cat("doing  vec x vec  ..\n")
  + 	stopifnot(identical(crossprod(2, v), t(2) %*% v),
  + 		  identical(5 %*% v, 5 %*% t(v)))
  +     }
  +     for(sp in c(FALSE, TRUE)) {
  +         m <- Matrix(1:2, 1,2, sparse=sp)
  +         cat(sprintf("class(m): '%s'\n", class(m)))
  + 	stopifnot(identical( crossprod(m, v), t(m) %*% v), # m'v gave *outer* prod wrongly!
  + 		  identical(tcrossprod(m, v2), m %*% v2))
  + 	assert.EQ.Mat(m %*% v2, m %*% 1:2, tol=0)
  +     }
  +     ## gave error "non-conformable arguments"
  + }
  sparseV: FALSE 
  doing  vec x vec  ..
  class(m): 'dgeMatrix'
  class(m): 'dgCMatrix'
  sparseV: TRUE 
  doing  vec x vec  ..
  class(m): 'dgeMatrix'
  class(m): 'dgCMatrix'
  > ##  crossprod(m, v)  t(1 x 2) * 3 ==> (2 x 1) *  (1 x 3) ==> 2 x 3
  > ## tcrossprod(m,v2)    1 x 2  * 2 ==> (1 x 2) * t(1 x 2) ==> 1 x 1
  > 
  > ### ------ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  > ### Sparse Matrix products
  > ### ------
  > ## solve() for dtC*
  > mc <- round(chol(crossprod(A)), 2)
  > B <- A[1:3,] # non-square on purpose
  > stopifnot(all.equal(sum(rowSums(B %*% mc)), 5.82424475145))
  > assert.EQ.mat(tcrossprod(B, mc), as.matrix(t(tcrossprod(mc, B))))
  > 
  > m <- kronecker(Diagonal(2), mc)
  > stopifnot(is(mc, "Cholesky"),
  + 	  is(m, "sparseMatrix"))
  > im <- solve(m)
  > round(im, 3)
  10 x 10 sparse Matrix of class "dtCMatrix"
                                                                           
   [1,] 0.515 -0.161 0.412 -0.095   2.126 .      .     .      .       .    
   [2,] .      0.763 0.896 -0.473  -9.344 .      .     .      .       .    
   [3,] .      .     0.645 -0.028 -10.444 .      .     .      .       .    
   [4,] .      .     .      0.862   1.533 .      .     .      .       .    
   [5,] .      .     .      .      11.111 .      .     .      .       .    
   [6,] .      .     .      .       .     0.515 -0.161 0.412 -0.095   2.126
   [7,] .      .     .      .       .     .      0.763 0.896 -0.473  -9.344
   [8,] .      .     .      .       .     .      .     0.645 -0.028 -10.444
   [9,] .      .     .      .       .     .      .     .      0.862   1.533
  [10,] .      .     .      .       .     .      .     .      .      11.111
  > itm <- solve(t(m))
  > iim <- solve(im) # should be ~= 'm' of course
  > iitm <- solve(itm)
  > I <- Diagonal(nrow(m))
  > (del <- c(mean(abs(as.numeric(im %*% m - I))),
  + 	  mean(abs(as.numeric(m %*% im - I))),
  + 	  mean(abs(as.numeric(im  - t(itm)))),
  + 	  mean(abs(as.numeric( m  - iim))),
  + 	  mean(abs(as.numeric(t(m)- iitm)))))
  [1] 2.479665e-17 8.147359e-17 1.476597e-16 3.275158e-17 1.998401e-17
  > stopifnot(is(m, "triangularMatrix"), is(m, "sparseMatrix"),
  + 	  is(im, "dtCMatrix"), is(itm, "dtCMatrix"), is(iitm, "dtCMatrix"),
  + 	  del < 1e-15)
  > 
  > ## crossprod(.,.) & tcrossprod(),  mixing dense & sparse
  > v <- c(0,0,2:0)
  > mv <- as.matrix(v) ## == cbind(v)
  > (V <- Matrix(v, 5,1, sparse=TRUE))
  5 x 1 sparse Matrix of class "dgCMatrix"
        
  [1,] .
  [2,] .
  [3,] 2
  [4,] 1
  [5,] .
  > sv <- as(v, "sparseVector")
  > a <- as.matrix(A)
  > cav <-	crossprod(a,v)
  > tva <- tcrossprod(v,a)
  > assert.EQ.mat(crossprod(A, V), cav) # gave infinite recursion
  > assert.EQ.mat(crossprod(A,sv), cav)
  > assert.EQ.mat(tcrossprod( sv, A), tva)
  > assert.EQ.mat(tcrossprod(t(V),A), tva)
  > 
  > ## [t]crossprod()  for  <sparsevector> . <sparsevector>  incl. one arg.:
  > stopifnotValid(s.s <-  crossprod(sv,sv), "Matrix")
  > stopifnotValid(ss. <- tcrossprod(sv,sv), "sparseMatrix")
  > stopifnot(identical(s.s,  crossprod(sv)),
  + 	  identical(ss., tcrossprod(sv)))
  > assert.EQ.mat(s.s,  crossprod(v,v))
  > assert.EQ.mat(ss., tcrossprod(v,v))
  > 
  > dm <- Matrix(v, sparse=FALSE)
  > sm <- Matrix(v, sparse=TRUE)
  > 
  > validObject(d.vvt <- as(as(vvt <- tcrossprod(v, v), "denseMatrix"), "dgeMatrix"))
  [1] TRUE
  > validObject(s.vvt <- as(as(vvt, "sparseMatrix"), "dgCMatrix"))
  [1] TRUE
  > 
  > stopifnot(
  +     identical4(tcrossprod(v, v), tcrossprod(mv, v),
  + 	       tcrossprod(v,mv), tcrossprod(mv,mv))## (base R)
  +     ,
  +     identical4(d.vvt,		 tcrossprod(dm, v),
  + 	       tcrossprod(v,dm), tcrossprod(dm,dm)) ## (v, dm) failed
  +     ,
  +     identical(s.vvt, tcrossprod(sm,sm))
  +     ,
  +     identical3(d.vvt, tcrossprod(sm, v),
  +                tcrossprod(v,sm)) ## both (sm,v) and (v,sm) failed
  +     )
  > assert.EQ.mat(d.vvt, vvt)
  > assert.EQ.mat(s.vvt, vvt)
  > 
  > M <- Matrix(0:5, 2,3) ; sM <- as(M, "sparseMatrix"); m <- as(M, "matrix")
  > v <- 1:3; v2 <- 2:1
  > sv  <- as( v, "sparseVector")
  > sv2 <- as(v2, "sparseVector")
  > tvm <- tcrossprod(v, m)
  > assert.EQ.mat(tcrossprod( v, M), tvm)
  > assert.EQ.mat(tcrossprod( v,sM), tvm)
  > assert.EQ.mat(tcrossprod(sv,sM), tvm)
  > assert.EQ.mat(tcrossprod(sv, M), tvm)
  > assert.EQ.mat(crossprod(M, sv2), crossprod(m, v2))
  > stopifnot(identical(tcrossprod(v, M), v %*% t(M)),
  + 	  identical(tcrossprod(v,sM), v %*% t(sM)),
  + 	  identical(tcrossprod(v, M), crossprod(v, t(M))),
  + 	  identical(tcrossprod(sv,sM), sv %*% t(sM)),
  + 	  identical(crossprod(sM, sv2), t(sM) %*% sv2),
  + 	  identical(crossprod(M, v2), t(M) %*% v2))
  > 
  > ## *unit* triangular :
  > t1 <- new("dtTMatrix", x= c(3,7), i= 0:1, j=3:2, Dim= as.integer(c(4,4)))
  > ## from	 0-diagonal to unit-diagonal {low-level step}:
  > tu <- t1 ; tu@diag <- "U"
  > cu <- as(tu, "dtCMatrix")
  > cl <- t(cu) # unit lower-triangular
  > cl10 <- cl %*% Diagonal(4, x=10)
  > assert.EQ.mat(cl10, as(cl, "matrix") %*% diag(4, x=10))
  > stopifnot(is(cl,"dtCMatrix"), cl@diag == "U")
  > (cu2 <- cu %*% cu)
  4 x 4 sparse Matrix of class "dtCMatrix" (unitriangular)
               
  [1,] I .  . 6
  [2,] . I 14 .
  [3,] . .  I .
  [4,] . .  . I
  > cl2 <- cl %*% cl
  > validObject(cl2)
  [1] TRUE
  > 
  > cu3 <- tu[-1,-1]
  > assert.EQ.mat(crossprod(tru, cu3),## <- "FIXME" should return triangular ...
  + 	      crossprod(trm, as.matrix(cu3)))
  > 
  > cl2
  4 x 4 sparse Matrix of class "dtCMatrix" (unitriangular)
               
  [1,] I  . . .
  [2,] .  I . .
  [3,] . 14 I .
  [4,] 6  . . I
  > mcu <- as.matrix(cu)
  > cu2. <- Diagonal(4) + Matrix(c(rep(0,9),14,0,0,6,0,0,0), 4,4)
  > D4 <- Diagonal(4, x=10:7); d4 <- as(D4, "matrix")
  > D.D4 <- crossprod(D4); assert.EQ.mat(D.D4, crossprod(d4))
  > stopifnotValid(D.D4, "ddiMatrix")
  > stopifnotValid(su <- crossprod(cu), "dsCMatrix")
  > stopifnot(
  +     all(cu2 == cu2.)
  +    ,# was wrong for ver. <= 0.999375-4
  +           identical(D.D4, tcrossprod(D4))
  + ,
  +           identical4(crossprod(d4, D4), crossprod(D4, d4), tcrossprod(d4, D4), D.D4)
  + ,
  + 	  is(cu2, "dtCMatrix"), is(cl2, "dtCMatrix"), # triangularity preserved
  + 	  cu2@diag == "U", cl2@diag == "U",# UNIT-triangularity preserved
  + 	  all.equal(D4 %*% cu, D4 %*% mcu)
  + ,
  + 	  all.equal(cu %*% D4, mcu %*% D4)
  + ,
  + 	  all(D4 %*% su == D4 %*% as.mat(su))
  + ,
  + 	  all(su %*% D4 == as.mat(su) %*% D4)
  + ,
  + 	  identical(t(cl2), cu2)
  +   , # !!
  +           identical ( crossprod(cu, D4),  crossprod(mcu, D4))
  + ,
  +           identical4(tcrossprod(cu, D4), tcrossprod(mcu, D4), cu %*% D4, mcu %*% D4)
  + ,
  +           identical4(tcrossprod(D4,cu), tcrossprod(D4,mcu), D4 %*% t(cu), D4 %*% t(mcu))
  + ,
  + 	  identical( crossprod(cu), Matrix( crossprod(mcu),sparse=TRUE))
  + ,
  + 	  identical(tcrossprod(cu), Matrix(tcrossprod(mcu),sparse=TRUE))
  +     )
  > assert.EQ.mat( crossprod(cu, D4),  crossprod(mcu, d4))
  > assert.EQ.mat(tcrossprod(cu, D4), tcrossprod(mcu, d4))
  > tr8 <- kronecker(rbind(c(2,0),c(1,4)), cl2)
  > T8 <- tr8 %*% (tr8/2) # triangularity preserved?
  > T8.2 <- (T8 %*% T8) / 4
  > stopifnot(is(T8, "triangularMatrix"), T8@uplo == "L", is(T8.2, "dtCMatrix"))
  > mr8 <- as(tr8,"matrix")
  > m8. <- (mr8 %*% mr8 %*% mr8 %*% mr8)/16
  > assert.EQ.mat(T8.2, m8.)
  > 
  > data(KNex); mm <- KNex$mm
  > M <- mm[1:500, 1:200]
  > MT <- as(M, "TsparseMatrix")
  > cpr   <- t(mm) %*% mm
  > cpr.  <- crossprod(mm)
  > cpr.. <- crossprod(mm, mm)
  > stopifnot(is(cpr., "symmetricMatrix"),
  + 	  identical3(cpr, as(cpr., class(cpr)), cpr..))
  > ## with dimnames:
  > m <- Matrix(c(0, 0, 2:0), 3, 5)
  > dimnames(m) <- list(LETTERS[1:3], letters[1:5])
  > m
  3 x 5 sparse Matrix of class "dgCMatrix"
    a b c d e
  A . 1 . . 2
  B . . 2 . 1
  C 2 . 1 . .
  > p1 <- t(m) %*% m
  > (p1. <- crossprod(m))
  5 x 5 sparse Matrix of class "dsCMatrix"
    a b c d e
  a 4 . 2 . .
  b . 1 . . 2
  c 2 . 5 . 2
  d . . . . .
  e . 2 2 . 5
  > t1 <- m %*% t(m)
  > (t1. <- tcrossprod(m))
  3 x 3 sparse Matrix of class "dsCMatrix"
    A B C
  A 5 2 .
  B 2 5 2
  C . 2 5
  > stopifnot(isSymmetric(p1.),
  + 	  isSymmetric(t1.),
  + 	  identical(p1, as(p1., class(p1))),
  + 	  identical(t1, as(t1., class(t1))),
  + 	  identical(dimnames(p1), dimnames(p1.)),
  + 	  identical(dimnames(p1), list(colnames(m), colnames(m))),
  + 	  identical(dimnames(t1), dimnames(t1.))
  + 	  )
  > 
  > showMethods("%*%", class=class(M))
  Function: %*% (package base)
  x="ddiMatrix", y="dgCMatrix"
      (inherited from: x="diagonalMatrix", y="CsparseMatrix")
  x="dgCMatrix", y="ddiMatrix"
      (inherited from: x="CsparseMatrix", y="diagonalMatrix")
  x="dgCMatrix", y="dgCMatrix"
      (inherited from: x="CsparseMatrix", y="CsparseMatrix")
  x="dgCMatrix", y="dgTMatrix"
      (inherited from: x="Matrix", y="TsparseMatrix")
  x="dgCMatrix", y="dtCMatrix"
      (inherited from: x="CsparseMatrix", y="CsparseMatrix")
  x="dgCMatrix", y="integer"
      (inherited from: x="CsparseMatrix", y="numLike")
  x="dgCMatrix", y="isparseVector"
      (inherited from: x="mMatrix", y="sparseVector")
  x="dgCMatrix", y="matrix"
      (inherited from: x="CsparseMatrix", y="matrix")
  x="dgTMatrix", y="dgCMatrix"
      (inherited from: x="TsparseMatrix", y="Matrix")
  x="dgeMatrix", y="dgCMatrix"
      (inherited from: x="ddenseMatrix", y="CsparseMatrix")
  x="dtCMatrix", y="dgCMatrix"
      (inherited from: x="CsparseMatrix", y="CsparseMatrix")
  x="integer", y="dgCMatrix"
      (inherited from: x="numLike", y="CsparseMatrix")
  x="isparseVector", y="dgCMatrix"
      (inherited from: x="sparseVector", y="mMatrix")
  x="matrix", y="dgCMatrix"
      (inherited from: x="matrix", y="CsparseMatrix")
  
  > 
  > v1 <- rep(1, ncol(M))
  > str(r <-  M %*% Matrix(v1))
  Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
    ..@ x       : num [1:500] 0.277 0.5 0.277 0.5 0.5 ...
    ..@ Dim     : int [1:2] 500 1
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ factors : list()
  > str(rT <- MT %*% Matrix(v1))
  Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
    ..@ x       : num [1:500] 0.277 0.5 0.277 0.5 0.5 ...
    ..@ Dim     : int [1:2] 500 1
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ factors : list()
  > stopifnot(identical(r, rT))
  > str(r. <- M %*% as.matrix(v1))
  Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
    ..@ x       : num [1:500] 0.277 0.5 0.277 0.5 0.5 ...
    ..@ Dim     : int [1:2] 500 1
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ factors : list()
  > stopifnot(identical4(r, r., rT, M %*% as(v1, "matrix")))
  > 
  > v2 <- rep(1,nrow(M))
  > r2 <- t(Matrix(v2)) %*% M
  > r2T <- v2 %*% MT
  > str(r2. <- v2 %*% M)
  Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
    ..@ x       : num [1:200] 1.66 2 1.33 1.15 1.79 ...
    ..@ Dim     : int [1:2] 1 200
    ..@ Dimnames:List of 2
    .. ..$ : NULL
    .. ..$ : NULL
    ..@ factors : list()
  > stopifnot(identical3(r2, r2., t(as(v2, "matrix")) %*% M))
  > 
  > 
  > ###------------------------------------------------------------------
  > ### Close to singular matrix W
  > ### (from  micEconAids/tests/aids.R ... priceIndex = "S" )
  > (load(system.file("external", "symW.rda", package="Matrix"))) # "symW"
  [1] "symW"
  > stopifnot(is(symW, "symmetricMatrix"))
  > n <- nrow(symW)
  > I <- .sparseDiagonal(n, shape="g")
  > S <- as(symW, "matrix")
  > sis <- solve(S,    S)
  > ## solve(<dsCMatrix>, <sparseMatrix>) when Cholmod fails was buggy for *long*:
  > SIS <- solve(symW, symW)
  > iw <- solve(symW)
  > iss <- iw %*% symW
  > ##                                     nb-mm3   openBLAS (Avi A.)
  > assert.EQ.(I, drop0(sis), tol = 1e-8)# 2.6e-10;  7.96e-9
  Mean relative difference: 2.611773e-10
  > assert.EQ.(I, SIS,        tol = 1e-7)# 8.2e-9
  Mean relative difference: 7.505174e-09
  > assert.EQ.(I, iss,        tol = 4e-4)# 3.3e-5
  Mean relative difference: 5.137495e-05
  > ## solve(<dsCMatrix>, <dense..>) :
  > I <- diag(nr=n)
  > SIS <- solve(symW, as(symW,"denseMatrix"))
  > iw  <- solve(symW, I)
  > iss <- iw %*% symW
  > assert.EQ.mat(SIS, I, tol = 1e-7, giveRE=TRUE)
  Mean relative difference: 4.054743e-09
  > assert.EQ.mat(iss, I, tol = 4e-4, giveRE=TRUE)
  Mean relative difference: 2.443711e-05
  > rm(SIS,iss)
  > 
  > WW <- as(symW, "generalMatrix") # the one that gave problems
  > IW <- solve(WW)
  > class(I1 <- IW %*% WW)# "dge" or "dgC" (!)
  [1] "dgCMatrix"
  attr(,"package")
  [1] "Matrix"
  > class(I2 <- WW %*% IW)
  [1] "dgCMatrix"
  attr(,"package")
  [1] "Matrix"
  > ## these two were wrong for for M.._1.0-13:
  > assert.EQ.(as(I1,"matrix"), I, tol = 1e-4)
  Mean relative difference: 5.137495e-05
  > assert.EQ.(as(I2,"matrix"), I, tol = 7e-7)
  Mean relative difference: 5.131231e-09
  > 
  > ## now slightly perturb WW (and hence break exact symmetry
  > set.seed(131); ii <- sample(length(WW), size= 100)
  > WW[ii] <- WW[ii] * (1 + 1e-7*runif(100))
  > SW. <- symmpart(WW)
  > SW2 <- Matrix:::forceSymmetric(WW)
  > stopifnot(all.equal(as(SW.,"matrix"),
  +                     as(SW2,"matrix"), tol = 1e-7))
  > (ch <- all.equal(WW, as(SW.,"dgCMatrix"), tolerance =0))
  [1] "Mean relative difference: 1.704172e-08"
  > stopifnot(is.character(ch), length(ch) == 1)## had length(.)  2  previously
  > IW <- solve(WW)
  > class(I1 <- IW %*% WW)# "dge" or "dgC" (!)
  [1] "dgeMatrix"
  attr(,"package")
  [1] "Matrix"
  > class(I2 <- WW %*% IW)
  [1] "dgeMatrix"
  attr(,"package")
  [1] "Matrix"
  > I <- diag(nr=nrow(WW))
  > stopifnot(all.equal(as(I1,"matrix"), I, check.attributes=FALSE, tolerance = 1e-4),
  +           ## "Mean relative difference: 3.296549e-05"  (or "1.999949" for Matrix_1.0-13 !!!)
  +           all.equal(as(I2,"matrix"), I, check.attributes=FALSE)) #default tol gives "1" for M.._1.0-13
  > 
  > if(doExtras) {
  +     print(kappa(WW)) ## [1] 5.129463e+12
  +     print(rcond(WW)) ## [1] 6.216103e-14
  +     ## Warning message: rcond(.) via sparse -> dense coercion
  + }
  > 
  > class(Iw. <- solve(SW.))# FIXME? should be "symmetric" but is not
  [1] "dgCMatrix"
  attr(,"package")
  [1] "Matrix"
  > class(Iw2 <- solve(SW2))# FIXME? should be "symmetric" but is not
  [1] "dgCMatrix"
  attr(,"package")
  [1] "Matrix"
  > class(IW. <- as(Iw., "denseMatrix"))
  [1] "dgeMatrix"
  attr(,"package")
  [1] "Matrix"
  > class(IW2 <- as(Iw2, "denseMatrix"))
  [1] "dgeMatrix"
  attr(,"package")
  [1] "Matrix"
  > 
  > ### The next two were wrong for very long, too
  > assert.EQ.(I, as.matrix(IW. %*% SW.), tol= 4e-4)
  Error in assert.EQ(target, current, tol = tol, giveRE = giveRE, showOnly = showOnly,  : 
    all.equal() |-> Mean relative difference: 0.00576082
  Calls: assert.EQ. -> assert.EQ
  Execution halted
Running the tests in 'tests/validObj.R' failed.
Complete output:
  > library(Matrix)
  > ### Do all kinds of object creation and coercion
  > 
  > source(system.file("test-tools.R", package = "Matrix"))
  Loading required package: tools
  > 
  > ## the empty ones:
  > checkMatrix(new("dgeMatrix"))
   norm(m [0 x 0]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: identical
  m >= m for all: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > checkMatrix(Matrix(,0,0))
   norm(m [0 x 0]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > 
  > ## "dge"
  > assertError( new("dgeMatrix", Dim = c(2,2), x= 1:4) )# double 'Dim'
  > assertError( new("dgeMatrix", Dim = as.integer(c(2,2)), x= 1:4) )# int 'x'
  > assertError( new("dgeMatrix", Dim = 2:2, x=as.double(1:4)) )# length(Dim) !=2
  > assertError( new("dgeMatrix", Dim = as.integer(c(2,2)), x= as.double(1:5)))
  > 
  > checkMatrix(m1 <- Matrix(1:6, ncol=2))
   norm(m [3 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  > checkMatrix(m2 <- Matrix(1:7 +0, ncol=3)) # a (desired) warning
   norm(m [3 x 3]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  Warning message:
  In Matrix(1:7 + 0, ncol = 3) :
    data length [7] is not a sub-multiple or multiple of the number of rows [3]
  > c("dgeMatrix", "ddenseMatrix", "generalMatrix", "geMatrix", "dMatrix",
  +   "denseMatrix", "compMatrix", "Matrix", "xMatrix", "mMatrix") -> m1.cl
  > stopifnot(!anyNA(match(m1.cl, is(m1))),
  + 	  dim(t(m1)) == 2:3, identical(m1, t(t(m1))))
  > c.nam <- paste("C",1:2, sep='')
  > dimnames(m1) <- list(NULL, c.nam)
  > checkMatrix(m1) # failed in 0.999375-10
   norm(m [3 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  > checkMatrix(tm1 <- t(m1))
   norm(m [2 x 3]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  > stopifnot(colnames(m1) == c.nam,
  + 	  identical(dimnames(tm1), list(c.nam, NULL)),
  + 	  identical(m1, t(tm1)))
  > 
  > ## an example of *named* dimnames
  > (t34N <- as(unclass(table(x = gl(3,4), y=gl(4,3))), "dgeMatrix"))
  3 x 4 Matrix of class "dgeMatrix"
     y
  x   1 2 3 4
    1 3 1 0 0
    2 0 2 2 0
    3 0 0 1 3
  > stopifnot(identical(dimnames(t34N),
  + 		    dimnames(as(t34N, "matrix"))),
  +           identical(t34N, t(t(t34N))))
  > 
  > ## "dpo"
  > checkMatrix(cm <- crossprod(m1))
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > checkMatrix(cp <- as(cm, "dppMatrix"))# 'dpp' + factors
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > checkMatrix(cs <- as(cm, "dsyMatrix"))# 'dsy' + factors
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > checkMatrix(dcm <- as(cm, "dgeMatrix"))#'dge'
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: identical
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > checkMatrix(mcm <- as(cm, "dMatrix")) # 'dsy' + factors -- buglet? rather == cm?
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > checkMatrix(mc. <- as(cm, "Matrix"))  # dpo --> dsy -- (as above)  FIXME? ??
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > stopifnot(identical(mc., mcm),
  + 	  identical(cm, (2*cm)/2),# remains dpo
  + 	  identical(cm + cp, cp + cs),# dge
  + 	  identical(mc., mcm),
  + 	  all(2*cm == mcm * 2))
  > 
  > checkMatrix(eq <- cm == cs)
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  > stopifnot(all(eq@x),
  + 	  identical3(pack(eq), cs == cp, cm == cp),
  + 	  as.logical(!(cs < cp)),
  + 	  identical4(!(cs < cp), !(cp > cs), cp <= cs, cs >= cp))
  > 
  > ## Coercion to 'dpo' should give an error if result would be invalid
  > M <- Matrix(diag(4) - 1)
  > assertError(as(M, "dpoMatrix"))
  > M. <- as(M, "dgeMatrix")
  > M.[1,2] <- 10 # -> not even symmetric anymore
  > assertError(as(M., "dpoMatrix"))
  > 
  > 
  > ## Cholesky
  > checkMatrix(ch <- chol(cm))
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  as(<triangular (ge)matrix>, Cholesky): valid: TRUE 
  > checkMatrix(ch2 <- chol(as(cm, "dsyMatrix")))
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  as(<triangular (ge)matrix>, Cholesky): valid: TRUE 
  > checkMatrix(ch3 <- chol(as(cm, "dgeMatrix")))
   norm(m [2 x 2]) : 1  I  F  M  ok
   Summary:  ok
  2*m =?= m+m: ok
  m >= m for all: ok
  m < m for none: ok
  symmpart(m) + skewpart(m) == m: ok;  determinant(): ok
  as(<triangular (ge)matrix>, Cholesky): valid: TRUE 
  > stopifnot(is.all.equal3(as(ch, "matrix"), as(ch2, "matrix"), as(ch3, "matrix")))
  > ### Very basic	triangular matrix stuff
  > 
  > assertError( new("dtrMatrix", Dim = c(2,2), x= 1:4) )# double 'Dim'
  > assertError( new("dtrMatrix", Dim = as.integer(c(2,2)), x= 1:4) )# int 'x'
  > ## This caused a segfault (before revision r1172 in ../src/dtrMatrix.c):
  > assertError( new("dtrMatrix", Dim = 2:2, x=as.double(1:4)) )# length(Dim) !=2
  > assertError( new("dtrMatrix", Dim = as.integer(c(2,2)), x= as.double(1:5)))
  > 
  > tr22 <- new("dtrMatrix", Dim = as.integer(c(2,2)), x=as.double(1:4))
  > tt22 <- t(tr22)
  > (tPt <- tr22 + tt22)
  2 x 2 Matrix of class "dgeMatrix"
       [,1] [,2]
  [1,]    2    3
  [2,]    3    8
  > stopifnot(identical(10 * tPt, tPt * 10),
  + 	  as.vector(t.22 <- (tr22 / .5)* .5) == c(1,0,3,4),
  + 	  TRUE) ## not yet: class(t.22) == "dtrMatrix")
  > 
  > ## non-square triagonal Matrices --- are forbidden ---
  > assertError(new("dtrMatrix", Dim = 2:3,
  + 		x=as.double(1:6), uplo="L", diag="U"))
  > 
  > n <- 3:3
  > assertError(new("dtCMatrix", Dim = c(n,n), diag = "U"))
  > validObject(T <- new("dtTMatrix", Dim = c(n,n), diag = "U"))
  [1] TRUE
  > validObject(M <- new("dtCMatrix", Dim = c(n,n), diag = "U",
  + 		     p = rep.int(0:0, n+1)))
  [1] TRUE
  > stopifnot(identical(as.mat(T), diag(n)))
  > 
  > set.seed(3) ; (p9 <- as(sample(9), "pMatrix"))
  9 x 9 sparse Matrix of class "pMatrix"
                         
   [1,] . . . . | . . . .
   [2,] . | . . . . . . .
   [3,] . . . | . . . . .
   [4,] . . . . . . | . .
   [5,] . . . . . . . | .
   [6,] . . . . . | . . .
   [7,] . . | . . . . . .
   [8,] | . . . . . . . .
   [9,] . . . . . . . . |
  > ## Check that the correct error message is triggered:
  > ind.try <- try(p9[1,1] <- 1, silent = TRUE)
  > np9 <- as(p9, "ngTMatrix")
  > stopifnot(grepl("replacing.*sensible", ind.try[1]),
  + 	  is.logical(p9[1,]),
  + 	  is(p9[2,, drop=FALSE], "indMatrix"),
  + 	  is(p9[9:1,], "indMatrix"),
  + 	  isTRUE(p9[-c(1:6, 8:9), 1]),
  + 	  identical(t(p9), solve(p9)),
  + 	  identical(p9[TRUE, ], p9),
  +           all.equal(p9[, TRUE], np9), # currently...
  + 	  identical(as(diag(9), "pMatrix"), as(1:9, "pMatrix"))
  + 	  )
  Error: isTRUE(p9[-c(1:6, 8:9), 1]) is not TRUE
  Execution halted