# SCM Repository

[matrix] Diff of /pkg/tests/factorizing.R
 [matrix] / pkg / tests / factorizing.R

# Diff of /pkg/tests/factorizing.R

revision 2136, Mon Mar 17 22:20:29 2008 UTC revision 2137, Mon Mar 17 22:21:24 2008 UTC
# Line 56  Line 56
56
57  ## FIXME:  expand(pmLU)  ## FIXME:  expand(pmLU)
58
59  ## Cholesky()  ###________ Cholesky() ________
60
61    ##--------  LDL' ---- small exact examples
62
63    set.seed(1)
64    for(n in c(5:12)) {
65        cat("\nn = ",n,"\n-------\n")
66        rr <- mkLDL(n)
67        ##    -------- from 'test-tools.R'
68        stopifnot(all(with(rr, A == as(L %*% D %*% t(L),
69                               "symmetricMatrix"))))
70        d <- rr\$d.half
71        A <- rr\$A
72        R <- chol(A)
73        print(d. <- diag(R))
74        D. <- Diagonal(x= d.^2)
75        L. <- t(R) %*% Diagonal(x = 1/d.)
76        stopifnot(all.equal(as.matrix(D.), as.matrix(rr\$ D)),
77                  all.equal(as.matrix(L.), as.matrix(rr\$ L)))
78        ##
79        CAp <- Cholesky(A)# perm=TRUE --> Permutation:
80        p <- CAp@perm + 1L
81        P <- as(p, "pMatrix")
82        ## the inverse permutation:
83        invP <- solve(P)@perm
84        lDet <- sum(2* log(d))# the "true" value
85        ldet <- Matrix:::.diag.dsC(Chx = CAp, res.kind = "sumLog")
86        ##
87        CA  <- Cholesky(A,perm=FALSE)
88        CA
89        ldet2 <- Matrix:::.diag.dsC(Chx = CA, res.kind = "sumLog")
90        ## not printing CAp : ends up non-integer for n >= 11
91        mCAp <- as(CAp,"sparseMatrix")
92        print(mCA  <- drop0(as(CA, "sparseMatrix")))
93        stopifnot(identical(A[p,p], as(P %*% A %*% t(P),
94                                       "symmetricMatrix")),
95                  all.equal(lDet, sum(log(Matrix:::.diag.dsC(Chx= CAp,res.kind="diag")))),
96                  relErr(d.^2, Matrix:::.diag.dsC(Chx= CA, res.kind="diag")) < 1e-14,
97                  all.equal(lDet, ldet),
98                  all.equal(lDet, ldet2),
99                  relErr(A[p,p], tcrossprod(mCAp)) < 1e-14)
100    }## for()
101
102    set.seed(17)
103    (rr <- mkLDL(4))
104    (CA <- Cholesky(rr\$A))
105    stopifnot(all.equal(determinant(rr\$A),
106                        determinant(as(rr\$A, "matrix"))))
107
108    ## --- now a "large" (712 x 712) real data example
109
110  data(KNex)  data(KNex)
111  mtm <- with(KNex, crossprod(mm))  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")
114    ld.4 <- .Call("dsCMatrix_LDL_D", mtm, perm=FALSE, "sumLog")# clearly slower
115
116  c1 <- Cholesky(mtm)  c1 <- Cholesky(mtm)
117  c2 <- Cholesky(mtm, super = TRUE)  c2 <- Cholesky(mtm, super = TRUE)
118  bv <- 1:nrow(mtm) # even integer  bv <- 1:nrow(mtm) # even integer
# Line 73  Line 127
127      stopifnot(dim(x) == c(712, 1),      stopifnot(dim(x) == c(712, 1),
128                identical(x, solve(c2, bv, system = sys)))                identical(x, solve(c2, bv, system = sys)))
129  }  }
130
131  ## log(|LL'|) - check if super = TRUE and simplicial give same determinant  ## log(|LL'|) - check if super = TRUE and simplicial give same determinant
132  all.equal(.Call("CHMfactor_ldetL2", c1), .Call("CHMfactor_ldetL2", c2))  ld1 <- .Call("CHMfactor_ldetL2", c1)
133  ## FIXME: Add other checks here for determinant(mtm) when available  ld2 <- .Call("CHMfactor_ldetL2", c2)
134    (ld1. <- determinant(mtm))
135    ## experimental
136    ld3 <- .Call("dsCMatrix_LDL_D", mtm, TRUE, "sumLog")
137    ld4 <- .Call("dsCMatrix_LDL_D", mtm, FALSE, "sumLog")
138    stopifnot(all.equal(ld1, ld2),
139              is.all.equal3(ld2, ld3, ld4),
140              all.equal(ld.3, ld3, tol = 1e-14),
141              all.equal(ld.4, ld4, tol = 1e-14),
142              ## must be identical, based on same CHMsimpl object:
143              identical(ld1, as.vector(ld1.\$modulus)))
144
145
146  ## Schur() ----------------------  ## Schur() ----------------------
147  checkSchur <- function(A, SchurA = Schur(A), tol = 1e-14) {  checkSchur <- function(A, SchurA = Schur(A), tol = 1e-14) {

Legend:
 Removed from v.2136 changed lines Added in v.2137