The R Project SVN R

R version 2.6.0 Under development (unstable) (2007-07-24 r42292)
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(signif(E$values, 5))
+     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
>