SCM

SCM Repository

[matrix] Annotation of /pkg/tests/simple.R
ViewVC logotype

Annotation of /pkg/tests/simple.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 961 - (view) (download)

1 : maechler 948 #### Currently a collection of simple tests
2 :     ## (since 'Matrix' takes long to load, rather have fewer source files!)
3 :    
4 : bates 346 library(Matrix)
5 : maechler 509
6 : maechler 948 source(system.file("test-tools.R", package = "Matrix"))# identical3() etc
7 :    
8 :     ###-- Sparse Triangular :
9 :    
10 :     (t1 <- new("dtTMatrix", x= c(3,7), i= 0:1, j=3:2,
11 :     Dim= as.integer(c(4,4))))
12 :     stopifnot(validObject(t1),
13 :     validObject(t1c <- as(t1, "dtCMatrix")))
14 :     assert.EQ.mat(t1, as(t1c, "matrix"))
15 :    
16 :     ## from 0-diagonal to unit-diagonal {low-level step}:
17 :     tu <- t1 ; tu@diag <- "U"
18 :     tu
19 :     stopifnot(validObject(cu <- as(tu, "dtCMatrix")),
20 :     validObject(t(cu)),
21 :     validObject(t(tu)))
22 :     assert.EQ.mat(cu, as(tu,"matrix"), tol=0)
23 :    
24 :    
25 :     ###-- Numeric Dense: Crossprod & Solve
26 :    
27 : maechler 509 set.seed(123)
28 : bates 346 mm <- Matrix(rnorm(500 * 150), nc = 150)
29 : maechler 509 stopifnot(validObject(mm))
30 :     xpx <- crossprod(mm)# alters mm !
31 : maechler 873 stopifnot(validObject(mm),
32 :     validObject(xpx))
33 :     str(mm) # 'dge*"
34 :     str(xpx)# 'dpo*"
35 : bates 346 xpy <- crossprod(mm, rnorm(500))
36 :     res <- solve(xpx, xpy)
37 : maechler 873 str(xpx)# now with Cholesky factor
38 :     stopifnot(validObject(xpx),
39 :     validObject(xpy),
40 :     validObject(res))
41 :     stopifnot(all.equal(xpx %*% res, xpy, tol= 1e-12))
42 : maechler 509
43 : maechler 948 ###-- more solve() methods {was ./solve.R }
44 :    
45 :     ## first for "dgeMatrix" and all kinds of RHS :
46 :     (m6 <- 1 + as(diag(0:5), "dgeMatrix"))
47 :     rcond(m6)
48 :     I6 <- as(diag(6), "dgeMatrix")
49 :     stopifnot(all.equal(I6, m6 %*% solve(m6)),
50 :     all.equal(I6, solve(m6) %*% m6) )
51 :    
52 :     (i6 <- solve(m6, Matrix(1:6)))
53 :     stopifnot(identical(i6, as(cbind(c(-4, rep(1,5))), "dgeMatrix")),
54 : maechler 949 identical(i6, solve(m6, 1:6)),
55 : maechler 948 identical(i6, solve(m6, matrix(1:6))),
56 : maechler 949 identical(i6, solve(m6, matrix(c(1,2,3,4,5,6))))
57 : maechler 948 )
58 :    
59 :     ###-- row- and column operations {was ./rowcolOps.R }
60 :    
61 :     set.seed(321)
62 : maechler 956 m1k <- Matrix(round(rnorm(1000), 2), 50, 20)
63 :     m.m <- as(m1k, "matrix")
64 :     stopifnot(all.equal(colMeans(m1k), colMeans(m.m)),
65 :     all.equal(colSums (m1k), colSums (m.m)),
66 :     all.equal(rowMeans(m1k), rowMeans(m.m)),
67 :     all.equal(rowSums (m1k), rowSums (m.m)))
68 : maechler 948
69 :     ###-- Testing expansions of factorizations {was ./expand.R }
70 :    
71 :     (m1 <- round(Matrix(rnorm(25), 5), 2))
72 :     (lul <- expand(lu(m1)))
73 :     stopifnot(all.equal(as(m1, "matrix"),
74 :     as(lul$P %*% (lul$L %*% lul$U), "matrix")))
75 :    
76 :    
77 : maechler 956 ###-- kronecker for nonsparse uses Matrix(.):
78 :     stopifnot(is(kr <- kronecker(m1, m6), "Matrix"))
79 :     assert.EQ.mat(kr,
80 :     kronecker(as(m1, "matrix"),
81 :     as(m6, "matrix")),
82 :     tol = 0)
83 :     ## sparse:
84 :     (kt1 <- kronecker(t1, tu))
85 :     kt2 <- kronecker(t1c, cu)
86 :     ktf <- kronecker(as.matrix(t1), as.matrix(tu))
87 :    
88 :     assert.EQ.mat(kt1, ktf, tol= 0)
89 :     assert.EQ.mat(kt2, ktf, tol= 0)
90 :     ## but kt1 and kt2, both "dgT" are different since entries are not ordered!
91 :    
92 :    
93 :     ###-- logical sparse : ----------
94 :    
95 :     (lkt <- as(as(kt1, "dgCMatrix"), "lgCMatrix"))# ok
96 :     (clt <- crossprod(lkt))
97 : maechler 961 if(FALSE)
98 : maechler 956 crossprod(clt)
99 :     ## CHOLMOD error: matrix cannot be symmetric
100 :    
101 : maechler 961
102 : maechler 956 ### "d" <-> "l" for (symmetric) sparse :
103 :     data(mm)
104 :     xpx <- crossprod(mm)
105 :     lxpx <- as(xpx, "lsCMatrix")
106 :     if(FALSE)
107 :     show(lxpx) ## gives error about "lsC" -> "lgT" coercion ..
108 :     ## The bug is actually from *subsetting* the large matrix:
109 :     if(FALSE) ## FIXME
110 :     r <- lxpx[1:2,]
111 :    
112 :     lmm <- as(mm, "lgCMatrix")
113 :     xlx <- crossprod(lmm)
114 : maechler 961 ## now A = lxpx and B = xlx should be close, but not quite the same
115 :     ## since <x,y> = 0 is well possible when x!=0 and y!=0 .
116 :     ## However, A[i,j] != 0 ==> B[i,j] != 0:
117 :     A <- as(as(lxpx, "lgCMatrix"), "lgTMatrix")
118 :     B <- as(as(xlx, "lgCMatrix"), "lgTMatrix")
119 :     ij <- function(a) a@i + ncol(a) * a@j
120 :     stopifnot(all(ij(A) %in% ij(B)))
121 : maechler 956
122 : maechler 509 proc.time()

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge