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
|
#### Triangular Sparse Matrices in compressed column-oriented format
setAs("dtCMatrix", "ltCMatrix",
function(from) new("ltCMatrix", i = from@i, p = from@p,
uplo = from@uplo, diag = from@diag,
x = as.logical(from@x),
## FIXME?: use from@factors smartly
Dim = from@Dim, Dimnames = from@Dimnames))
setAs("dtCMatrix", "ntCMatrix", # just drop 'x' slot:
function(from) new("ntCMatrix", i = from@i, p = from@p,
uplo = from@uplo, diag = from@diag,
## FIXME?: use from@factors smartly
Dim = from@Dim, Dimnames = from@Dimnames))
setAs("matrix", "dtCMatrix",
function(from) as(as(from, "dtTMatrix"), "dtCMatrix"))
setAs("dtCMatrix", "dgCMatrix",
function(from) {
if (from@diag == "U")
from <- .Call(Csparse_diagU2N, from)
## new("dgCMatrix",
## i = from@i, p = from@p, x = from@x,
## Dim = from@Dim, Dimnames = from@Dimnames)
## ---> Rather faster, no checking:
copyClass(from, "dgCMatrix",
sNames = c("i", "p", "x", "Dim", "Dimnames"), check = FALSE)
})
setAs("dtCMatrix", "dsCMatrix", function(from) as(from, "symmetricMatrix"))
setAs("dtCMatrix", "dgTMatrix",
function(from) {
if (from@diag == "U") from <- .Call(Csparse_diagU2N, from)
## ignore triangularity in conversion to TsparseMatrix
.Call(Csparse_to_Tsparse, from, FALSE)
})
## FIXME: make more efficient
## ----- and as(., "triangularMatrix") is even worse via as_Sp()
setAs("dgCMatrix", "dtCMatrix", # to triangular, needed for triu,..
function(from) as(.Call(Csparse_to_Tsparse, from, FALSE), "dtCMatrix"))
setAs("dtCMatrix", "dgeMatrix",
function(from) as(as(from, "dgTMatrix"), "dgeMatrix"))
## These are all needed because cholmod doesn't support triangular:
## (see end of ./Csparse.R ), e.g. for triu()
setAs("dtCMatrix", "dtTMatrix",
function(from) .Call(Csparse_to_Tsparse, from, TRUE))
## {# and this is not elegant:
## x <- as(from, "dgTMatrix")
## if (from@diag == "U") { ## drop diagonal entries '1':
## i <- x@i; j <- x@j
## nonD <- i != j
## xx <- x@x[nonD] ; i <- i[nonD] ; j <- j[nonD]
## } else {
## xx <- x@x; i <- x@i; j <- x@j
## }
## new("dtTMatrix", x = xx, i = i, j = j, Dim = x@Dim,
## Dimnames = x@Dimnames, uplo = from@uplo, diag = from@diag)
## })
## Now that we support triangular matrices use the inherited method.
## setAs("dtCMatrix", "TsparseMatrix", function(from) as(from, "dtTMatrix"))
setAs("dtCMatrix", "dtrMatrix",
function(from) as(as(from, "dtTMatrix"), "dtrMatrix"))
setMethod("determinant", signature(x = "dtCMatrix", logarithm = "logical"),
function(x, logarithm = TRUE, ...) {
if(x@diag == "N")
mkDet(diag(x), logarithm)
else
structure(list(modulus = structure(if (logarithm) 0 else 1,
"logarithm" = logarithm),
sign = 1L),
class = "det")
})
setMethod("solve", signature(a = "dtCMatrix", b = "missing"),
function(a, b, ...) {
stopifnot((n <- nrow(a)) == ncol(a))
as(.Call(dtCMatrix_sparse_solve, a, .trDiagonal(n, unitri=FALSE)),
"dtCMatrix")
}, valueClass = "dtCMatrix")
setMethod("solve", signature(a = "dtCMatrix", b = "dgeMatrix"),
function(a, b, ...) .Call(dtCMatrix_matrix_solve, a, b, TRUE),
valueClass = "dgeMatrix")
setMethod("solve", signature(a = "dtCMatrix", b = "CsparseMatrix"),
function(a, b, ...) .sortCsparse(.Call(dtCMatrix_sparse_solve, a, b)),
## ------------ TODO: both in C code
valueClass = "dgCMatrix")
setMethod("solve", signature(a = "dtCMatrix", b = "matrix"),
function(a, b, ...) {
storage.mode(b) <- "double"
.Call(dtCMatrix_matrix_solve, a, b, FALSE)
}, valueClass = "dgeMatrix")
## Isn't this case handled by the method for (a = "Matrix', b =
## "numeric") in ./Matrix.R? Or is this method defined here for
## the as.double coercion?
setMethod("solve", signature(a = "dtCMatrix", b = "numeric"),
function(a, b, ...) .Call(dtCMatrix_matrix_solve, a,
cbind(as.double(b), deparse.level=0L), FALSE),
valueClass = "dgeMatrix")
if(FALSE)## still not working
setMethod("diag", "dtCMatrix",
function(x, nrow, ncol) .Call(diag_tC, x, "diag"))
## no pivoting here, use L or U
setMethod("lu", "dtCMatrix",
function(x, ...) {
n <- (d <- x@Dim)[1L]
p <- 0:(n-1L)
if(x@uplo == "U")
new("sparseLU",
L = .trDiagonal(n, uplo="L"),
U = x,
p = p, q = p, Dim = d)
else { ## "L" : x = L = L I
d <- diag(x)
new("sparseLU",
L = x %*% Diagonal(n, 1/d),
U = .trDiagonal(n, x = d),
p = p, q = p, Dim = d)
}
})
|
