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

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

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

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

> ### triangular packed
> 
> ## for R_DEFAULT_PACKAGES=NULL :
> library(stats)
> library(utils)
> 
> library(Matrix)
> source(system.file("test-tools.R", package = "Matrix"))# identical3() etc
Loading required package: tools
> 
> cp6 <- as(Cholesky(H6 <- Hilbert(6), perm = FALSE), "dtrMatrix")
> (tp6 <- as(cp6, "packedMatrix"))
6 x 6 Matrix of class "dtpMatrix"
     [,1]        [,2]        [,3]        [,4]        [,5]        [,6]       
[1,] 1.000000000 0.500000000 0.333333333 0.250000000 0.200000000 0.166666667
[2,]           . 0.288675135 0.288675135 0.259807621 0.230940108 0.206196525
[3,]           .           . 0.074535599 0.111803399 0.127775313 0.133099284
[4,]           .           .           . 0.018898224 0.037796447 0.052495066
[5,]           .           .           .           . 0.004761905 0.011904762
[6,]           .           .           .           .           . 0.001196474
> round(tp6, 3)## round() is "Math2" group method
6 x 6 Matrix of class "dtpMatrix"
     [,1]  [,2]  [,3]  [,4]  [,5]  [,6] 
[1,] 1.000 0.500 0.333 0.250 0.200 0.167
[2,]     . 0.289 0.289 0.260 0.231 0.206
[3,]     .     . 0.075 0.112 0.128 0.133
[4,]     .     .     . 0.019 0.038 0.052
[5,]     .     .     .     . 0.005 0.012
[6,]     .     .     .     .     . 0.001
> 1/tp6        ## "Arith" group : gives 'dgeMatrix'
6 x 6 Matrix of class "dgeMatrix"
     [,1]     [,2]      [,3]      [,4]       [,5]       [,6]
[1,]    1 2.000000  3.000000  4.000000   5.000000   6.000000
[2,]  Inf 3.464102  3.464102  3.849002   4.330127   4.849742
[3,]  Inf      Inf 13.416408  8.944272   7.826238   7.513188
[4,]  Inf      Inf       Inf 52.915026  26.457513  19.049409
[5,]  Inf      Inf       Inf       Inf 210.000000  84.000000
[6,]  Inf      Inf       Inf       Inf        Inf 835.789447
> str(tp6)
Formal class 'dtpMatrix' [package "Matrix"] with 5 slots
  ..@ uplo    : chr "U"
  ..@ Dim     : int [1:2] 6 6
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:21] 1 0.5 0.289 0.333 0.289 ...
  ..@ diag    : chr "N"
> ## arithmetic with a mix of dimnames / no dimnames
> tp <- tp6; dimnames(tp) <- list(LETTERS[1:6], letters[11:16])
> ## as.matrix() --> "logical" matrix
> stopifnot(as.matrix(tp - tp6 == tp6 - tp),
+ 	  as.matrix(0 == tp - tp6),
+ 	  identical(as(tp6,"CsparseMatrix"),
+ 		    as(cp6,"CsparseMatrix")))
> 
> stopifnot(validObject(tp6),
+           all.equal(tp6 %*% diag(6), as(tp6, "generalMatrix")),
+           validObject(tp6. <- diag(6) %*% tp6),
+           identical(t(tt6 <- t(tp6)), tp6),
+           tp6@uplo == "U" && tt6@uplo == "L")
> 
> all.equal(as(tp6.,"matrix"),
+           as(tp6, "matrix"), tolerance= 1e-15)
[1] TRUE
> dH6 <- determinant(H6)
> D. <- determinant(tp6)
> rc <- rcond(tp6)
> stopifnot(all.equal(dH6$modulus, determinant(as.matrix(H6))$modulus),
+           is.all.equal3(c(D.$modulus), c(dH6$modulus) / 2, -19.883103353),
+           all.equal(rc, 1.791511257e-4),
+           all.equal(norm(tp6, "I") , 2.45),
+           all.equal(norm(tp6, "1") , 1),
+           all.equal(norm(tp6, "F") , 1.37047826623)
+           )
> object.size(tp6)
1632 bytes
> object.size(as(tp6, "unpackedMatrix"))
1752 bytes
> object.size(as(tp6, "matrix"))
504 bytes
> D6 <- as(diag(6), "generalMatrix")
> ge6 <- as(tp6, "generalMatrix")
> stopifnot(all.equal(D6 %*% tp6, ge6),
+           all.equal(tp6 %*% D6, ge6))
> 
> ## larger case
> RNGversion("3.6.0")# future proof
> set.seed(123)
> rl <- new("dtpMatrix", uplo="L", diag="N", Dim = c(1000L, 1000L),
+           x = rnorm(500*1001))
> validObject(rl)
[1] TRUE
> str(rl)
Formal class 'dtpMatrix' [package "Matrix"] with 5 slots
  ..@ uplo    : chr "L"
  ..@ Dim     : int [1:2] 1000 1000
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:500500] -0.5605 -0.2302 1.5587 0.0705 0.1293 ...
  ..@ diag    : chr "N"
> sapply(c("I", "1", "F"), function(type) norm(rl, type=type))
       I        1        F 
822.0263 824.2083 708.1537 
> rcond(rl)# 0 !
[1] 0
> stopifnot(all.equal(as(rl %*% diag(1000),"matrix"),
+                     as(rl, "matrix")))
> object.size(rl) ## 4 MB
4005464 bytes
> object.size(as(rl, "unpackedMatrix"))# 8 MB
8001464 bytes
> object.size(as(rl, "matrix"))# ditto
8000216 bytes
> print(drl <- determinant(rl), digits = 12)
$modulus
[1] -638.257312422
attr(,"logarithm")
[1] TRUE

$sign
[1] 1

attr(,"class")
[1] "det"
> stopifnot(all.equal(c(drl$modulus), -638.257312422))
> 
> 
> cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
Time elapsed:  1.95 0.17 2.1 NA NA 
> 
> proc.time()
   user  system elapsed 
   1.95    0.17    2.10