File: triangularMatrix.R

package info (click to toggle)
rmatrix 1.3-2-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 7,024 kB
  • sloc: ansic: 42,435; makefile: 330; sh: 180
file content (69 lines) | stat: -rw-r--r-- 2,425 bytes parent folder | download
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
#### Methods for the virtual class 'triangularMatrix' of triangular matrices
#### Note that specific methods are in (8 different) ./?t?Matrix.R

setAs("triangularMatrix", "symmetricMatrix",
      function(from) as(as(from, "generalMatrix"), "symmetricMatrix"))

setAs("dgeMatrix", "triangularMatrix", function(from) asTri(from, "dtrMatrix"))
setAs("lgeMatrix", "triangularMatrix", function(from) asTri(from, "ltrMatrix"))
setAs("ngeMatrix", "triangularMatrix", function(from) asTri(from, "ntrMatrix"))

setAs("matrix", "triangularMatrix", function(from) mat2tri(from))

.tril.tr <- function(x, k = 0, ...) {  # are always square
    k <- as.integer(k[1])
    dd <- dim(x)
    stopifnot(-dd[1] <= k, k <= dd[1])  # had k <= 0
    if(k == 0 && x@uplo == "L") x
    else { ## more to do
        if(x@diag == "U") x <- .diagU2N(x, class(x), checkDense = TRUE)
        callNextMethod()
    }
}

.triu.tr <- function(x, k = 0, ...) {  # are always square
    k <- as.integer(k[1])
    dd <- dim(x)
    stopifnot(-dd[1] <= k, k <= dd[1])  # had k >= 0
    if(k == 0 && x@uplo == "U") x
    else { ## more to do
        if(x@diag == "U") x <- .diagU2N(x, class(x), checkDense = TRUE)
        callNextMethod()
    }
}

## In order to evade method dispatch ambiguity (with [CTR]sparse* and ddense*),
## but still remain "general"
## we use this hack instead of signature  x = "triangularMatrix" :

trCls <- names(getClass("triangularMatrix")@subclasses)
trCls. <- trCls[grep(".t.Matrix", trCls)]  # not "*Cholesky", "*Kaufman" ..
for(cls in trCls.) {
    setMethod("tril", cls, .tril.tr)
    setMethod("triu", cls, .triu.tr)
}

## ditto here:

isTriTri <- function(x, upper=NA) {
    if(is.na(upper)) structure(TRUE, kind=x@uplo)
    else if(upper) x@uplo == "U"
    else           x@uplo == "L"
}
for(cls in trCls)
    setMethod("isTriangular", signature(object = cls),
	      function(object, upper=NA, ...) isTriTri(object, upper))
## instead of just for ....   signature(object = "triangularMatrix")

rm(trCls, trCls., cls)

cholTrimat <- function(x, ...) {
    if(isDiagonal(x))
        cholDiag(as(x, "diagonalMatrix"), ...)
    else stop("'x' is not symmetric -- chol() undefined.")
}
setMethod("chol", signature(x = "dtCMatrix"), cholTrimat)
setMethod("chol", signature(x = "dtTMatrix"), cholTrimat)
setMethod("chol", signature(x = "dtRMatrix"), cholTrimat)
## setMethod("chol", signature(x = "triangularMatrix"), cholTrimat)