R Under development (unstable) (2023-11-05 r85474 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## for R_DEFAULT_PACKAGES=NULL : > library(stats) > library(utils) > > library(Matrix) > source(system.file("test-tools.R", package = "Matrix")) Loading required package: tools > > data(KNex, package = "Matrix") > mm <- KNex$mm > stopifnot(##is(mm) == c("dgCMatrix", "dMatrix", "Matrix"), + dim(mm) == (dm <- c(1850, 712)), + identical(dimnames(mm), list(NULL,NULL))) > str(mm) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:8755] 0 2 25 27 163 165 1258 1261 1276 1278 ... ..@ p : int [1:713] 0 13 17 26 38 43 52 56 61 67 ... ..@ Dim : int [1:2] 1850 712 ..@ Dimnames:List of 2 .. ..$ : NULL .. ..$ : NULL ..@ x : num [1:8755] 0.277 0.277 0.277 0.277 0.277 ... ..@ factors : list() > tmm <- t(mm) > str(tmm) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:8755] 0 257 427 549 697 1 257 427 697 0 ... ..@ p : int [1:1851] 0 5 9 14 18 22 27 31 36 41 ... ..@ Dim : int [1:2] 712 1850 ..@ Dimnames:List of 2 .. ..$ : NULL .. ..$ : NULL ..@ x : num [1:8755] 0.2774 0.5 -0.574 -0.2995 0.0298 ... ..@ factors : list() > > str(mTm <- crossprod(mm)) Formal class 'dsCMatrix' [package "Matrix"] with 7 slots ..@ i : int [1:4918] 0 1 2 3 4 5 6 7 8 9 ... ..@ p : int [1:713] 0 1 2 3 4 5 6 7 8 9 ... ..@ Dim : int [1:2] 712 712 ..@ Dimnames:List of 2 .. ..$ : NULL .. ..$ : NULL ..@ x : num [1:4918] 1 1 1 1 1 ... ..@ uplo : chr "U" ..@ factors : list() > mmT <- crossprod(tmm) > mmT. <- tcrossprod(mm) > stopifnot(all.equal(mmT, mmT.)) > ## Previously these were not the same > ## Should be the same but not quite: even length( * @ x ) differs! > ##str(mmT, max=2)# much larger than mTm (i.e less sparse) > ##str(mmT., max=2)# x slot is currently slightly larger --> improve tcrossprod()? > ##system.time(ae <- all.equal(as(mmT.,"matrix"), as(mmT,"matrix"), tolerance = 1e-14)) > ## 4-5 seconds on a 850 MHz, P III > ##stopifnot(ae) > > stopifnot(validObject(tmm), dim(tmm) == dm[2:1], + validObject(mTm), dim(mTm) == dm[c(2,2)], + validObject(mmT), dim(mmT) == dm[c(1,1)], + identical(as(tmm, "matrix"), t(as(mm, "matrix")))) > > ## from a bug report by Guissepe Ragusa > RNGversion("3.6.0")# future proof > set.seed(101) > for(i in 1:10) { + A <- matrix(rnorm(400), nrow = 100, ncol = 4) + A[A < +1] <- 0 ; Am <- A + Acsc <- as(Am, "CsparseMatrix") + A <- as(Am, "unpackedMatrix") + b <- matrix(rnorm(400), nrow = 4, ncol = 100) + B <- as(b, "unpackedMatrix") + assert.EQ.mat(A %*% B, Am %*% b) + assert.EQ.mat(B %*% A, b %*% Am) + stopifnot(identical(A, as(Acsc, "unpackedMatrix")), + identical(Acsc, as(A, "CsparseMatrix")), + is.all.equal4(A %*% B, Acsc %*% B, + A %*% b, Acsc %*% b), + is.all.equal4(b %*% A, b %*% Acsc, + B %*% A, B %*% Acsc)) + } > > ###--- dgTMatrix {was ./dgTMatrix.R } ------- > > ### Use ``non-unique'' versions of dgTMatrix objects > > N <- 200 > set.seed(1) > i <- as.integer(round(runif (N, 0, 100))) > j <- as.integer(3* rpois (N, lambda=15)) > x <- round(rnorm(N), 2) > which(duplicated(cbind(i,j))) # 8 index pairs are duplicated [1] 28 43 75 103 139 160 165 193 > > m1 <- new("dgTMatrix", Dim = c(max(i)+1:1, max(j)+1:1), i = i, j = j, x = x) > mc <- as(m1, "CsparseMatrix") > m2 <- as(mc, "TsparseMatrix")## the same as 'm1' but without duplicates > > stopifnot(!isTRUE(all.equal.default(m1, m2)), + all.equal(as(m1,"matrix"), as(m2,"matrix"), tolerance =1e-15), + all.equal(crossprod(m1), crossprod(m2), tolerance =1e-15), + identical(mc, as(m2, "CsparseMatrix"))) > > ### -> uniq* functions now in ../R/Auxiliaries.R > (t2 <- system.time(um2 <- asUniqueT(m1))) user system elapsed 0 0 0 > stopifnot(identical(m2,um2)) > > ### -> error/warning condition for solve() of a singular matrix (Barry Rowlingson) > (M <- Matrix(0+ 1:16, ncol = 4)) 4 x 4 Matrix of class "dgeMatrix" [,1] [,2] [,3] [,4] [1,] 1 5 9 13 [2,] 2 6 10 14 [3,] 3 7 11 15 [4,] 4 8 12 16 > assertError(solve(M), verbose=TRUE)## ".. computationally singular" + warning + caches LU Asserted error: 'a' is computationally singular, rcond(a)=0 > assertError(solve(t(M))) > options(warn=2) # no more warnings allowed from here > lum <- lu(M, warnSing=FALSE) > stopifnot(is(fLU <- M@factors[["denseLU"]], "MatrixFactorization"), + identical(lum, fLU)) > (e.lu <- expand(fLU)) $L 4 x 4 Matrix of class "dtrMatrix" (unitriangular) [,1] [,2] [,3] [,4] [1,] 1.0000000 . . . [2,] 0.2500000 1.0000000 . . [3,] 0.7500000 0.3333333 1.0000000 . [4,] 0.5000000 0.6666667 0.0000000 1.0000000 $U 4 x 4 Matrix of class "dtrMatrix" [,1] [,2] [,3] [,4] [1,] 4 8 12 16 [2,] . 3 6 9 [3,] . . 0 0 [4,] . . . 0 $P 4 x 4 sparse Matrix of class "pMatrix" [1,] . | . . [2,] . . . | [3,] . . | . [4,] | . . . > M2 <- with(e.lu, P %*% L %*% U) > assert.EQ.mat(M2, as(M, "matrix")) > ## now the sparse LU : > M. <- as(M,"sparseMatrix") > tt <- try(solve(M.)) # less nice: factor is *not* cached Error in .local(x, ...) : LU factorization of .gCMatrix failed: out of memory or near-singular > ## use a non-singular one: > M1 <- M. + 0.5*Diagonal(nrow(M.)) > luM1 <- lu(M1) > d1 <- determinant(as(M1,"denseMatrix")) > stopifnot(identical(luM1, M1@factors[["sparseLU~"]]), + diag(luM1@L) == 1,# L is *unit*-triangular + all.equal(log(-prod(diag(luM1@U))), c(d1$modulus))) > > cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons'' Time elapsed: 1.34 0.15 1.48 NA NA > > proc.time() user system elapsed 1.34 0.15 1.48