The R Project SVN R

Rev

Rev 17084 | Rev 30396 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 17084 Rev 17161
Line 8... Line 8...
8 
X <- hilbert(9)[,1:6]
 8 
X <- hilbert(9)[,1:6]
9 
str(s <- La.svd(X)); D <- diag(s$d)
 9 
str(s <- La.svd(X)); D <- diag(s$d)
10 
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)#  X = U D V'
 10 
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)#  X = U D V'
11 
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)#  D = U' X V
 11 
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)#  D = U' X V
12 
 
 12 
 
13 
str(s <- La.svd(X, method = "dgesdd")); D <- diag(s$d)
 13 
str(s <- La.svd(X, method = "dgesvd")); D <- diag(s$d)
14 
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)#  X = U D V'
 14 
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)#  X = U D V'
15 
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)#  D = U' X V
 15 
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)#  D = U' X V
16 
 
 16 
 
17 
X <- cbind(1, 1:7)
 17 
X <- cbind(1, 1:7)
18 
str(s <- La.svd(X)); D <- diag(s$d)
 18 
str(s <- La.svd(X)); D <- diag(s$d)
19 
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)#  X = U D V'
 19 
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)#  X = U D V'
20 
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)#  D = U' X V
 20 
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)#  D = U' X V
21 
 
 21 
 
22 
X <- cbind(1, 1:7)
 22 
X <- cbind(1, 1:7)
23 
str(s <- La.svd(X, method = "dgesdd")); D <- diag(s$d)
 23 
str(s <- La.svd(X, method = "dgesvd")); D <- diag(s$d)
24 
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)#  X = U D V'
 24 
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)#  X = U D V'
25 
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)#  D = U' X V
 25 
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)#  D = U' X V
26 
 
 26 
 
27 
# test nu and nv
 27 
# test nu and nv
28 
La.svd(X, nu = 0)
 28 
La.svd(X, nu = 0)
29 
(s <- La.svd(X, nu = 7))
 29 
(s <- La.svd(X, nu = 7))
30 
stopifnot(dim(s$u) == c(7,7))
 30 
stopifnot(dim(s$u) == c(7,7))
31 
La.svd(X, nv = 0)
 31 
La.svd(X, nv = 0)
32 
 
 32 
 
33 
La.svd(X, nu = 0, method = "dgesdd")
 33 
La.svd(X, nu = 0, method = "dgesvd")
34 
(s <- La.svd(X, nu = 7, method = "dgesdd"))
 34 
(s <- La.svd(X, nu = 7, method = "dgesvd"))
35 
stopifnot(dim(s$u) == c(7,7))
 35 
stopifnot(dim(s$u) == c(7,7))
36 
La.svd(X, nv = 0, method = "dgesdd")
 36 
La.svd(X, nv = 0, method = "dgesvd")
37 
 
 37 
 
38 
# test of complex case
 38 
# test of complex case
39 
 
 39 
 
40 
X <- cbind(1, 1:7+(-3:3)*1i)
 40 
X <- cbind(1, 1:7+(-3:3)*1i)
41 
str(s <- La.svd(X)); D <- diag(s$d)
 41 
str(s <- La.svd(X)); D <- diag(s$d)
Line 58... Line 58...
58 
La.eigen(print(cbind(c(0,1i), c(-1i,0))))# Hermite ==> real eigenvalues
 58 
La.eigen(print(cbind(c(0,1i), c(-1i,0))))# Hermite ==> real eigenvalues
59 
## 3 x 3:
 59 
## 3 x 3:
60 
La.eigen(cbind( 1,3:1,1:3))
 60 
La.eigen(cbind( 1,3:1,1:3))
61 
La.eigen(cbind(-1,c(1:2,0),0:2)) # complex values
 61 
La.eigen(cbind(-1,c(1:2,0),0:2)) # complex values
62 
 
 62 
 
63 
La.eigen(cbind(c(1,-1),c(-1,1)), method = "dsyevr")
 63 
La.eigen(cbind(c(1,-1),c(-1,1)), method = "dsyev")
64 
 
 64 
 
65 
 
 65 
 
66 
set.seed(1234)
 66 
set.seed(1234)
67 
Meps <- .Machine$double.eps
 67 
Meps <- .Machine$double.eps
68 
m <- matrix(round(rnorm(25),3), 5,5)
 68 
m <- matrix(round(rnorm(25),3), 5,5)
Line 72... Line 72...
72 
 
 72 
 
73 
stopifnot(
 73 
stopifnot(
74 
 abs(sm %*% V - V %*% diag(lam))          < 60*Meps,
 74 
 abs(sm %*% V - V %*% diag(lam))          < 60*Meps,
75 
 abs(sm       - V %*% diag(lam) %*% t(V)) < 60*Meps)
 75 
 abs(sm       - V %*% diag(lam) %*% t(V)) < 60*Meps)
76 
 
 76 
 
77 
em <- La.eigen(sm, method = "dsyevr"); V <- em$vect
 77 
em <- La.eigen(sm, method = "dsyev"); V <- em$vect
78 
print(lam <- em$values) # ordered DEcreasingly
 78 
print(lam <- em$values) # ordered DEcreasingly
79 
 
 79 
 
80 
stopifnot(
 80 
stopifnot(
81 
 abs(sm %*% V - V %*% diag(lam))          < 60*Meps,
 81 
 abs(sm %*% V - V %*% diag(lam))          < 60*Meps,
82 
 abs(sm       - V %*% diag(lam) %*% t(V)) < 60*Meps)
 82 
 abs(sm       - V %*% diag(lam) %*% t(V)) < 60*Meps)
83 
 
 83 
 
84 
# check only.values = TRUE too
 84 
# check only.values = TRUE too
85 
La.eigen(sm, only.values = TRUE)
 85 
La.eigen(sm, only.values = TRUE)
86 
La.eigen(sm, only.values = TRUE, method = "dsyevr")
 86 
La.eigen(sm, only.values = TRUE, method = "dsyev")
87 
 
 87 
 
88 
 
 88 
 
89 
## symmetric = FALSE
 89 
## symmetric = FALSE
90 
 
 90 
 
91 
em <- La.eigen(sm, symmetric = FALSE); V2 <- em$vect
 91 
em <- La.eigen(sm, symmetric = FALSE); V2 <- em$vect
Line 121... Line 121...
121 
set.seed(123)
 121 
set.seed(123)
122 
sm <- matrix(rnorm(25), 5, 5)
 122 
sm <- matrix(rnorm(25), 5, 5)
123 
sm <- 0.5 * (sm + t(sm))
 123 
sm <- 0.5 * (sm + t(sm))
124 
eigenok(sm, eigen(sm))
 124 
eigenok(sm, eigen(sm))
125 
eigenok(sm, La.eigen(sm))
 125 
eigenok(sm, La.eigen(sm))
126 
eigenok(sm, La.eigen(sm, method="dsyevr"))
 126 
eigenok(sm, La.eigen(sm, method="dsyev"))
127 
eigenok(sm, La.eigen(sm, sym=FALSE))
 127 
eigenok(sm, La.eigen(sm, sym=FALSE))
128 
 
 128 
 
129 
sm[] <- as.complex(sm)
 129 
sm[] <- as.complex(sm)
130 
Ceigenok(sm, eigen(sm))
 130 
Ceigenok(sm, eigen(sm))
131 
Ceigenok(sm, La.eigen(sm))
 131 
Ceigenok(sm, La.eigen(sm))