1## tests of R functions based on the lapack module
2
3## NB: the signs of singular and eigenvectors are arbitrary,
4## so there may be differences from the reference ouptut,
5## especially when alternative BLAS are used.
6
7options(digits = 4L)
8
9##    -------  examples from ?svd ---------
10
11hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
12Eps <- 100 * .Machine$double.eps
13
14## The signs of the vectors are not determined here, so don't print
15X <- hilbert(9L)[, 1:6]
16s <- svd(X); D <- diag(s$d)
17stopifnot(abs(X - s$u %*% D %*% t(s$v)) < Eps)#  X = U D V'
18stopifnot(abs(D - t(s$u) %*% X %*% s$v) < Eps)#  D = U' X V
19
20## ditto
21X <- cbind(1, 1:7)
22s <- svd(X); D <- diag(s$d)
23stopifnot(abs(X - s$u %*% D %*% t(s$v)) < Eps)#  X = U D V'
24stopifnot(abs(D - t(s$u) %*% X %*% s$v) < Eps)#  D = U' X V
25
26# test nu and nv
27s <- svd(X, nu = 0L)
28s <- svd(X, nu = 7L) # the last 5 columns are not determined here
29stopifnot(dim(s$u) == c(7L,7L))
30s <- svd(X, nv = 0L)
31
32# test of complex case
33
34X <- cbind(1, 1:7+(-3:3)*1i)
35s <- svd(X); D <- diag(s$d)
36stopifnot(abs(X - s$u %*% D %*% Conj(t(s$v))) < Eps)
37stopifnot(abs(D - Conj(t(s$u)) %*% X %*% s$v) < Eps)
38
39
40
41##  -------  tests of random real and complex matrices ------
42fixsign <- function(A) {
43    A[] <- apply(A, 2L, function(x) x*sign(Re(x[1L])))
44    A
45}
46##			       100  may cause failures here.
47eigenok <- function(A, E, Eps=1000*.Machine$double.eps)
48{
49    print(fixsign(E$vectors))
50    print(zapsmall(E$values))
51    V <- E$vectors; lam <- E$values
52    stopifnot(abs(A %*% V - V %*% diag(lam)) < Eps,
53              abs(lam[length(lam)]/lam[1]) < Eps | # this one not for singular A :
54              abs(A - V %*% diag(lam) %*% t(V)) < Eps)
55}
56
57Ceigenok <- function(A, E, Eps=1000*.Machine$double.eps)
58{
59    print(fixsign(E$vectors))
60    print(signif(E$values, 5))
61    V <- E$vectors; lam <- E$values
62    stopifnot(Mod(A %*% V - V %*% diag(lam)) < Eps,
63              Mod(A - V %*% diag(lam) %*% Conj(t(V))) < Eps)
64}
65
66## failed for some 64bit-Lapack-gcc combinations:
67sm <- cbind(1, 3:1, 1:3)
68eigenok(sm, eigen(sm))
69eigenok(sm, eigen(sm, sym=FALSE))
70
71set.seed(123)
72sm <- matrix(rnorm(25), 5, 5)
73sm <- 0.5 * (sm + t(sm))
74eigenok(sm, eigen(sm))
75eigenok(sm, eigen(sm, sym=FALSE))
76
77sm[] <- as.complex(sm)
78Ceigenok(sm, eigen(sm))
79Ceigenok(sm, eigen(sm, sym=FALSE))
80
81sm[] <- sm + rnorm(25) * 1i
82sm <- 0.5 * (sm + Conj(t(sm)))
83Ceigenok(sm, eigen(sm))
84Ceigenok(sm, eigen(sm, sym=FALSE))
85
86
87##  -------  tests of integer matrices -----------------
88
89set.seed(123)
90A <- matrix(rpois(25, 5), 5, 5)
91
92A %*% A
93crossprod(A)
94tcrossprod(A)
95
96solve(A)
97qr(A)
98determinant(A, log = FALSE)
99
100rcond(A)
101rcond(A, "I")
102rcond(A, "1")
103
104eigen(A)
105svd(A)
106La.svd(A)
107
108As <- crossprod(A)
109E <- eigen(As)
110E$values
111abs(E$vectors) # signs vary
112chol(As)
113backsolve(As, 1:5)
114
115##  -------  tests of logical matrices -----------------
116
117set.seed(123)
118A <- matrix(runif(25) > 0.5, 5, 5)
119
120A %*% A
121crossprod(A)
122tcrossprod(A)
123
124Q <- qr(A)
125zapsmall(Q$qr)
126zapsmall(Q$qraux)
127determinant(A, log = FALSE) # 0
128
129rcond(A)
130rcond(A, "I")
131rcond(A, "1")
132
133E <- eigen(A)
134zapsmall(E$values)
135zapsmall(Mod(E$vectors))
136S <- svd(A)
137zapsmall(S$d)
138S <- La.svd(A)
139zapsmall(S$d)
140
141As <- A
142As[upper.tri(A)] <- t(A)[upper.tri(A)]
143det(As)
144E <- eigen(As)
145E$values
146## The eigenvectors are of arbitrary sign, so we fix the first element to
147## be positive for cross-platform comparisons.
148Ev <- E$vectors
149zapsmall(Ev * rep(sign(Ev[1, ]), each = 5))
150solve(As)
151
152## quite hard to come up with an example where this might make sense.
153Ac <- A; Ac[] <- as.logical(diag(5))
154chol(Ac)
155
156
157