1#  File src/library/base/R/eigen.R
2#  Part of the R package, https://www.R-project.org
3#
4#  Copyright (C) 1995-2019 The R Core Team
5#
6#  This program is free software; you can redistribute it and/or modify
7#  it under the terms of the GNU General Public License as published by
8#  the Free Software Foundation; either version 2 of the License, or
9#  (at your option) any later version.
10#
11#  This program is distributed in the hope that it will be useful,
12#  but WITHOUT ANY WARRANTY; without even the implied warranty of
13#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14#  GNU General Public License for more details.
15#
16#  A copy of the GNU General Public License is available at
17#  https://www.R-project.org/Licenses/
18
19
20isSymmetric <- function(object, ...) UseMethod("isSymmetric")
21
22isSymmetric.matrix <- function(object, tol = 100*.Machine$double.eps, tol1 = 8*tol, ...)
23{
24    if(!is.matrix(object)) return(FALSE) ## we test for  symmetric *matrix*
25    ## cheap pretest: is it square?
26    d <- dim(object)
27    if((n <- d[1L]) != d[2L]) return(FALSE)
28    iCplx <- is.complex(object)
29    if(n > 1L && length(tol1)) {
30	## initial pre-tests, fast for large non-symmetric:
31	Cj <- if(iCplx) Conj else identity
32	for(i in unique(c(1L, 2L, n-1L, n)))
33	    if(is.character(all.equal(object[i, ], Cj(object[, i]), tolerance = tol1, ...)))
34		return(FALSE)
35    }
36    test <-
37        if(iCplx)
38            all.equal.numeric(object, Conj(t(object)), tolerance = tol, ...)
39        else # numeric, character, ..
40            all.equal(object, t(object), tolerance = tol, ...)
41    isTRUE(test)
42}
43
44eigen <- function(x, symmetric, only.values = FALSE, EISPACK = FALSE)
45{
46    x <- unname(as.matrix(x))
47    n <- nrow(x)
48    if (!n) stop("0 x 0 matrix")
49    if (n != ncol(x)) stop("non-square matrix in 'eigen'")
50    n <- as.integer(n)
51    if(is.na(n)) stop("invalid nrow(x)")
52
53    complex.x <- is.complex(x)
54    if (!all(is.finite(x))) stop("infinite or missing values in 'x'")
55
56    if(missing(symmetric)) symmetric <- isSymmetric.matrix(x)
57
58    if (symmetric) {
59        z <- if(!complex.x) .Internal(La_rs(x, only.values))
60        else .Internal(La_rs_cmplx(x, only.values))
61        ord <- rev(seq_along(z$values))
62    } else {
63        z <- if(!complex.x) .Internal(La_rg(x, only.values))
64        else .Internal(La_rg_cmplx(x, only.values))
65        ord <- sort.list(Mod(z$values), decreasing = TRUE)
66    }
67    if(only.values)
68	list(values = z$values[ord], vectors = NULL)
69    else
70	structure(class = "eigen",
71		  list(values = z$values[ord],
72		       vectors = z$vectors[, ord, drop = FALSE]))
73}
74
75print.eigen <- function (x, ...) {
76    cat("eigen() decomposition\n")
77    print(unclass(x), ...)
78    invisible(x)
79}
80