R version 2.6.0 Under development (unstable) (2007-05-04 r41441) Copyright (C) 2007 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## tests of R functions based on the lapack module > > options(digits=4) > > ## ------- examples from ?svd --------- > > hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") } > Eps <- 100 * .Machine$double.eps > > X <- hilbert(9)[,1:6] > (s <- svd(X)); D <- diag(s$d) $d [1] 1.668e+00 2.774e-01 2.224e-02 1.085e-03 3.244e-05 5.235e-07 $u [,1] [,2] [,3] [,4] [,5] [,6] [1,] -0.7245 0.6266 0.27350 -0.08527 0.02074 -0.004025 [2,] -0.4282 -0.1299 -0.64294 0.55047 -0.27253 0.092816 [3,] -0.3122 -0.2804 -0.33633 -0.31418 0.61632 -0.440904 [4,] -0.2479 -0.3142 -0.06931 -0.44667 0.02945 0.530120 [5,] -0.2064 -0.3141 0.10786 -0.30242 -0.35567 0.237038 [6,] -0.1771 -0.3027 0.22106 -0.09042 -0.38879 -0.260449 [7,] -0.1553 -0.2877 0.29281 0.11551 -0.19286 -0.420945 [8,] -0.1384 -0.2722 0.33784 0.29313 0.11633 -0.160790 [9,] -0.1249 -0.2571 0.36543 0.43885 0.46497 0.434600 $v [,1] [,2] [,3] [,4] [,5] [,6] [1,] -0.7365 0.6225 0.2550 -0.06976 0.01328 -0.001588 [2,] -0.4433 -0.1819 -0.6867 0.50860 -0.19627 0.041117 [3,] -0.3275 -0.3509 -0.2611 -0.50474 0.61606 -0.259216 [4,] -0.2626 -0.3922 0.1044 -0.43748 -0.40834 0.638902 [5,] -0.2204 -0.3946 0.3510 0.01612 -0.46428 -0.675827 [6,] -0.1904 -0.3832 0.5111 0.53856 0.44664 0.257249 > stopifnot(abs(X - s$u %*% D %*% t(s$v)) < Eps)# X = U D V' > stopifnot(abs(D - t(s$u) %*% X %*% s$v) < Eps)# D = U' X V > > X <- cbind(1, 1:7) > (s <- svd(X)); D <- diag(s$d) $d [1] 12.07 1.16 $u [,1] [,2] [1,] 0.09762 0.6744 [2,] 0.17884 0.5037 [3,] 0.26006 0.3331 [4,] 0.34128 0.1624 [5,] 0.42250 -0.0082 [6,] 0.50372 -0.1788 [7,] 0.58494 -0.3495 $v [,1] [,2] [1,] 0.1979 0.9802 [2,] 0.9802 -0.1979 > stopifnot(abs(X - s$u %*% D %*% t(s$v)) < Eps)# X = U D V' > stopifnot(abs(D - t(s$u) %*% X %*% s$v) < Eps)# D = U' X V > > # test nu and nv > svd(X, nu = 0) $d [1] 12.07 1.16 $u NULL $v [,1] [,2] [1,] 0.1979 0.9802 [2,] 0.9802 -0.1979 > (s <- svd(X, nu = 7)) $d [1] 12.07 1.16 $u [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0.09762 0.6744 -0.38314 -0.35392 -0.32471 -0.29549 -0.26628 [2,] 0.17884 0.5037 -0.13857 0.05843 0.25543 0.45242 0.64942 [3,] 0.26006 0.3331 0.89399 -0.09282 -0.07963 -0.06644 -0.05326 [4,] 0.34128 0.1624 -0.10083 0.88314 -0.13289 -0.14891 -0.16494 [5,] 0.42250 -0.0082 -0.09566 -0.14090 0.81386 -0.23139 -0.27663 [6,] 0.50372 -0.1788 -0.09048 -0.16494 -0.23940 0.68614 -0.38831 [7,] 0.58494 -0.3495 -0.08531 -0.18898 -0.29265 -0.39633 0.50000 $v [,1] [,2] [1,] 0.1979 0.9802 [2,] 0.9802 -0.1979 > stopifnot(dim(s$u) == c(7,7)) > svd(X, nv = 0) $d [1] 12.07 1.16 $u [,1] [,2] [1,] 0.09762 0.6744 [2,] 0.17884 0.5037 [3,] 0.26006 0.3331 [4,] 0.34128 0.1624 [5,] 0.42250 -0.0082 [6,] 0.50372 -0.1788 [7,] 0.58494 -0.3495 $v NULL > > # test of complex case > > X <- cbind(1, 1:7+(-3:3)*1i) > 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 ------ > fixsign <- function(A) { + A[] <- apply(A, 2, function(x) x*sign(Re(x[1]))) + A + } > ## 100 may cause failures here. > eigenok <- function(A, E, Eps=1000*.Machine$double.eps) + { + print(fixsign(E$vectors)) + print(zapsmall(E$values)) + V <- E$vectors; 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) + { + print(fixsign(E$vectors)) + print(E$values) + V <- E$vectors; 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)) [,1] [,2] [,3] [1,] 0.5774 0.8452 0.9428 [2,] 0.5774 0.1690 -0.2357 [3,] 0.5774 -0.5071 -0.2357 [1] 5 1 0 > eigenok(sm, eigen(sm, sym=FALSE)) [,1] [,2] [,3] [1,] 0.5774 0.8452 0.9428 [2,] 0.5774 0.1690 -0.2357 [3,] 0.5774 -0.5071 -0.2357 [1] 5 1 0 > > set.seed(123) > sm <- matrix(rnorm(25), 5, 5) > sm <- 0.5 * (sm + t(sm)) > eigenok(sm, eigen(sm)) [,1] [,2] [,3] [,4] [,5] [1,] 0.5899 0.1683 0.02315 0.471808 0.6329 [2,] 0.1936 0.2931 0.89217 -0.009784 -0.2838 [3,] 0.6627 -0.4812 -0.15825 0.082550 -0.5454 [4,] 0.1404 0.7985 -0.41848 0.094314 -0.3983 [5,] -0.3946 -0.1285 0.05768 0.872692 -0.2507 [1] 1.7814 1.5184 0.5833 -1.0148 -2.4908 > eigenok(sm, eigen(sm, sym=FALSE)) [,1] [,2] [,3] [,4] [,5] [1,] 0.6329 0.5899 0.1683 0.471808 0.02315 [2,] -0.2838 0.1936 0.2931 -0.009784 0.89217 [3,] -0.5454 0.6627 -0.4812 0.082550 -0.15825 [4,] -0.3983 0.1404 0.7985 0.094314 -0.41848 [5,] -0.2507 -0.3946 -0.1285 0.872692 0.05768 [1] -2.4908 1.7814 1.5184 -1.0148 0.5833 > > sm[] <- as.complex(sm) > Ceigenok(sm, eigen(sm)) [,1] [,2] [,3] [,4] [,5] [1,] 0.5899+0i 0.1683+0i 0.02315+0i 0.471808+0i 0.6329+0i [2,] 0.1936+0i 0.2931+0i 0.89217+0i -0.009784+0i -0.2838+0i [3,] 0.6627+0i -0.4812+0i -0.15825+0i 0.082550+0i -0.5454+0i [4,] 0.1404+0i 0.7985+0i -0.41848+0i 0.094314+0i -0.3983+0i [5,] -0.3946+0i -0.1285+0i 0.05768+0i 0.872692+0i -0.2507+0i [1] 1.7814 1.5184 0.5833 -1.0148 -2.4908 > Ceigenok(sm, eigen(sm, sym=FALSE)) [,1] [,2] [,3] [,4] [,5] [1,] 0.6329+0i 0.5899+0i 0.1683+0i 0.471808+0i 0.02315+0i [2,] -0.2838+0i 0.1936+0i 0.2931+0i -0.009784+0i 0.89217+0i [3,] -0.5454+0i 0.6627+0i -0.4812+0i 0.082550+0i -0.15825+0i [4,] -0.3983+0i 0.1404+0i 0.7985+0i 0.094314+0i -0.41848+0i [5,] -0.2507+0i -0.3946+0i -0.1285+0i 0.872692+0i 0.05768+0i [1] -2.4908+0i 1.7814+0i 1.5184+0i -1.0148+0i 0.5833+0i > > sm[] <- sm + rnorm(25) * 1i > sm <- 0.5 * (sm + Conj(t(sm))) > Ceigenok(sm, eigen(sm)) [,1] [,2] [,3] [,4] [1,] 0.5373+0.0000i 0.3338+0.0000i 0.02834+0.0000i 0.4378+0.0000i [2,] 0.3051+0.0410i -0.0264-0.1175i -0.43963+0.7256i -0.0474+0.2975i [3,] 0.3201-0.3756i 0.3379+0.4760i -0.09325-0.3281i 0.0536+0.2447i [4,] 0.3394+0.2330i -0.1044-0.6839i 0.09966-0.3629i 0.1894+0.1979i [5,] -0.2869+0.3483i -0.0766+0.2210i -0.14602+0.0132i 0.7449-0.1576i [,5] [1,] 0.6383+0.0000i [2,] -0.1909-0.2093i [3,] -0.4788-0.0861i [4,] -0.3654+0.0418i [5,] -0.2229-0.3012i [1] 2.4043 1.3934 0.7854 -1.4050 -2.8006 > Ceigenok(sm, eigen(sm, sym=FALSE)) [,1] [,2] [,3] [,4] [1,] 0.6383+0.0000i 0.5373+0.0000i 0.4283+0.0907i 0.0504-0.3300i [2,] -0.1909-0.2093i 0.3051+0.0410i -0.1080+0.2813i -0.1201+0.0084i [3,] -0.4788-0.0861i 0.3201-0.3756i 0.0018+0.2505i 0.5216-0.2622i [4,] -0.3654+0.0418i 0.3394+0.2330i 0.1443+0.2329i -0.6918+0.0000i [5,] -0.2229-0.3012i -0.2869+0.3483i 0.7614+0.0000i 0.2069+0.1091i [,5] [1,] 0.01468+0.02424i [2,] -0.84836+0.00000i [3,] 0.23232-0.24980i [4,] 0.36200-0.10282i [5,] -0.08698-0.11804i [1] -2.8006+0i 2.4043-0i -1.4050+0i 1.3934+0i 0.7854+0i >