Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
branch: tags/R-uthread…
Fetching contributors…

Cannot retrieve contributors at this time

file 138 lines (104 sloc) 4.115 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 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
## 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)


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

La.eigen(cbind(c(1,-1),c(-1,1)))
La.eigen(cbind(c(1,-1),c(-1,1)), symmetric = FALSE)

La.eigen(cbind(1,c(1,-1)), only.values = TRUE)
La.eigen(cbind(-1,2:1)) # complex values
La.eigen(print(cbind(c(0,1i), c(-1i,0))))# Hermite ==> real eigenvalues
## 3 x 3:
La.eigen(cbind( 1,3:1,1:3))
La.eigen(cbind(-1,c(1:2,0),0:2)) # complex values

La.eigen(cbind(c(1,-1),c(-1,1)), method = "dsyev")


set.seed(1234)
Meps <- .Machine$double.eps
m <- matrix(round(rnorm(25),3), 5,5)
sm <- m + t(m) #- symmetric matrix
em <- La.eigen(sm); V <- em$vect
print(lam <- em$values) # ordered DEcreasingly

stopifnot(
 abs(sm %*% V - V %*% diag(lam)) < 60*Meps,
 abs(sm - V %*% diag(lam) %*% t(V)) < 60*Meps)

em <- La.eigen(sm, method = "dsyev"); V <- em$vect
print(lam <- em$values) # ordered DEcreasingly

stopifnot(
 abs(sm %*% V - V %*% diag(lam)) < 60*Meps,
 abs(sm - V %*% diag(lam) %*% t(V)) < 60*Meps)

# check only.values = TRUE too
La.eigen(sm, only.values = TRUE)
La.eigen(sm, only.values = TRUE, method = "dsyev")


## symmetric = FALSE

em <- La.eigen(sm, symmetric = FALSE); V2 <- em$vect
print(lam2 <- em$values) # ordered decreasingly in ABSolute value !
print(i <- rev(order(lam2)))
stopifnot(abs(lam - lam2[i]) < 60 * Meps)

zapsmall(Diag <- t(V2) %*% V2) # orthonormal
print(norm2V <- apply(V2 * V2, 2, sum))
stopifnot( abs(1- norm2V / diag(Diag)) < 60*Meps)

stopifnot(abs(sm %*% V2 - V2 %*% diag(lam2)) < 60*Meps,
          abs(sm - V2 %*% diag(lam2) %*% t(V2)) < 60*Meps)


## ------- 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(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)
}

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

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

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