R Under development (unstable) (2023-10-11 r85316) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu 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 > > ## NB: the signs of singular and eigenvectors are arbitrary, > ## so there may be differences from the reference ouptut, > ## especially when alternative BLAS are used. > > options(digits = 4L) > tryCmsg <- function(expr) tryCatch(expr, error = conditionMessage) # typically == *$message > > ## ------- examples from ?svd --------- > > hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") } > Eps <- 100 * .Machine$double.eps > > ## The signs of the vectors are not determined here, so don't print > X <- hilbert(9L)[, 1:6] > s <- svd(X); D <- diag(s$d) > 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 > > ## ditto > X <- cbind(1, 1:7) > s <- svd(X); D <- diag(s$d) > 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 > s <- svd(X, nu = 0L) > s <- svd(X, nu = 7L) # the last 5 columns are not determined here > stopifnot(dim(s$u) == c(7L,7L)) > s <- svd(X, nv = 0L) > > # 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, 2L, function(x) x*sign(Re(x[1L]))) + 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.00000i 0.33381+0.0000i 0.02834+0.00000i 0.43783+0.0000i [2,] 0.3051+0.04099i -0.02643-0.1175i -0.43963+0.72556i -0.04739+0.2975i [3,] 0.3201-0.37556i 0.33790+0.4760i -0.09325-0.32814i 0.05364+0.2447i [4,] 0.3394+0.23303i -0.10443-0.6839i 0.09966-0.36289i 0.18940+0.1979i [5,] -0.2869+0.34830i -0.07660+0.2210i -0.14602+0.01322i 0.74490-0.1576i [,5] [1,] 0.6383+0.00000i [2,] -0.1909-0.20935i [3,] -0.4788-0.08610i [4,] -0.3654+0.04183i [5,] -0.2229-0.30121i [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.00000i 0.5373+0.00000i 0.428339+0.09065i 0.05039-0.329984i [2,] -0.1909-0.20935i 0.3051+0.04099i -0.107969+0.28126i -0.12013+0.008395i [3,] -0.4788-0.08610i 0.3201-0.37556i 0.001812+0.25051i 0.52156-0.262169i [4,] -0.3654+0.04183i 0.3394+0.23303i 0.144306+0.23287i -0.69180+0.000000i [5,] -0.2229-0.30121i -0.2869+0.34830i 0.761400+0.00000i 0.20693+0.109088i [,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 > > > ## ------- tests of integer matrices ----------------- > > set.seed(123) > A <- matrix(rpois(25, 5), 5, 5) > > A %*% A [,1] [,2] [,3] [,4] [,5] [1,] 202 170 156 160 234 [2,] 161 124 145 147 185 [3,] 166 136 134 130 174 [4,] 218 156 169 204 234 [5,] 205 134 175 181 249 > crossprod(A) [,1] [,2] [,3] [,4] [,5] [1,] 226 160 153 174 240 [2,] 160 143 126 112 179 [3,] 153 126 171 137 205 [4,] 174 112 137 174 192 [5,] 240 179 205 192 293 > tcrossprod(A) [,1] [,2] [,3] [,4] [,5] [1,] 229 155 150 207 184 [2,] 155 144 140 184 161 [3,] 150 140 156 176 142 [4,] 207 184 176 251 209 [5,] 184 161 142 209 227 > > solve(A) [,1] [,2] [,3] [,4] [,5] [1,] -0.048676 0.3390 -0.15756 -0.05892 -0.00854 [2,] -0.058711 -0.1262 0.19812 -0.03160 0.06426 [3,] 0.092656 0.2319 -0.02904 -0.12468 -0.09778 [4,] 0.062553 -0.1637 0.03800 -0.04270 0.12062 [5,] -0.002775 -0.2351 0.02391 0.22032 -0.02242 > qr(A) $qr [,1] [,2] [,3] [,4] [,5] [1,] -15.0333 -10.64304 -10.1774 -11.5743 -15.965 [2,] 0.4656 -5.45212 -3.2430 2.0516 -1.667 [3,] 0.2661 0.97998 -7.5434 -3.4278 -4.920 [4,] 0.5322 -0.05761 -0.3972 4.9068 -1.269 [5,] 0.5987 -0.17944 -0.9104 -0.9086 3.088 $rank [1] 5 $qraux [1] 1.266 1.064 1.116 1.418 3.088 $pivot [1] 1 2 3 4 5 attr(,"class") [1] "qr" > determinant(A, log = FALSE) $modulus [1] 9368 attr(,"logarithm") [1] FALSE $sign [1] -1 attr(,"class") [1] "det" > > rcond(A) [1] 0.02466 > rcond(A, "I") [1] 0.06007 > rcond(A, "1") [1] 0.02466 > > eigen(A) eigen() decomposition $values [1] 29.660+0.000i -4.631+0.000i 4.556+0.000i -2.292+3.117i -2.292-3.117i $vectors [,1] [,2] [,3] [,4] [,5] [1,] 0.4698+0i 0.34581+0i 0.07933+0i 0.72463+0.0000i 0.72463+0.0000i [2,] 0.3885+0i -0.30397+0i 0.31927+0i -0.24806-0.2591i -0.24806+0.2591i [3,] 0.3760+0i -0.03566+0i 0.70330+0i 0.04957+0.3697i 0.04957-0.3697i [4,] 0.5029+0i -0.74239+0i -0.32651+0i -0.22618+0.1063i -0.22618-0.1063i [5,] 0.4839+0i 0.48539+0i -0.53901+0i -0.33752-0.1752i -0.33752+0.1752i > ## The signs of the 'u' and 'v/vt' components can vary in the next two > A0 <- svd(A) > A1 <- La.svd(A) > ## OK to test == as these are the same Fortran calls. > stopifnot(A1$d == A0$d, A1$u == A0$u, A1$vt == t(A0$v)) > ## Fix the signs before printing. > s <- rep(sign(A0$u[1,]), each=5); A0$u <- s * A0$u; A0$v <- s * A0$v > A0 $d [1] 29.929 6.943 6.668 3.960 1.707 $u [,1] [,2] [,3] [,4] [,5] [1,] 0.4659 0.04141 0.8795 8.440e-02 0.02379 [2,] 0.3931 -0.15031 -0.2249 1.824e-05 0.87881 [3,] 0.3811 -0.62205 -0.2170 5.570e-01 -0.33239 [4,] 0.5166 -0.14296 -0.1852 -7.659e-01 -0.30290 [5,] 0.4652 0.75386 -0.3073 3.097e-01 -0.15780 $v [,1] [,2] [,3] [,4] [,5] [1,] 0.4831 0.3264 -0.47567 -0.1955 0.6289 [2,] 0.3627 -0.3731 -0.53450 0.5920 -0.3051 [3,] 0.3995 -0.4779 0.59210 0.2252 0.4590 [4,] 0.3983 0.6985 0.36298 0.3821 -0.2752 [5,] 0.5628 -0.1948 0.07549 -0.6439 -0.4743 > > > As <- crossprod(A) > E <- eigen(As) > E$values [1] 895.737 48.201 44.468 15.678 2.915 > abs(E$vectors) # signs vary [,1] [,2] [,3] [,4] [,5] [1,] 0.4831 0.3264 0.47567 0.1955 0.6289 [2,] 0.3627 0.3731 0.53450 0.5920 0.3051 [3,] 0.3995 0.4779 0.59210 0.2252 0.4590 [4,] 0.3983 0.6985 0.36298 0.3821 0.2752 [5,] 0.5628 0.1948 0.07549 0.6439 0.4743 > chol(As) [,1] [,2] [,3] [,4] [,5] [1,] 15.03 10.643 10.177 11.574 15.965 [2,] 0.00 5.452 3.243 -2.052 1.667 [3,] 0.00 0.000 7.543 3.428 4.920 [4,] 0.00 0.000 0.000 4.907 -1.269 [5,] 0.00 0.000 0.000 0.000 3.088 > backsolve(As, 1:5) [1] -0.009040 -0.005129 -0.006246 0.004158 0.017065 > > ## ------- tests of logical matrices ----------------- > > set.seed(123) > A <- matrix(runif(25) > 0.5, 5, 5) > > A %*% A [,1] [,2] [,3] [,4] [,5] [1,] 2 2 2 1 3 [2,] 2 1 1 2 3 [3,] 2 2 1 1 3 [4,] 2 2 2 2 4 [5,] 2 1 2 2 3 > crossprod(A) [,1] [,2] [,3] [,4] [,5] [1,] 3 2 1 1 3 [2,] 2 3 2 0 3 [3,] 1 2 3 1 3 [4,] 1 0 1 2 2 [5,] 3 3 3 2 5 > tcrossprod(A) [,1] [,2] [,3] [,4] [,5] [1,] 3 1 2 2 2 [2,] 1 3 2 3 2 [3,] 2 2 3 3 1 [4,] 2 3 3 4 2 [5,] 2 2 1 2 3 > > Q <- qr(A) > zapsmall(Q$qr) [,1] [,2] [,3] [,4] [,5] [1,] -1.7321 -1.1547 -0.5774 -0.5774 -1.7321 [2,] 0.5774 -1.2910 -1.0328 0.5164 -0.7746 [3,] 0.0000 0.7746 -1.2649 -0.9487 -0.9487 [4,] 0.5774 0.2582 0.0508 0.7071 0.7071 [5,] 0.5774 -0.5164 -0.6803 -0.3136 0.0000 > zapsmall(Q$qraux) [1] 1.000 1.258 1.731 1.950 0.000 > determinant(A, log = FALSE) # 0 $modulus [1] 0 attr(,"logarithm") [1] FALSE $sign [1] 1 attr(,"class") [1] "det" > > rcond(A) [1] 0 > rcond(A, "I") [1] 0 > rcond(A, "1") [1] 0 > > E <- eigen(A) > zapsmall(E$values) [1] 3.163+0.000i 0.271+0.908i 0.271-0.908i -0.705+0.000i 0.000+0.000i > zapsmall(Mod(E$vectors)) [,1] [,2] [,3] [,4] [,5] [1,] 0.4358 0.3604 0.3604 0.6113 0.0000 [2,] 0.4087 0.4495 0.4495 0.3771 0.5774 [3,] 0.3962 0.5870 0.5870 0.2028 0.0000 [4,] 0.5340 0.2792 0.2792 0.6649 0.5774 [5,] 0.4483 0.4955 0.4955 0.0314 0.5774 > S <- svd(A) > zapsmall(S$d) [1] 3.379 1.536 1.414 0.472 0.000 > S <- La.svd(A) > zapsmall(S$d) [1] 3.379 1.536 1.414 0.472 0.000 > > As <- A > As[upper.tri(A)] <- t(A)[upper.tri(A)] > det(As) [1] 2 > E <- eigen(As) > E$values [1] 3.465 1.510 0.300 -1.000 -1.275 > ## The eigenvectors are of arbitrary sign, so we fix the first element to > ## be positive for cross-platform comparisons. > Ev <- E$vectors > zapsmall(Ev * rep(sign(Ev[1, ]), each = 5)) [,1] [,2] [,3] [,4] [,5] [1,] 0.4023 0.2877 0.3638 0.5 0.6108 [2,] 0.5338 -0.3474 0.5742 -0.5 -0.1207 [3,] 0.4177 -0.5380 -0.6384 0.0 0.3585 [4,] 0.4959 0.0733 -0.1273 0.5 -0.6946 [5,] 0.3644 0.7084 -0.3378 -0.5 0.0369 > solve(As) [,1] [,2] [,3] [,4] [,5] [1,] 0 1 -1.0 0.0 0.0 [2,] 1 1 -1.0 0.0 -1.0 [3,] -1 -1 1.5 0.5 0.5 [4,] 0 0 0.5 -0.5 0.5 [5,] 0 -1 0.5 0.5 0.5 > > ## quite hard to come up with an example where this might make sense. > Ac <- A; Ac[] <- as.logical(diag(5)) > chol(Ac) [,1] [,2] [,3] [,4] [,5] [1,] 1 0 0 0 0 [2,] 0 1 0 0 0 [3,] 0 0 1 0 0 [4,] 0 0 0 1 0 [5,] 0 0 0 0 1 > > ## ------- tests of non-finite values ----------------- > > a <- matrix(NaN, 3, 3,, list(one=1:3, two=letters[1:3])) > b <- cbind(1:3, NA) > dimnames(b) <- list(One=4:6, Two=11:12) > bb <- 1:3; names(bb) <- 11:12 > ## gave error with LAPACK 3.11.0 > ## names(dimnames(.)), ("two", "Two") are lost {FIXME?}: > ## IGNORE_RDIFF_BEGIN > stopifnot(is.na(print(solve(a, b )))) # is.na(): NA *or* NaN 11 12 a NaN NA b NaN NA c NaN NA > ## IGNORE_RDIFF_END > stopifnot(is.na(print(solve(a, bb)))) # all NaN a b c NaN NaN NaN > > A <- a + 0i > A_b <- solve(A, b) # platform dependent result (e.g. OPENBLAS ..) > stopifnot(is.na(A_b)) > ## IGNORE_RDIFF_BEGIN > A_b 11 12 a NaN+NaNi NA b NaN+NaNi NA c NaN+NaNi NA > rbind(re = Re(A_b[,2]), im = Im(A_b[,2])) # often was "all NA", now typically "re=NA, im=NaN" a b c re NA NA NA im NaN NaN NaN > ## IGNORE_RDIFF_END > > > ## PR#18541 by Mikael Jagan -- chol() error & warning message: > x <- diag(-1, 5L) > (chF <- tryCmsg(chol(x, pivot = FALSE))) # dpotrf [1] "the leading minor of order 1 is not positive" > (chT <- withCallingHandlers(warning = function(w) ..W <<- conditionMessage(w), + chol(x, pivot = TRUE ))) # dpstrf [,1] [,2] [,3] [,4] [,5] [1,] -1 0 0 0 0 [2,] 0 -1 0 0 0 [3,] 0 0 -1 0 0 [4,] 0 0 0 -1 0 [5,] 0 0 0 0 -1 attr(,"pivot") [1] 1 2 3 4 5 attr(,"rank") [1] 0 Warning message: In chol.default(x, pivot = TRUE) : the matrix is either rank-deficient or not positive definite > stopifnot(exprs = { + grepl(" minor .* not positive$", chF) # was "not positive *definite* + grepl("rank-deficient or not positive definite$", ..W) # was "indefinite* + ## platform dependent, Mac has several NaN's chT == -diag(5) + attr(chT, "rank") %in% 0:1 + }) >