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