# SCM Repository

[matrix] Diff of /pkg/Matrix/R/dgCMatrix.R
 [matrix] / pkg / Matrix / R / dgCMatrix.R

# Diff of /pkg/Matrix/R/dgCMatrix.R

revision 2888, Wed Aug 7 13:56:43 2013 UTC revision 2889, Thu Aug 8 21:06:22 2013 UTC
# Line 79  Line 79
79  setMethod("lu", signature(x = "dgCMatrix"), LU.dgC)  setMethod("lu", signature(x = "dgCMatrix"), LU.dgC)
80
81  setMethod("lu", signature(x = "sparseMatrix"),  setMethod("lu", signature(x = "sparseMatrix"),
82            function(x, ...) # "FIXME": do in C, so can cache 'x@factors\$LU'            function(x, ...)
83            lu(as(as(as(x, "CsparseMatrix"), "dsparseMatrix"), "dgCMatrix"), ...))        {
84              ## "FIXME": do in C, so can cache 'x@factors\$LU'
85              warning(gettextf(
86                  "lu(<%s>) cannot cache the LU decomposition; consider working with \"dgCMatrix\"",
87                  class(x)),
88                      call. = FALSE, domain=NA)
89              lu(as(as(as(x, "CsparseMatrix"), "dsparseMatrix"), "dgCMatrix"), ...)
90    })
91
92
93    ## MM: see solveSparse() in  ~/R/MM/Pkg-ex/Matrix/Doran-A.R
94    .solve.sparse.dgC <- function(a, b) {
95        lu.a <- lu(a)
96        ## per default:  b = Identity = Diagonal(nrow(a)), however more efficiently
97        b.miss <- missing(b)
98        b.isMat <- b.miss || !is.null(dim(b))
99        if(b.miss) b <- .sparseDiagonal(dim(a)[1L])
100        ## bp := P %*% b
101        bp <- if(b.isMat) b[lu.a@p+1L, ] else b[lu.a@p+1L]
102        ## R:= U^{-1} L^{-1} P b
103        R <- solve(lu.a@U, solve(lu.a@L, bp))
104        ## result = Q'R = Q' U^{-1} L^{-1} P  b  = A^{-1} b,  as  A = P'LUQ
105        R[, invPerm(lu.a@q, zero.p=TRUE)]
106    }
107
108    ## FIXME: workaround, till  .Call(dgCMatrix_matrix_solve, a, b, sparse=TRUE)  works:
109    .solve.dgC <- function(a, b, sparse)
110        if(sparse) .solve.sparse.dgC(a, b) else .Call(dgCMatrix_matrix_solve, a, b, FALSE)
111
112    .solve.dgC.mat <- function(a, b, sparse=FALSE, ...) {
113        chk.s(...)
114        if(sparse) .solve.sparse.dgC(a, b) else .Call(dgCMatrix_matrix_solve, a, b, FALSE)
115    }
116
117  setMethod("solve", signature(a = "dgCMatrix", b = "matrix"),  setMethod("solve", signature(a = "dgCMatrix", b = "matrix"), .solve.dgC.mat)
function(a, b, ...) .Call(dgCMatrix_matrix_solve, a, b),
valueClass = "dgeMatrix")
118
119  setMethod("solve", signature(a = "dgCMatrix", b = "ddenseMatrix"),  setMethod("solve", signature(a = "dgCMatrix", b = "ddenseMatrix"), .solve.dgC.mat)
function(a, b, ...) .Call(dgCMatrix_matrix_solve, a, b),
valueClass = "dgeMatrix")
120
121  setMethod("solve", signature(a = "dgCMatrix", b = "dsparseMatrix"),  setMethod("solve", signature(a = "dgCMatrix", b = "dsparseMatrix"),
122            function(a, b, ...) {            function(a, b, sparse=FALSE, ...) { ## TODO: or rather TRUE [not back compatible] ??
123                  chk.s(...)
124                if(isSymmetric(a)) # TODO: fast cholmod_symmetric() for Cholesky                if(isSymmetric(a)) # TODO: fast cholmod_symmetric() for Cholesky
125                    solve(forceCspSymmetric(a, isTri=FALSE), b) #-> sparse result                    solve(forceCspSymmetric(a, isTri=FALSE), b) #-> sparse result
126                else .Call(dgCMatrix_matrix_solve, a, as(b, "denseMatrix"))                else ## FIXME: be better when sparse=TRUE (?)
127                      .solve.dgC(a, as(b, "denseMatrix"), sparse)
128            })            })
129
130  ## This is a really dumb method but some people apparently want it  ## This is a really dumb method but some people apparently want it
131  ## (MM: a bit less dumb now with possibility of staying sparse)  ## (MM: a bit less dumb now with possibility of staying sparse)
132  setMethod("solve", signature(a = "dgCMatrix", b = "missing"),  setMethod("solve", signature(a = "dgCMatrix", b = "missing"),
133            function(a, b, ...) {            function(a, b, sparse=FALSE, ...) { ## TODO: or rather TRUE [not back compatible] ??
134                  chk.s(...)
135                if(isSymmetric(a)) # TODO: fast cholmod_symmetric() for Cholesky                if(isSymmetric(a)) # TODO: fast cholmod_symmetric() for Cholesky
136                    solve(forceCspSymmetric(a, isTri=FALSE),                    solve(forceCspSymmetric(a, isTri=FALSE),
137                          b = Diagonal(nrow(a))) #-> sparse result                          b = Diagonal(nrow(a))) #-> sparse result
138                else .Call(dgCMatrix_matrix_solve, a, b = diag(nrow(a)))                else ## FIXME: be better when sparse=TRUE (?)
139                      ## .solve.dgC(a, b=diag(nrow(a)), sparse)
140                      if(sparse) .solve.sparse.dgC(a) # -> "smart" diagonal b
141                  else .Call(dgCMatrix_matrix_solve, a, b=diag(nrow(a)), FALSE)
142            })            })
143

Legend:
 Removed from v.2888 changed lines Added in v.2889