R : Copyright 2006, The R Foundation for Statistical Computing
Version 2.3.0 Under development (unstable) (2006-01-02 r36956)
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ### Testing positive definite matrices
> 
> library(Matrix)
> 
> 
> h9 <- Hilbert(9)
> stopifnot(c(0,0) == dim(Hilbert(0)),
+           c(9,9) == dim(h9))
> str(h9)
Formal class 'dpoMatrix' [package "Matrix"] with 6 slots
  ..@ rcond   : num(0) 
  ..@ x       : num [1:81] 1.000 0.500 0.333 0.250 0.200 ...
  ..@ Dim     : int [1:2] 9 9
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ factors : list()
  ..@ uplo    : chr "U"
> all.equal(c(determinant(h9)$modulus), -96.7369456, tol= 2e-8)
[1] TRUE
> stopifnot(0 == length(h9@factors))# nothing yet
> round(ch9 <- chol(h9), 3) ## round() preserves 'triangular' !
9 x 9 Matrix of class "Cholesky"
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] 
 [1,] 1.000 0.500 0.333 0.250 0.200 0.167 0.143 0.125 0.111
 [2,]     . 0.289 0.289 0.260 0.231 0.206 0.186 0.168 0.154
 [3,]     .     . 0.075 0.112 0.128 0.133 0.133 0.130 0.126
 [4,]     .     .     . 0.019 0.038 0.052 0.063 0.070 0.075
 [5,]     .     .     .     . 0.005 0.012 0.019 0.027 0.033
 [6,]     .     .     .     .     . 0.001 0.004 0.007 0.010
 [7,]     .     .     .     .     .     . 0.000 0.001 0.002
 [8,]     .     .     .     .     .     .     . 0.000 0.000
 [9,]     .     .     .     .     .     .     .     . 0.000
> str(f9 <- as(chol(h9), "dtrMatrix"))
Formal class 'dtrMatrix' [package "Matrix"] with 7 slots
  ..@ rcond   : num(0) 
  ..@ x       : num [1:81] 1 0 0 0 0 0 0 0 0 0.5 ...
  ..@ Dim     : int [1:2] 9 9
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ factors : list()
  ..@ uplo    : chr "U"
  ..@ diag    : chr "N"
> ## h9 now has factorization
> stopifnot(names(h9@factors) == "Cholesky",
+           all.equal(rcond(h9), 9.0938e-13),
+           all.equal(rcond(f9), 9.1272e-7, tol = 1e-6))# more precision fails
> str(h9)# has 'rcond' and 'factors'
Formal class 'dpoMatrix' [package "Matrix"] with 6 slots
  ..@ rcond   : Named num 9.1e-13
  .. ..- attr(*, "names")= chr "O"
  ..@ x       : num [1:81] 1.000 0.500 0.333 0.250 0.200 ...
  ..@ Dim     : int [1:2] 9 9
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ factors :List of 1
  .. ..$ Cholesky:Formal class 'Cholesky' [package "Matrix"] with 7 slots
  .. .. .. ..@ rcond   : num(0) 
  .. .. .. ..@ x       : num [1:81] 1 0 0 0 0 0 0 0 0 0.5 ...
  .. .. .. ..@ Dim     : int [1:2] 9 9
  .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : NULL
  .. .. .. ..@ factors : list()
  .. .. .. ..@ uplo    : chr "U"
  .. .. .. ..@ diag    : chr "N"
  ..@ uplo    : chr "U"
> options(digits=4)
> (cf9 <- crossprod(f9))# looks the same as  h9 :
9 x 9 Matrix of class "dpoMatrix"
        [,1]   [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
 [1,] 1.0000 0.5000 0.33333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111
 [2,] 0.5000 0.3333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000
 [3,] 0.3333 0.2500 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091
 [4,] 0.2500 0.2000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333
 [5,] 0.2000 0.1667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692
 [6,] 0.1667 0.1429 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143
 [7,] 0.1429 0.1250 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667
 [8,] 0.1250 0.1111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250
 [9,] 0.1111 0.1000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250 0.05882
> stopifnot(all.equal(as.matrix(h9),
+                     as.matrix(cf9), tol= 1e-15))
> 
> str(hp9 <- as(h9, "dppMatrix"))
Formal class 'dppMatrix' [package "Matrix"] with 6 slots
  ..@ rcond   : Named num 9.1e-13
  .. ..- attr(*, "names")= chr "O"
  ..@ x       : num [1:45] 1.000 0.500 0.333 0.333 0.250 ...
  ..@ Dim     : int [1:2] 9 9
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ factors :List of 1
  .. ..$ pCholesky:Formal class 'pCholesky' [package "Matrix"] with 7 slots
  .. .. .. ..@ rcond   : num(0) 
  .. .. .. ..@ x       : num [1:45] 1.000 0.500 0.289 0.333 0.289 ...
  .. .. .. ..@ Dim     : int [1:2] 9 9
  .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : NULL
  .. .. .. ..@ factors : list()
  .. .. .. ..@ uplo    : chr "U"
  .. .. .. ..@ diag    : chr "N"
  ..@ uplo    : chr "U"
> 
> s9 <- solve(hp9, seq(nrow(hp9)))
> signif(t(s9)/10000, 4)# only rounded numbers are platform-independent
1 x 9 Matrix of class "dgeMatrix"
       [,1]  [,2]  [,3]   [,4] [,5]  [,6] [,7]  [,8] [,9]
[1,] 0.0729 -5.76 109.5 -864.9 3468 -7668 9459 -6095 1597
> (I9 <- hp9 %*% s9)
9 x 1 Matrix of class "dgeMatrix"
      [,1]
 [1,]    1
 [2,]    2
 [3,]    3
 [4,]    4
 [5,]    5
 [6,]    6
 [7,]    7
 [8,]    8
 [9,]    9
> m9 <- matrix(1:9, dimnames = list(NULL,NULL))
> stopifnot(all.equal(m9, as.matrix(I9), tol = 2e-9))
> 
> cat('Time elapsed: ', proc.time(),'\n') # "stats"
Time elapsed:  8.44 0.06 8.5 0 0 
> 
>