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