1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
#### eigen() , Schur() etc
#### ===== =====
## eigen() is not even generic, and we haven't any C code,
## NOTE base::eigen() "magically" can work via as.matrix()
if(.Matrix.avoiding.as.matrix) {
## ---- IFF as.matrix(.) <==> as(., "matrix") [which we consider _deprecating_]
## FIXME: Code for *sparse* !! [RcppEigen ~??~]
setMethod("eigen", signature(x = "Matrix", only.values = "missing"),
function(x, symmetric, only.values, EISPACK) # << must match generic
base::eigen(as(x,"matrix"), symmetric, FALSE))
setMethod("eigen", signature(x = "Matrix", only.values = "logical"),
function(x, symmetric, only.values, EISPACK)
base::eigen(as(x,"matrix"), symmetric, only.values))
## base::svd() using as.matrix() := asRbasematrix()
if(getRversion() < "3.5.0") # svd not yet implicit generic
setGeneric("svd", function(x, ...) base::svd(x, ...))
setMethod("svd", "Matrix",
function (x, ...) base::svd(as(x,"matrix"), ...))
}
.dgeSchur <- function(x, vectors, ...) {
cl <- .Call(dgeMatrix_Schur, x, TRUE, TRUE)
realEV <- all(cl$WI == 0)
## TODO: do all this in C
new("Schur", Dim = x@Dim,
Q = as(cl$Z, "dgeMatrix"),
T = as(cl$T, if(realEV)"dtrMatrix" else "dgeMatrix"),
EValues = if(realEV) cl$WR else complex(real = cl$WR, imaginary = cl$WI))
}
setMethod("Schur", signature(x = "dgeMatrix", vectors = "missing"),
.dgeSchur)
setMethod("Schur", signature(x = "dgeMatrix", vectors = "logical"),
function(x, vectors, ...) {
if(vectors) .dgeSchur(x)
else {
cl <- .Call(dgeMatrix_Schur, x, FALSE, TRUE)
realEV <- all(cl$WI == 0)
list(T = as(cl$T, if(realEV) "dtrMatrix" else "dgeMatrix"),
EValues =
if(realEV) cl$WR else complex(real = cl$WR, imaginary = cl$WI))
}})
## Ok, for the faint of heart, also provide "matrix" methods :
.mSchur <- function(x, vectors, ...) {
cl <- .Call(dgeMatrix_Schur, x, TRUE, FALSE)
list(Q = cl$Z,
T = cl$T,
EValues = if(all(cl$WI == 0)) cl$WR
else complex(real = cl$WR, imaginary = cl$WI))
}
setMethod("Schur", signature(x = "matrix", vectors = "missing"), .mSchur)
setMethod("Schur", signature(x = "matrix", vectors = "logical"),
function(x, vectors, ...) {
if(vectors) .mSchur(x)
else {
cl <- .Call(dgeMatrix_Schur, x, FALSE, FALSE)
EV <- if(all(cl$WI == 0)) cl$WR
else complex(real = cl$WR, imaginary = cl$WI)
cl$WR <- cl$WI <- NULL
cl$EValues <- EV
cl
}})
Schur.dsy <- function(x, vectors, ...)
{
if(missing(vectors)) vectors <- TRUE
## TODO: do all this in C
## Should directly call LAPACK dsyev()
evl <- eigen(x, only.values = !vectors)
eVals <- evl$values
if(vectors)
new("Schur", Dim = x@Dim,
Q = as(evl$vectors, "dgeMatrix"),
T = Diagonal(x = eVals),
EValues = eVals)
else
list(T = Diagonal(x = eVals), EValues = eVals)
}
setMethod("Schur", signature(x = "dsyMatrix", vectors = "ANY"), Schur.dsy)
## FIXME(?) these coerce from sparse to *dense*
setMethod("Schur", signature(x = "generalMatrix", vectors = "missing"),
function(x, vectors, ...) callGeneric(as(x, "dgeMatrix")))
setMethod("Schur", signature(x = "generalMatrix", vectors = "logical"),
function(x, vectors, ...) callGeneric(as(x, "dgeMatrix"), vectors))
setMethod("Schur", signature(x = "symmetricMatrix", vectors = "missing"),
function(x, vectors, ...) Schur.dsy(as(x, "dsyMatrix")))
setMethod("Schur", signature(x = "symmetricMatrix", vectors = "logical"),
function(x, vectors, ...) Schur.dsy(as(x, "dsyMatrix"), vectors))
## Schur(<diagonal>) : {Note that the Schur decomposition is not unique here}
.simpleSchur <- function(x, vectors, ...) {
x <- as(x, "dMatrix")
d <- dim(x)
new("Schur", Dim = d, Q = Diagonal(d[1]), T = x, EValues = diag(x))
}
setMethod("Schur", signature(x = "diagonalMatrix", vectors = "missing"),
.simpleSchur)
setMethod("Schur", signature(x = "diagonalMatrix", vectors = "logical"),
function(x, vectors, ...) {
if(vectors) .simpleSchur(x)
else {
x <- as(x, "dMatrix")
list(T = x, EValues = x@x)
}})
.triSchur <- function(x, vectors, ...) {
x <- as(x, "dMatrix")
d <- dim(x)
n <- d[1]
if(x@uplo == "U" || n == 0)
new("Schur", Dim = d, Q = Diagonal(n), T = x, EValues = diag(x))
else {
i <- n:1
new("Schur", Dim = d, Q = as(i, "pMatrix"),
T = t(t(x)[i,i]), EValues = diag(x)[i])
}
}
setMethod("Schur", signature(x = "triangularMatrix", vectors = "missing"),
.triSchur)
setMethod("Schur", signature(x = "triangularMatrix", vectors = "logical"),
function(x, vectors, ...) {
if(vectors) .triSchur(x)
else {
x <- as(x, "dMatrix")
list(T = x, EValues = x@x)
}})
|
