#### Currently a collection of simple tests ## (since 'Matrix' takes long to load, rather have fewer source files!) library(Matrix) source(system.file("test-tools.R", package = "Matrix"))# identical3() etc ###-- Sparse Triangular : (t1 <- new("dtTMatrix", x= c(3,7), i= 0:1, j=3:2, Dim= as.integer(c(4,4)))) stopifnot(validObject(t1), validObject(t1c <- as(t1, "dtCMatrix"))) assert.EQ.mat(t1, as(t1c, "matrix")) ## from 0-diagonal to unit-diagonal {low-level step}: tu <- t1 ; tu@diag <- "U" tu stopifnot(validObject(cu <- as(tu, "dtCMatrix")), validObject(t(cu)), validObject(t(tu))) assert.EQ.mat(cu, as(tu,"matrix"), tol=0) ###-- Numeric Dense: Crossprod & Solve set.seed(123) mm <- Matrix(rnorm(500 * 150), nc = 150) stopifnot(validObject(mm)) xpx <- crossprod(mm)# alters mm ! stopifnot(validObject(mm), validObject(xpx)) str(mm) # 'dge*" str(xpx)# 'dpo*" xpy <- crossprod(mm, rnorm(500)) res <- solve(xpx, xpy) str(xpx)# now with Cholesky factor stopifnot(validObject(xpx), validObject(xpy), validObject(res)) stopifnot(all.equal(xpx %*% res, xpy, tol= 1e-12)) ###-- more solve() methods {was ./solve.R } ## first for "dgeMatrix" and all kinds of RHS : (m6 <- 1 + as(diag(0:5), "dgeMatrix")) rcond(m6) I6 <- as(diag(6), "dgeMatrix") stopifnot(all.equal(I6, m6 %*% solve(m6)), all.equal(I6, solve(m6) %*% m6) ) (i6 <- solve(m6, Matrix(1:6))) stopifnot(identical(i6, as(cbind(c(-4, rep(1,5))), "dgeMatrix")), identical(i6, solve(m6, 1:6)), identical(i6, solve(m6, matrix(1:6))), identical(i6, solve(m6, matrix(c(1,2,3,4,5,6)))) ) ###-- row- and column operations {was ./rowcolOps.R } set.seed(321) m1k <- Matrix(round(rnorm(1000), 2), 50, 20) m.m <- as(m1k, "matrix") stopifnot(all.equal(colMeans(m1k), colMeans(m.m)), all.equal(colSums (m1k), colSums (m.m)), all.equal(rowMeans(m1k), rowMeans(m.m)), all.equal(rowSums (m1k), rowSums (m.m))) ###-- Testing expansions of factorizations {was ./expand.R } (m1 <- round(Matrix(rnorm(25), 5), 2)) (lul <- expand(lu(m1))) stopifnot(all.equal(as(m1, "matrix"), as(lul\$P %*% (lul\$L %*% lul\$U), "matrix"))) ###-- kronecker for nonsparse uses Matrix(.): stopifnot(is(kr <- kronecker(m1, m6), "Matrix")) assert.EQ.mat(kr, kronecker(as(m1, "matrix"), as(m6, "matrix")), tol = 0) ## sparse: (kt1 <- kronecker(t1, tu)) kt2 <- kronecker(t1c, cu) ktf <- kronecker(as.matrix(t1), as.matrix(tu)) assert.EQ.mat(kt1, ktf, tol= 0) assert.EQ.mat(kt2, ktf, tol= 0) ## but kt1 and kt2, both "dgT" are different since entries are not ordered! ###-- logical sparse : ---------- (lkt <- as(as(kt1, "dgCMatrix"), "lgCMatrix"))# ok (clt <- crossprod(lkt)) if(FALSE) crossprod(clt) ## CHOLMOD error: matrix cannot be symmetric ### "d" <-> "l" for (symmetric) sparse : data(mm) xpx <- crossprod(mm) lxpx <- as(xpx, "lsCMatrix") if(FALSE) show(lxpx) ## gives error about "lsC" -> "lgT" coercion .. ## The bug is actually from *subsetting* the large matrix: if(FALSE) ## FIXME r <- lxpx[1:2,] lmm <- as(mm, "lgCMatrix") xlx <- crossprod(lmm) ## now A = lxpx and B = xlx should be close, but not quite the same ## since = 0 is well possible when x!=0 and y!=0 . ## However, A[i,j] != 0 ==> B[i,j] != 0: A <- as(as(lxpx, "lgCMatrix"), "lgTMatrix") B <- as(as(xlx, "lgCMatrix"), "lgTMatrix") ij <- function(a) a@i + ncol(a) * a@j stopifnot(all(ij(A) %in% ij(B))) proc.time()