SCM

SCM Repository

[matrix] Annotation of /pkg/Matrix/R/diagMatrix.R
ViewVC logotype

Annotation of /pkg/Matrix/R/diagMatrix.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2123 - (view) (download)
Original Path: pkg/R/diagMatrix.R

1 : maechler 1617 #### All methods for "diagonalMatrix" and its subclasses,
2 :     #### currently "ddiMatrix", "ldiMatrix"
3 :    
4 : maechler 1109 ## Purpose: Constructor of diagonal matrices -- ~= diag() ,
5 :     ## but *not* diag() extractor!
6 :     Diagonal <- function(n, x = NULL)
7 :     {
8 : maechler 1575 ## Allow Diagonal(4) and Diagonal(x=1:5)
9 : maechler 1109 if(missing(n))
10 : maechler 1575 n <- length(x)
11 : maechler 1109 else {
12 : maechler 1575 stopifnot(length(n) == 1, n == as.integer(n), n >= 0)
13 :     n <- as.integer(n)
14 : maechler 1109 }
15 :    
16 : maechler 1654 if(missing(x)) ## unit diagonal matrix
17 : maechler 1575 new("ddiMatrix", Dim = c(n,n), diag = "U")
18 : maechler 1109 else {
19 : maechler 1575 stopifnot(length(x) == n)
20 :     if(is.logical(x))
21 :     cl <- "ldiMatrix"
22 :     else if(is.numeric(x)) {
23 :     cl <- "ddiMatrix"
24 :     x <- as.numeric(x)
25 :     }
26 :     else if(is.complex(x)) {
27 :     cl <- "zdiMatrix" # will not yet work
28 :     } else stop("'x' has invalid data type")
29 :     new(cl, Dim = c(n,n), diag = "N", x = x)
30 : maechler 1109 }
31 :     }
32 :    
33 : maechler 1617 ### This is modified from a post of Bert Gunter to R-help on 1 Sep 2005.
34 :     ### Bert's code built on a post by Andy Liaw who most probably was influenced
35 :     ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
36 :     ### who posted his bdiag() function written in December 1995.
37 :    
38 :     bdiag <- function(...) {
39 :     if(nargs() == 0) return(new("dgCMatrix"))
40 :     ## else :
41 :     mlist <- if (nargs() == 1) as.list(...) else list(...)
42 :     dims <- sapply(mlist, dim)
43 :     ## make sure we had all matrices:
44 :     if(!(is.matrix(dims) && nrow(dims) == 2))
45 :     stop("some arguments are not matrices")
46 : maechler 1845 csdim <- rbind(rep.int(0L, 2),
47 : maechler 1617 apply(sapply(mlist, dim), 1, cumsum))
48 :     ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))
49 :     add1 <- matrix(1:0, 2,2)
50 : maechler 1654 for(i in seq_along(mlist)) {
51 : maechler 1617 indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
52 :     if(is.null(dim(indx))) ## non-square matrix
53 :     ret[indx[[1]],indx[[2]]] <- mlist[[i]]
54 :     else ## square matrix
55 :     ret[indx[,1],indx[,2]] <- mlist[[i]]
56 :     }
57 :     ## slightly debatable if we really should return Csparse.. :
58 :     as(ret, "CsparseMatrix")
59 :     }
60 :    
61 : maechler 1845 diag2tT <- function(from) {
62 :     i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
63 : maechler 1654 new(paste(.M.kind(from), "tTMatrix", sep=''),
64 :     diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
65 :     x = from@x, # <- ok for diag = "U" and "N" (!)
66 :     i = i, j = i)
67 :     }
68 : maechler 1109
69 : maechler 1845 diag2sT <- function(from) { # to symmetric Tsparse
70 :     i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
71 :     new(paste(.M.kind(from), "sTMatrix", sep=''),
72 :     Dim = from@Dim, Dimnames = from@Dimnames,
73 :     x = from@x, i = i, j = i)
74 :     }
75 :    
76 :     setAs("diagonalMatrix", "triangularMatrix", diag2tT)
77 :     setAs("diagonalMatrix", "sparseMatrix", diag2tT)
78 : maechler 1805 ## needed too (otherwise <dense> -> Tsparse is taken):
79 : maechler 1845 setAs("diagonalMatrix", "TsparseMatrix", diag2tT)
80 : maechler 1654 ## is better than this:
81 :     ## setAs("diagonalMatrix", "sparseMatrix",
82 :     ## function(from)
83 :     ## as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))
84 :     setAs("diagonalMatrix", "CsparseMatrix",
85 : maechler 1845 function(from) as(diag2tT(from), "CsparseMatrix"))
86 : maechler 1654
87 : maechler 1845 setAs("diagonalMatrix", "symmetricMatrix", diag2sT)
88 :    
89 : maechler 1109 setAs("diagonalMatrix", "matrix",
90 :     function(from) {
91 :     n <- from@Dim[1]
92 :     diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE
93 :     } else from@x,
94 :     nrow = n, ncol = n)
95 :     })
96 :    
97 : maechler 2098 setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
98 :     function(x, mode) {
99 :     n <- x@Dim[1]
100 :     mod <- mode(x@x)
101 :     r <- vector(mod, length = n^2)
102 :     if(n)
103 :     r[1 + 0:(n - 1) * (n + 1)] <-
104 :     if(x@diag == "U")
105 :     switch(mod, "integer"= 1L,
106 :     "numeric"= 1, "logical"= TRUE)
107 :     else x@x
108 :     r
109 :     })
110 :    
111 : maechler 1654 setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
112 :     function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))
113 : maechler 1174
114 : maechler 1655 .diag.x <- function(m) {
115 :     if(m@diag == "U")
116 :     rep.int(if(is.numeric(m@x)) 1. else TRUE,
117 :     m@Dim[1])
118 :     else m@x
119 :     }
120 :    
121 :     .diag.2N <- function(m) {
122 :     if(m@diag == "U") m@diag <- "N"
123 :     m
124 :     }
125 :    
126 : maechler 1654 ## given the above, the following 4 coercions should be all unneeded;
127 :     ## we prefer triangular to general:
128 : maechler 1295 setAs("ddiMatrix", "dgTMatrix",
129 :     function(from) {
130 : maechler 1654 .Deprecated("as(, \"sparseMatrix\")")
131 : maechler 1295 n <- from@Dim[1]
132 : maechler 1845 i <- seq_len(n) - 1L
133 : maechler 1655 new("dgTMatrix", i = i, j = i, x = .diag.x(from),
134 : maechler 1295 Dim = c(n,n), Dimnames = from@Dimnames) })
135 :    
136 :     setAs("ddiMatrix", "dgCMatrix",
137 : maechler 1655 function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))
138 : maechler 1295
139 :     setAs("ldiMatrix", "lgTMatrix",
140 :     function(from) {
141 : maechler 1654 .Deprecated("as(, \"sparseMatrix\")")
142 : maechler 1295 n <- from@Dim[1]
143 : maechler 1575 if(from@diag == "U") { # unit-diagonal
144 :     x <- rep.int(TRUE, n)
145 : maechler 1845 i <- seq_len(n) - 1L
146 : maechler 1575 } else { # "normal"
147 :     nz <- nz.NA(from@x, na. = TRUE)
148 :     x <- from@x[nz]
149 : maechler 1845 i <- which(nz) - 1L
150 : maechler 1575 }
151 :     new("lgTMatrix", i = i, j = i, x = x,
152 : maechler 1295 Dim = c(n,n), Dimnames = from@Dimnames) })
153 :    
154 :     setAs("ldiMatrix", "lgCMatrix",
155 :     function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))
156 :    
157 :    
158 : maechler 1447 if(FALSE) # now have faster "ddense" -> "dge"
159 : maechler 1174 setAs("ddiMatrix", "dgeMatrix",
160 :     function(from) as(as(from, "matrix"), "dgeMatrix"))
161 :    
162 : maechler 1109 setAs("matrix", "diagonalMatrix",
163 :     function(from) {
164 : maechler 1295 d <- dim(from)
165 : maechler 1109 if(d[1] != (n <- d[2])) stop("non-square matrix")
166 :     if(any(from[row(from) != col(from)] != 0))
167 :     stop("has non-zero off-diagonal entries")
168 : maechler 1295 x <- diag(from)
169 :     if(is.logical(x)) {
170 :     cl <- "ldiMatrix"
171 :     uni <- all(x)
172 :     } else {
173 :     cl <- "ddiMatrix"
174 :     uni <- all(x == 1)
175 :     storage.mode(x) <- "double"
176 : maechler 1575 } ## TODO: complex
177 : maechler 1295 new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
178 :     x = if(uni) x[FALSE] else x)
179 : maechler 1109 })
180 :    
181 :     ## ``generic'' coercion to diagonalMatrix : build on isDiagonal() and diag()
182 :     setAs("Matrix", "diagonalMatrix",
183 :     function(from) {
184 :     d <- dim(from)
185 :     if(d[1] != (n <- d[2])) stop("non-square matrix")
186 :     if(!isDiagonal(from)) stop("matrix is not diagonal")
187 :     ## else:
188 :     x <- diag(from)
189 :     if(is.logical(x)) {
190 :     cl <- "ldiMatrix"
191 :     uni <- all(x)
192 :     } else {
193 :     cl <- "ddiMatrix"
194 :     uni <- all(x == 1)
195 :     storage.mode(x) <- "double"
196 :     }
197 :     new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
198 :     x = if(uni) x[FALSE] else x)
199 :     })
200 :    
201 : maechler 1708
202 :     setMethod("diag", signature(x = "diagonalMatrix"),
203 : maechler 2052 function(x = 1, nrow, ncol) .diag.x(x))
204 : maechler 1708
205 : maechler 1799
206 : maechler 2098 subDiag <- function(x, i, j, ..., drop) {
207 : maechler 1799 x <- as(x, "sparseMatrix")
208 :     x <- if(missing(i))
209 :     x[, j, drop=drop]
210 :     else if(missing(j))
211 :     x[i, , drop=drop]
212 :     else
213 :     x[i,j, drop=drop]
214 : maechler 2120 if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
215 : maechler 1799 }
216 :    
217 :     setMethod("[", signature(x = "diagonalMatrix", i = "index",
218 :     j = "index", drop = "logical"), subDiag)
219 :     setMethod("[", signature(x = "diagonalMatrix", i = "index",
220 :     j = "missing", drop = "logical"),
221 : maechler 2098 function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))
222 : maechler 1799 setMethod("[", signature(x = "diagonalMatrix", i = "missing",
223 :     j = "index", drop = "logical"),
224 : maechler 2098 function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
225 : maechler 1799
226 : maechler 1617 ## When you assign to a diagonalMatrix, the result should be
227 : maechler 1708 ## diagonal or sparse ---
228 :     ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
229 : maechler 2098 ## Only(?) current bug: x[i] <- value is wrong when i is *vector*
230 :     replDiag <- function(x, i, j, ..., value) {
231 : maechler 1710 x <- as(x, "sparseMatrix")
232 :     if(missing(i))
233 :     x[, j] <- value
234 : maechler 2098 else if(missing(j)) { ## x[i , ] <- v *OR* x[i] <- v
235 :     na <- nargs()
236 :     ## message("diagnosing replDiag() -- nargs()= ", na)
237 :     if(na == 4)
238 :     x[i, ] <- value
239 :     else if(na == 3)
240 :     x[i] <- value
241 :     else stop("Internal bug: nargs()=",na,"; please report")
242 :     } else
243 : maechler 1710 x[i,j] <- value
244 :     if(isDiagonal(x)) as(x, "diagonalMatrix") else x
245 :     }
246 : maechler 1617
247 : maechler 1710 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
248 :     j = "index", value = "replValue"), replDiag)
249 : maechler 2098
250 : maechler 1710 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
251 :     j = "missing", value = "replValue"),
252 : maechler 2098 function(x,i,j, ..., value) {
253 :     ## message("before replDiag() -- nargs()= ", nargs())
254 :     if(nargs() == 3)
255 :     replDiag(x, i=i, value=value)
256 :     else ## nargs() == 4 :
257 :     replDiag(x, i=i, , value=value)
258 :     })
259 :    
260 : maechler 2096 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix
261 :     j = "missing", value = "replValue"),
262 : maechler 2098 function(x,i,j, ..., value) {
263 : maechler 2096 if(ncol(i) == 2) {
264 :     if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
265 :     x@x[ii] <- value
266 :     x
267 :     } else { ## no longer diagonal, but remain sparse:
268 :     x <- as(x, "sparseMatrix")
269 :     x[i] <- value
270 :     x
271 :     }
272 :     }
273 :     else { # behave as "base R": use as if vector
274 :     x <- as(x, "matrix")
275 :     x[i] <- value
276 :     Matrix(x)
277 :     }
278 :     })
279 :    
280 : maechler 1710 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
281 :     j = "index", value = "replValue"),
282 : maechler 2098 function(x,i,j, ..., value) replDiag(x, j=j, value=value))
283 : maechler 1710
284 :    
285 : maechler 1109 setMethod("t", signature(x = "diagonalMatrix"),
286 :     function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
287 :    
288 : maechler 1331 setMethod("isDiagonal", signature(object = "diagonalMatrix"),
289 :     function(object) TRUE)
290 :     setMethod("isTriangular", signature(object = "diagonalMatrix"),
291 :     function(object) TRUE)
292 : maechler 1109 setMethod("isSymmetric", signature(object = "diagonalMatrix"),
293 : maechler 2112 function(object, ...) TRUE)
294 : maechler 1109
295 : maechler 2112 setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
296 :     setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
297 :    
298 : maechler 1654 setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"
299 :     function(x, pivot) {
300 :     if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")
301 :     x@x <- sqrt(x@x)
302 :     x
303 :     })
304 :     ## chol(L) is L for logical diagonal:
305 :     setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)
306 :    
307 : maechler 1109 ## Basic Matrix Multiplication {many more to add}
308 : maechler 1654 ## ---------------------
309 :     ## Note that "ldi" logical are treated as numeric
310 : maechler 1109 diagdiagprod <- function(x, y) {
311 :     if(any(dim(x) != dim(y))) stop("non-matching dimensions")
312 :     if(x@diag != "U") {
313 : maechler 1654 if(y@diag != "U") {
314 :     nx <- x@x * y@x
315 :     if(is.numeric(nx) && !is.numeric(x@x))
316 :     x <- as(x, "dMatrix")
317 :     x@x <- as.numeric(nx)
318 :     }
319 :     return(x)
320 : maechler 1109 } else ## x is unit diagonal
321 :     return(y)
322 :     }
323 :    
324 : maechler 1654 setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
325 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix")
326 :    
327 : maechler 1654 formals(diagdiagprod) <- alist(x=, y=x)
328 :     setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
329 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix")
330 : maechler 1654 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
331 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix")
332 : maechler 1654 setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),
333 :     diagdiagprod, valueClass = "ddiMatrix")
334 :     setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),
335 :     diagdiagprod, valueClass = "ddiMatrix")
336 : maechler 1109
337 :    
338 :     diagmatprod <- function(x, y) {
339 :     dx <- dim(x)
340 :     dy <- dim(y)
341 :     if(dx[2] != dy[1]) stop("non-matching dimensions")
342 :     n <- dx[1]
343 :     as(if(x@diag == "U") y else x@x * y, "Matrix")
344 :     }
345 :    
346 :     setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
347 : maechler 1654 diagmatprod)
348 : maechler 1109 formals(diagmatprod) <- alist(x=, y=NULL)
349 :     setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
350 : maechler 1654 diagmatprod)
351 :     setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),
352 :     diagmatprod)
353 : maechler 1109
354 :     diagdgeprod <- function(x, y) {
355 :     dx <- dim(x)
356 :     dy <- dim(y)
357 :     if(dx[2] != dy[1]) stop("non-matching dimensions")
358 :     if(x@diag != "U")
359 :     y@x <- x@x * y@x
360 :     y
361 :     }
362 :     setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),
363 :     diagdgeprod, valueClass = "dgeMatrix")
364 :     formals(diagdgeprod) <- alist(x=, y=NULL)
365 :     setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
366 :     diagdgeprod, valueClass = "dgeMatrix")
367 :    
368 :     setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
369 :     function(x, y) {
370 : maechler 1635 dx <- dim(x)
371 :     dy <- dim(y)
372 :     if(dx[2] != dy[1]) stop("non-matching dimensions")
373 :     as(if(y@diag == "U") x else x * rep(y@x, each = dx[1]), "Matrix")
374 :     })
375 : maechler 1109
376 :     setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),
377 :     function(x, y) {
378 : maechler 1635 dx <- dim(x)
379 :     dy <- dim(y)
380 :     if(dx[2] != dy[1]) stop("non-matching dimensions")
381 :     if(y@diag == "N")
382 :     x@x <- x@x * rep(y@x, each = dx[1])
383 :     x
384 :     })
385 : maechler 1109
386 : maechler 1295 ## crossprod {more of these}
387 : maechler 1109
388 : maechler 1295 ## tcrossprod --- all are not yet there: do the dense ones here:
389 : maechler 1109
390 : maechler 1295 ## FIXME:
391 :     ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
392 :     ## function(x, y = NULL) {
393 :     ## })
394 : maechler 1109
395 : maechler 1295 ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),
396 :     ## function(x, y = NULL) {
397 :     ## })
398 : maechler 1109
399 : maechler 1799 setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
400 : maechler 2123 function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))
401 : maechler 1295
402 : maechler 1799 setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
403 : maechler 2123 function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))
404 : maechler 1799
405 :     setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
406 : maechler 2123 function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))
407 : maechler 1799
408 :     setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
409 : maechler 2123 function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))
410 : maechler 1799
411 :    
412 :     ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()
413 :     setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
414 :     function(x, y) as(x, "sparseMatrix") %*% y)
415 :     ## NB: The previous is *not* triggering for "ddi" o "dgC" (= distance 3)
416 :     ## since there's a "ddense" o "Csparse" at dist. 2 => triggers first.
417 :     ## ==> do this:
418 :     setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
419 :     function(x, y) as(x, "CsparseMatrix") %*% y)
420 :     setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
421 :     function(x, y) x %*% as(y, "CsparseMatrix"))
422 :     ## NB: this is *not* needed for Tsparse & Rsparse
423 :     ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
424 :     ## do indeed work by going through sparse (and *not* ddense)!
425 :    
426 :     setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
427 :     function(x, y) x %*% as(y, "sparseMatrix"))
428 :    
429 :    
430 :     setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
431 :     function(a, b, ...) {
432 :     a@x <- 1/ a@x
433 :     a@Dimnames <- a@Dimnames[2:1]
434 :     a
435 :     })
436 :    
437 :     solveDiag <- function(a, b, ...) {
438 :     if((n <- a@Dim[1]) != nrow(b))
439 :     stop("incompatible matrix dimensions")
440 :     ## trivially invert a 'in place' and multiply:
441 :     a@x <- 1/ a@x
442 :     a@Dimnames <- a@Dimnames[2:1]
443 :     a %*% b
444 :     }
445 :     setMethod("solve", signature(a = "diagonalMatrix", b = "matrix"),
446 :     solveDiag)
447 :     setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
448 :     solveDiag)
449 :    
450 : maechler 2106 ## Schur() ---> ./eigen.R
451 : maechler 1799
452 :    
453 :    
454 : maechler 1654 ### ---------------- diagonal o sparse -----------------------------
455 : maechler 1295
456 : maechler 1655
457 :     ## Use function for several signatures, in order to evade
458 :     ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
459 :     diagOdiag <- function(e1,e2) { # result should also be diagonal
460 :     r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
461 :     if(is.numeric(r)) {
462 :     if(is.numeric(e2@x)) {
463 :     e2@x <- r; return(.diag.2N(e2)) }
464 :     if(!is.numeric(e1@x))
465 :     ## e.g. e1, e2 are logical;
466 :     e1 <- as(e1, "dMatrix")
467 :     }
468 :     else if(is.logical(r))
469 :     e1 <- as(e1, "lMatrix")
470 :     else stop("intermediate 'r' is of type", typeof(r))
471 :     e1@x <- r
472 :     .diag.2N(e1)
473 :     }
474 :    
475 :     setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
476 :     diagOdiag)
477 :     ## These two are just for method disambiguation:
478 :     setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),
479 :     diagOdiag)
480 :     setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),
481 :     diagOdiag)
482 :    
483 : maechler 1845 ## FIXME: diagonal o triangular |--> triangular
484 :     ## ----- diagonal o symmetric |--> symmetric
485 :     ## {also when other is sparse: do these "here" --
486 :     ## before conversion to sparse, since that loses "diagonality"}
487 :    
488 : maechler 1655 ## For almost everything else, diag* shall be treated "as sparse" :
489 : maechler 1295 ## These are cheap implementations via coercion
490 :    
491 : maechler 1655 ## for disambiguation
492 : maechler 1654 setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),
493 :     function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))
494 :     setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),
495 :     function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))
496 : maechler 1655 ## in general:
497 :     setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),
498 :     function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))
499 :     setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),
500 :     function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))
501 : maechler 1654
502 : maechler 1655
503 :    
504 : maechler 1109 ## similar to prTriang() in ./Auxiliaries.R :
505 :     prDiag <-
506 :     function(x, digits = getOption("digits"), justify = "none", right = TRUE)
507 :     {
508 :     cf <- array(".", dim = x@Dim, dimnames = x@Dimnames)
509 :     cf[row(cf) == col(cf)] <-
510 :     sapply(diag(x), format, digits = digits, justify = justify)
511 :     print(cf, quote = FALSE, right = right)
512 :     invisible(x)
513 :     }
514 :    
515 :     setMethod("show", signature(object = "diagonalMatrix"),
516 : maechler 1592 function(object) {
517 :     d <- dim(object)
518 :     cl <- class(object)
519 :     cat(sprintf('%d x %d diagonal matrix of class "%s"\n',
520 :     d[1], d[2], cl))
521 :     prDiag(object)
522 :     })

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge