SCM

SCM Repository

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

Diff of /pkg/tests/factorizing.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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