Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
branch: tags/R-2-1-1
Fetching contributors…

Cannot retrieve contributors at this time

file 88 lines (68 sloc) 2.629 kb
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
## tests of R functions based on the lapack module

## ------- examples from ?svd using La.svd ---------

hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
Eps <- 100 * .Machine$double.eps

X <- hilbert(9)[,1:6]
str(s <- La.svd(X)); D <- diag(s$d)
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)# X = U D V'
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)# D = U' X V

str(s <- La.svd(X, method = "dgesvd")); D <- diag(s$d)
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)# X = U D V'
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)# D = U' X V

X <- cbind(1, 1:7)
str(s <- La.svd(X)); D <- diag(s$d)
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)# X = U D V'
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)# D = U' X V

X <- cbind(1, 1:7)
str(s <- La.svd(X, method = "dgesvd")); D <- diag(s$d)
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)# X = U D V'
stopifnot(abs(D - t(s$u) %*% X %*% t(s$vt)) < Eps)# D = U' X V

# test nu and nv
La.svd(X, nu = 0)
(s <- La.svd(X, nu = 7))
stopifnot(dim(s$u) == c(7,7))
La.svd(X, nv = 0)

La.svd(X, nu = 0, method = "dgesvd")
(s <- La.svd(X, nu = 7, method = "dgesvd"))
stopifnot(dim(s$u) == c(7,7))
La.svd(X, nv = 0, method = "dgesvd")

# test of complex case

X <- cbind(1, 1:7+(-3:3)*1i)
str(s <- La.svd(X)); D <- diag(s$d)
stopifnot(abs(X - s$u %*% D %*% s$vt) < Eps)
stopifnot(abs(D - Conj(t(s$u)) %*% X %*% Conj(t(s$vt))) < Eps)

# in this case svd calls La.svd
str(s <- svd(X)); D <- diag(s$d)
stopifnot(abs(X - s$u %*% D %*% Conj(t(s$v))) < Eps)
stopifnot(abs(D - Conj(t(s$u)) %*% X %*% s$v) < Eps)



## ------- tests of random real and complex matrices ------

## 100 may cause failures here.
eigenok <- function(A, E, Eps=1000*.Machine$double.eps)
{
    V <- E$vect; lam <- E$values
    stopifnot(abs(A %*% V - V %*% diag(lam)) < Eps,
              abs(lam[length(lam)]/lam[1]) < Eps || # this one not for singular A :
              abs(A - V %*% diag(lam) %*% t(V)) < Eps)
}

Ceigenok <- function(A, E, Eps=1000*.Machine$double.eps)
{
    V <- E$vect; lam <- E$values
    stopifnot(Mod(A %*% V - V %*% diag(lam)) < Eps,
              Mod(A - V %*% diag(lam) %*% Conj(t(V))) < Eps)
}

## failed for some 64bit-Lapack-gcc combinations:
sm <- cbind(1, 3:1, 1:3)
eigenok(sm, eigen(sm))
eigenok(sm, eigen(sm, sym=FALSE))

set.seed(123)
sm <- matrix(rnorm(25), 5, 5)
sm <- 0.5 * (sm + t(sm))
eigenok(sm, eigen(sm))
eigenok(sm, eigen(sm, sym=FALSE))

sm[] <- as.complex(sm)
Ceigenok(sm, eigen(sm))
Ceigenok(sm, eigen(sm, sym=FALSE))

sm[] <- sm + rnorm(25) * 1i
sm <- 0.5 * (sm + Conj(t(sm)))
Ceigenok(sm, eigen(sm))
Ceigenok(sm, eigen(sm, sym=FALSE))
Something went wrong with that request. Please try again.