File: dtpMatrix.R

package info (click to toggle)
rmatrix 0.95.5-1
  • links: PTS
  • area: main
  • in suites: sarge
  • size: 4,732 kB
  • ctags: 2,028
  • sloc: ansic: 22,357; makefile: 74; sh: 28
file content (55 lines) | stat: -rw-r--r-- 1,741 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
### triangular packed
library(Matrix)

cp6 <- chol(Hilbert(6))
tp6 <- as(cp6,"dtpMatrix")
round(tp6, 3)## round() is "Math2" group method
1/tp6        ## "Arith" group : gives 'dgeMatrix'
str(tp6)
stopifnot(validObject(tp6),
          all.equal(tp6 %*% diag(6), as(tp6, "dgeMatrix")),
          validObject(tp6. <- diag(6) %*% tp6),
          class((tt6 <- t(tp6))) == "dtpMatrix",
          identical(t(tt6), tp6),
          tp6@uplo == "U" && tt6@uplo == "L")

all.equal(as(tp6.,"matrix"),
          as(tp6, "matrix"), tol= 1e-15)
(tr6 <- as(tp6, "dtrMatrix")) ## prints using wrong class name
D. <- determinant(tp6)
rc <- rcond(tp6)
stopifnot(all.equal(D.$modulus, -6.579251212),
          all.equal(rc, 1.791511257e-4),
          rc == tp6@rcond,
          all.equal(norm(tp6, "I") , 2.45),
          all.equal(norm(tp6, "1") , 1),
          all.equal(norm(tp6, "F") , 1.37047826623)
          )
object.size(tp6)
object.size(as(tp6, "dtrMatrix"))
object.size(as(tp6, "matrix"))
D6 <- as(diag(6), "dgeMatrix")
ge6 <- as(tp6, "dgeMatrix")
## Direct all.equal() fails, because ge6 has 'rcond', but the product not:
mge6 <- as(ge6, "matrix")
stopifnot(all.equal(as(D6 %*% tp6,"matrix"), mge6),
          all.equal(as(tp6 %*% D6,"matrix"), mge6)
          )

## larger case
set.seed(123)
rl <- new("dtpMatrix", uplo="L", diag="N", Dim = rep.int(1000:1000,2),
          x = rnorm(500*1001))
validObject(rl)
str(rl)
sapply(c("I", "1", "F"), function(type) norm(rl, type=type))
rcond(rl)# 0 !
stopifnot(all.equal(as(rl %*% diag(1000),"matrix"),
                    as(rl, "matrix")))
object.size(rl) ## 4 mio
object.size(as(rl, "dtrMatrix"))# 8 mio
object.size(as(rl, "matrix"))# dito
determinant(rl)


proc.time() # for ``statistical reasons''