The R Project SVN R-packages

Rev

Rev 4998 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 4998 Rev 5012
Line 107... Line 107...
107 
 
 107 
 
108 
## --- now a "large" (712 x 712) real data example
 108 
## --- now a "large" (712 x 712) real data example
109 
 
 109 
 
110 
data(KNex)
 110 
data(KNex)
111 
mtm <- with(KNex, crossprod(mm))
 111 
mtm <- with(KNex, crossprod(mm))
112 
## BEFORE Cholesky() is called and any factors are cashed inside mtm:
 - 
 
113 
ld.3 <- .Call("dsCMatrix_LDL_D", mtm, perm=TRUE,  "sumLog")
 112 
ld.3 <- .Call("dsCMatrix_LDL_D", mtm, perm=TRUE,  "sumLog")
- 
 
 113 
stopifnot(names(mtm@factors) == "sPDCholesky")
114 
ld.4 <- .Call("dsCMatrix_LDL_D", mtm, perm=FALSE, "sumLog")# clearly slower
 114 
ld.4 <- .Call("dsCMatrix_LDL_D", mtm, perm=FALSE, "sumLog")# clearly slower
115 
stopifnot(length(mtm@factors) == 0)
 115 
stopifnot(names(mtm@factors) == paste(c("sPD", "spD"),"Cholesky", sep=''))
116 
 
 116 
 
- 
 
 117 
## these are now taken from cache
117 
c1 <- Cholesky(mtm)
 118 
c1 <- Cholesky(mtm)
118 
stopifnot(length(mtm@factors) == 1)
 - 
 
119 
c2 <- Cholesky(mtm, super = TRUE)
 119 
c2 <- Cholesky(mtm, super = TRUE)
120 
stopifnot(length(mtm@factors) == 2)
 - 
 
- 
 
 120 
 
121 
bv <- 1:nrow(mtm) # even integer
 121 
bv <- 1:nrow(mtm) # even integer
122 
b <- matrix(bv)
 122 
b <- matrix(bv)
123 
## solve(c2, b) by default solves  Ax = b, where A = c2'c2 !
 123 
## solve(c2, b) by default solves  Ax = b, where A = c2'c2 !
124 
x <- solve(c2,b)
 124 
x <- solve(c2,b)
125 
stopifnot(identical3(x, solve(c2, bv), solve(c2, b, system = "A")),
 125 
stopifnot(identical3(x, solve(c2, bv), solve(c2, b, system = "A")),
Line 145... Line 145...
145 
	  ## must be identical, based on same CHMsimpl object:
 145 
	  ## must be identical, based on same CHMsimpl object:
146 
	  identical(ld1, as.vector(ld1.$modulus)))
 146 
	  identical(ld1, as.vector(ld1.$modulus)))
147 
 
 147 
 
148 
## Some timing measurements
 148 
## Some timing measurements
149 
mtm <- with(KNex, crossprod(mm))
 149 
mtm <- with(KNex, crossprod(mm))
150 
I <- as(as(Diagonal(n=nrow(mtm)), "CsparseMatrix"), "symmetricMatrix")
 150 
I <- .symDiagonal(n=nrow(mtm))
151 
set.seed(101); r <- runif(100)
 151 
set.seed(101); r <- runif(100)
152 
 
 152 
 
153 
system.time(D1 <- sapply(r, function(rho) Matrix:::ldet1.dsC(mtm + (1/rho) * I)))
 153 
system.time(D1 <- sapply(r, function(rho) Matrix:::ldet1.dsC(mtm + (1/rho) * I)))
154 
## 1.099 on fast cmath-5
 154 
## 0.842 on fast cmath-5
155 
system.time(D2 <- sapply(r, function(rho) Matrix:::ldet2.dsC(mtm + (1/rho) * I)))
 155 
system.time(D2 <- sapply(r, function(rho) Matrix:::ldet2.dsC(mtm + (1/rho) * I)))
156 
## 1.167
 156 
## 0.819
157 
system.time(D3 <- sapply(r, function(rho) Matrix:::ldet3.dsC(mtm + (1/rho) * I)))
 157 
system.time(D3 <- sapply(r, function(rho) Matrix:::ldet3.dsC(mtm + (1/rho) * I)))
158 
## 1.007
 158 
## 0.810
159 
stopifnot(is.all.equal3(D1,D2,D3, tol = 1e-13))
 159 
stopifnot(is.all.equal3(D1,D2,D3, tol = 1e-13))
160 
 
 160 
 
161 
## Schur() ----------------------
 161 
## Schur() ----------------------
162 
checkSchur <- function(A, SchurA = Schur(A), tol = 1e-14) {
 162 
checkSchur <- function(A, SchurA = Schur(A), tol = 1e-14) {
163 
    stopifnot(is(SchurA, "Schur"),
 163 
    stopifnot(is(SchurA, "Schur"),