Revision 064f27d90fcbda7122d5b74508651e2917cc3feb authored by Doug and Martin on 17 October 2011, 00:00:00 UTC, committed by Gabor Csardi on 17 October 2011, 00:00:00 UTC
1 parent a707be9
Csparse.R
#### Methods for the virtual class 'CsparseMatrix' of sparse matrices stored in
#### "column compressed" format.
#### -- many more specific things are e.g. in ./dgCMatrix.R
setAs("CsparseMatrix", "TsparseMatrix",
function(from)
## |-> cholmod_C -> cholmod_T -> chm_triplet_to_SEXP
## modified to support triangular (../src/Csparse.c)
.Call(Csparse_to_Tsparse, from, is(from, "triangularMatrix")))
## special cases (when a specific "to" class is specified)
setAs("dgCMatrix", "dgTMatrix",
function(from) .Call(Csparse_to_Tsparse, from, FALSE))
setAs("dsCMatrix", "dsTMatrix",
function(from) .Call(Csparse_to_Tsparse, from, FALSE))
setAs("dsCMatrix", "dgCMatrix",
function(from) .Call(Csparse_symmetric_to_general, from))
for(prefix in c("d", "l", "n"))
setAs(paste(prefix,"sCMatrix",sep=''), "generalMatrix",
function(from) .Call(Csparse_symmetric_to_general, from))
setAs("dtCMatrix", "dtTMatrix",
function(from) .Call(Csparse_to_Tsparse, from, TRUE))
setAs("CsparseMatrix", "denseMatrix",
function(from) {
## |-> cholmod_C -> cholmod_dense -> chm_dense_to_dense
cld <- getClassDef(class(from))
if (extends(cld, "generalMatrix"))
.Call(Csparse_to_dense, from)
else {
## Csparse_to_dense loses symmetry and triangularity properties.
## With suitable changes to chm_dense_to_SEXP (../src/chm_common.c)
## we could do this in C code -- or do differently in C {FIXME!}
if (extends(cld, "triangularMatrix") && from@diag == "U")
from <- .Call(Csparse_diagU2N, from)
as(.Call(Csparse_to_dense, from), # -> "[dln]geMatrix"
paste(.M.kind(from, cld),
.dense.prefixes[.M.shape(from, cld)], "Matrix", sep=''))
}
})
## special cases (when a specific "to" class is specified)
setAs("dgCMatrix", "dgeMatrix",
function(from) .Call(Csparse_to_dense, from))
## cholmod_sparse_to_dense converts symmetric storage to general
## storage so symmetric classes are ok for conversion to matrix.
## unit triangular needs special handling
setAs("CsparseMatrix", "matrix",
function(from) {
## |-> cholmod_C -> cholmod_dense -> chm_dense_to_matrix
if (is(from, "triangularMatrix") && from@diag == "U")
from <- .Call(Csparse_diagU2N, from)
.Call(Csparse_to_matrix, from)
})
setAs("CsparseMatrix", "symmetricMatrix",
function(from) {
if(isSymmetric(from)) {
isTri <- is(from, "triangularMatrix")# i.e. effectively *diagonal*
if (isTri && from@diag == "U")
from <- .Call(Csparse_diagU2N, from)
.Call(Csparse_general_to_symmetric, from,
uplo = if(isTri) from@uplo else "U")
} else
stop("not a symmetric matrix; consider forceSymmetric() or symmpart()")
})
.validateCsparse <- function(x, sort.if.needed = FALSE)
.Call(Csparse_validate2, x, sort.if.needed)
##-> to be used in sparseMatrix(.), e.g. --- but is unused currently
## NB: 'sort.if.needed' is called 'maybe_modify' in C -- so be careful
### Some group methods:
setMethod("Math",
signature(x = "CsparseMatrix"),
function(x) {
f0 <- callGeneric(0.)
if(is0(f0)) {
## sparseness, symm., triang.,... preserved
cl <- class(x)
has.x <- !extends(cl, "nsparseMatrix")
## has.x <==> *not* nonzero-pattern == "nMatrix"
if(has.x) {
type <- storage.mode(x@x)
r <- callGeneric(x@x)
} else { ## nsparseMatrix
type <- ""
r <- rep.int(as.double(callGeneric(TRUE)),
switch(.sp.class(cl),
CsparseMatrix = length(x@i),
TsparseMatrix = length(x@i),
RsparseMatrix = length(x@j)))
}
if(type == storage.mode(r)) {
x@x <- r
x
} else { ## e.g. abs( <lgC> ) --> integer Csparse
## FIXME: when we have 'i*' classes, use them here:
rx <- new(sub("^.", "d", cl))
rx@x <- as.double(r)
## result is "same"
sNams <- slotNames(cl)
for(nm in sNams[sNams != "x"])
slot(rx, nm) <- slot(x, nm)
rx
}
} else { ## no sparseness
callGeneric(as_dense(x))
}
}) ## {Math}
### Subsetting -- basic things (drop = "missing") are done in ./Matrix.R
### ---------- "[" and (currently) also ./sparseMatrix.R
subCsp_cols <- function(x, j, drop)
{
## x[ , j, drop=drop] where we know that x is Csparse*
dn <- x@Dimnames
jj <- intI(j, n = x@Dim[2], dn[[2]], give.dn = FALSE)
r <- .Call(Csparse_submatrix, x, NULL, jj)
if(!is.null(n <- dn[[1]])) r@Dimnames[[1]] <- n
if(!is.null(n <- dn[[2]])) r@Dimnames[[2]] <- n[jj+1L]
if(drop && any(r@Dim == 1L)) drop(as(r, "matrix")) else r
}
subCsp_rows <- function(x, i, drop)# , cl = getClassDef(class(x))
{
## x[ i, drop=drop] where we know that x is Csparse*
dn <- x@Dimnames
ii <- intI(i, n = x@Dim[1], dn[[1]], give.dn = FALSE)
r <- .Call(Csparse_submatrix, x, ii, NULL)
if(!is.null(n <- dn[[1]])) r@Dimnames[[1]] <- n[ii+1L]
if(!is.null(n <- dn[[2]])) r@Dimnames[[2]] <- n
if(drop && any(r@Dim == 1L)) drop(as(r, "matrix")) else r
}
subCsp_ij <- function(x, i, j, drop)
{
## x[i, j, drop=drop] where we know that x is Csparse*
d <- x@Dim
dn <- x@Dimnames
## Take care that x[i,i] for symmetricM* stays symmetric
i.eq.j <- identical(i,j) # < want fast check
ii <- intI(i, n = d[1], dn[[1]], give.dn = FALSE)
jj <- if(i.eq.j && d[1] == d[2]) ii
else intI(j, n = d[2], dn[[2]], give.dn = FALSE)
r <- .Call(Csparse_submatrix, x, ii, jj)
if(!is.null(n <- dn[[1]])) r@Dimnames[[1]] <- n[ii + 1L]
if(!is.null(n <- dn[[2]])) r@Dimnames[[2]] <- n[jj + 1L]
if(!i.eq.j) {
if(drop && any(r@Dim == 1L)) drop(as(r, "matrix")) else r
} else { ## i == j
if(drop) drop <- any(r@Dim == 1L)
if(drop)
drop(as(r, "matrix"))
else if(extends((cx <- getClassDef(class(x))),
"symmetricMatrix")) ## TODO? make more efficient:
.gC2sym(r, uplo = x@uplo)## preserve uplo !
else if(extends(cx, "triangularMatrix") && !is.unsorted(ii))
as(r, "triangularMatrix")
else r
}
}
setMethod("[", signature(x = "CsparseMatrix", i = "index", j = "missing",
drop = "logical"),
function (x, i,j, ..., drop) {
na <- nargs()
Matrix.msg("Csp[i,m,l] : nargs()=",na, .M.level = 2)
if(na == 4)
subCsp_rows(x, i, drop=drop)
else if(na == 3)
.M.vectorSub(x, i) # as(x, "TsparseMatrix")[i, drop=drop]
else ## should not happen
stop("Matrix-internal error in <CsparseM>[i,,d]; please report")
})
setMethod("[", signature(x = "CsparseMatrix", i = "missing", j = "index",
drop = "logical"),
function (x,i,j, ..., drop) {
Matrix.msg("Csp[m,i,l] : nargs()=",nargs(), .M.level = 2)
subCsp_cols(x, j, drop=drop)
})
setMethod("[", signature(x = "CsparseMatrix",
i = "index", j = "index", drop = "logical"),
function (x, i, j, ..., drop) {
Matrix.msg("Csp[i,i,l] : nargs()=",nargs(), .M.level = 2)
subCsp_ij(x, i, j, drop=drop)
})
## workhorse for "[<-" -- both for d* and l* C-sparse matrices :
## --------- -----
replCmat <- function (x, i, j, ..., value)
{
di <- dim(x)
dn <- dimnames(x)
iMi <- missing(i)
jMi <- missing(j)
spV <- is(value, "sparseVector")
na <- nargs()
Matrix.msg("replCmat[x,i,j,.., val] : nargs()=", na,"; ",
if(iMi | jMi) sprintf("missing (i,j) = (%d,%d)", iMi,jMi),
.M.level = 2)
if(na == 3) { ## "vector (or 2-col) indexing" M[i] <- v
x <- as(x, "TsparseMatrix")
x[i] <- value # may change class e.g. from dtT* to dgT*
clx <- sub(".Matrix$", "CMatrix", class(x))
if(any(is0(x@x))) ## drop all values that "happen to be 0"
drop0(x, is.Csparse=FALSE) else as_CspClass(x, clx)
}
else ## nargs() == 4 :
replCmat4(x,
i1 = if(iMi) 0:(di[1] - 1L) else .ind.prep2(i, 1, di, dn),
i2 = if(jMi) 0:(di[2] - 1L) else .ind.prep2(j, 2, di, dn),
iMi=iMi, jMi=jMi, value=value, spV=spV)
} ## replCmat
replCmat4 <- function(x, i1, i2, iMi, jMi, value, spV = is(value,"sparseVector"))
{
dind <- c(length(i1), length(i2)) # dimension of replacement region
lenRepl <- prod(dind)
lenV <- length(value)
if(lenV == 0) {
if(lenRepl != 0)
stop("nothing to replace with")
else return(x)
}
## else: lenV := length(value) is > 0
if(lenRepl %% lenV != 0)
stop("number of items to replace is not a multiple of replacement length")
if(lenV > lenRepl)
stop("too many replacement values")
clx <- class(x)
clDx <- getClassDef(clx) # extends() , is() etc all use the class definition
## keep "symmetry" if changed here:
x.sym <- extends(clDx, "symmetricMatrix")
if(x.sym) { ## only half the indices are there..
## using array() for large dind is a disaster...
mkArray <- if(spV) # TODO: room for improvement
function(v, dim) spV2M(v, dim[1],dim[2]) else array
x.sym <-
(dind[1] == dind[2] && all(i1 == i2) &&
(lenRepl == 1 || lenV == 1 ||
isSymmetric(mkArray(value, dim=dind))))
## x.sym : result is *still* symmetric
x <- .Call(Csparse_symmetric_to_general, x) ## but do *not* redefine clx!
}
else if(extends(clDx, "triangularMatrix")) {
xU <- x@uplo == "U"
r.tri <- ((any(dind == 1) || dind[1] == dind[2]) &&
if(xU) max(i1) <= min(i2) else max(i2) <= min(i1))
if(r.tri) { ## result is *still* triangular
if(any(i1 == i2)) # diagonal will be changed
x <- diagU2N(x) # keeps class (!)
}
else { # go to "generalMatrix" and continue
x <- as(x, paste(.M.kind(x), "gCMatrix", sep='')) ## & do not redefine clx!
}
}
## Temporary hack for debugging --- remove eventually -- FIXME :
## see also MATRIX_SUBASSIGN_VERBOSE in ../src/t_Csparse_subassign.c
if(!is.null(v <- getOption("Matrix.subassign.verbose")) && v) {
op <- options(Matrix.verbose = 2); on.exit(options(op))
## the "hack" to signal "verbose" to the C code:
i1[1] <- -i1[1]
if(i1[1] == 0)
warning("i1[1] == 0 ==> C-level verbosity will not happen!")
}
if(extends(clDx, "dMatrix")) {
has.x <- TRUE
x <- .Call(dCsparse_subassign,
if(clx %in% c("dgCMatrix", "dtCMatrix"))x
else as(x, "dgCMatrix"),
i1, i2,
as(value, "sparseVector"))
}
else if(extends(clDx, "lMatrix")) {
has.x <- TRUE
x <- .Call(lCsparse_subassign,
if(clx %in% c("lgCMatrix", "ltCMatrix"))x
else as(x, "lgCMatrix"),
i1, i2,
as(value, "sparseVector"))
}
else if(extends(clDx, "nMatrix")) {
has.x <- FALSE
x <- .Call(nCsparse_subassign,
if(clx %in% c("ngCMatrix", "ntCMatrix"))x
else as(x, "ngCMatrix"),
i1, i2,
as(value, "sparseVector"))
}
else if(extends(clDx, "iMatrix")) {
has.x <- TRUE
x <- .Call(iCsparse_subassign,
if(clx %in% c("igCMatrix", "itCMatrix"))x
else as(x, "igCMatrix"),
i1, i2,
as(value, "sparseVector"))
}
else if(extends(clDx, "zMatrix")) {
has.x <- TRUE
x <- .Call(zCsparse_subassign,
if(clx %in% c("zgCMatrix", "ztCMatrix"))x
else as(x, "zgCMatrix"),
i1, i2,
## here we only want zsparseVector {to not have to do this in C}:
as(value, "zsparseVector"))
}
else { ## use "old" code ...
## does this happen ? ==>
if(identical(Sys.getenv("USER"),"maechler"))## does it still happen? __ FIXME __
stop("using \"old code\" part in Csparse subassignment")
## else
warning("using \"old code\" part in Csparse subassignment\n",
" >>> please report to Matrix-authors@r-project.org",
immediate. = TRUE, domain = NA)
xj <- .Call(Matrix_expand_pointers, x@p)
sel <- (!is.na(match(x@i, i1)) &
!is.na(match( xj, i2)))
has.x <- "x" %in% slotNames(clDx)# === slotNames(x),
## has.x <==> *not* nonzero-pattern == "nMatrix"
if(has.x && sum(sel) == lenRepl) { ## all entries to be replaced are non-zero:
## need indices instead of just 'sel', for, e.g., A[2:1, 2:1] <- v
non0 <- cbind(match(x@i[sel], i1),
match(xj [sel], i2)) - 1L
iN0 <- 1L + .Call(m_encodeInd, non0, di = dind, checkBounds = FALSE)
has0 <-
if(spV) length(value@i) < lenV else any(value[!is.na(value)] == 0)
if(lenV < lenRepl)
value <- rep(value, length = lenRepl)
## Ideally we only replace them where value != 0 and drop the value==0
## ones; FIXME: see Davis(2006) "2.7 Removing entries", p.16, e.g. use cs_dropzeros()
## but really could be faster and write something like cs_drop_k(A, k)
## v0 <- 0 == value
## if (lenRepl == 1) and v0 is TRUE, the following is not doing anything
##- --> ./dgTMatrix.R and its replTmat()
## x@x[sel[!v0]] <- value[!v0]
x@x[sel] <- as.vector(value[iN0])
if(extends(clDx, "compMatrix") && length(x@factors)) # drop cashed ones
x@factors <- list()
if(has0) x <- .Call(Csparse_drop, x, 0)
return(if(x.sym) as_CspClass(x, clx) else x)
}
## else go via Tsparse.. {FIXME: a waste! - we already have 'xj' ..}
## and inside Tsparse... the above i1, i2,..., sel are *all* redone!
## Happens too often {not anymore, I hope!}
##
Matrix.msg("wasteful C -> T -> C in replCmat(x,i,j,v) for <sparse>[i,j] <- v")
x <- as(x, "TsparseMatrix")
if(iMi)
x[ ,i2+1L] <- value
else if(jMi)
x[i1+1L, ] <- value
else
x[i1+1L,i2+1L] <- value
if(extends(clDx, "compMatrix") && length(x@factors)) # drop cashed ones
x@factors <- list()
}# else{ not using new memory-sparse code
if(has.x && any(is0(x@x))) ## drop all values that "happen to be 0"
as_CspClass(drop0(x), clx)
else as_CspClass(x, clx)
} ## replCmat4
setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "missing",
value = "replValue"),
replCmat)
setReplaceMethod("[", signature(x = "CsparseMatrix", i = "missing", j = "index",
value = "replValue"),
replCmat)
setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "index",
value = "replValue"),
replCmat)
### When the RHS 'value' is a sparseVector, now can use replCmat as well
setReplaceMethod("[", signature(x = "CsparseMatrix", i = "missing", j = "index",
value = "sparseVector"),
replCmat)
setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "missing",
value = "sparseVector"),
replCmat)
setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "index",
value = "sparseVector"),
replCmat)
## A[ ij ] <- value, where ij is (i,j) 2-column matrix
setReplaceMethod("[", signature(x = "CsparseMatrix", i = "matrix", j = "missing",
value = "replValue"),
function(x, i, value)
## goto Tsparse modify and convert back:
as(.TM.repl.i.mat(as(x, "TsparseMatrix"), i=i, value=value),
"CsparseMatrix"))
## more in ./sparseMatrix.R (and ./Matrix.R )
setMethod("t", signature(x = "CsparseMatrix"),
function(x) .Call(Csparse_transpose, x, is(x, "triangularMatrix")))
## NB: have extra tril(), triu() methods for symmetric ["dsC" and "lsC"] and
## NB: for all triangular ones, where the latter may 'callNextMethod()' these:
setMethod("tril", "CsparseMatrix",
function(x, k = 0, ...) {
k <- as.integer(k[1])
dd <- dim(x); sqr <- dd[1] == dd[2]
stopifnot(-dd[1] <= k, k <= dd[1]) # had k <= 0
r <- .Call(Csparse_band, x, -dd[1], k)
## return "lower triangular" if k <= 0
if(sqr && k <= 0)
as(r, paste(.M.kind(x), "tCMatrix", sep='')) else r
})
setMethod("triu", "CsparseMatrix",
function(x, k = 0, ...) {
k <- as.integer(k[1])
dd <- dim(x); sqr <- dd[1] == dd[2]
stopifnot(-dd[1] <= k, k <= dd[1]) # had k >= 0
r <- .Call(Csparse_band, x, k, dd[2])
## return "upper triangular" if k >= 0
if(sqr && k >= 0)
as(r, paste(.M.kind(x), "tCMatrix", sep='')) else r
})
setMethod("band", "CsparseMatrix",
function(x, k1, k2, ...) {
k1 <- as.integer(k1[1])
k2 <- as.integer(k2[1])
dd <- dim(x); sqr <- dd[1] == dd[2]
stopifnot(-dd[1] <= k1, k1 <= k2, k2 <= dd[2])
r <- .Call(Csparse_band, diagU2N(x), k1, k2)
if(sqr && k1 * k2 >= 0) ## triangular
as(r, paste(.M.kind(x), "tCMatrix", sep=''))
else if (k1 < 0 && k1 == -k2 && isSymmetric(x)) ## symmetric
as(r, paste(.M.kind(x), "sCMatrix", sep=''))
else
r
})
setMethod("diag", "CsparseMatrix",
function(x, nrow, ncol) {
## "FIXME": could be more efficient; creates new ..CMatrix:
dm <- .Call(Csparse_band, diagU2N(x), 0, 0)
dlen <- min(dm@Dim)
ind1 <- dm@i + 1L # 1-based index vector
if (is(dm, "nMatrix")) {
val <- rep.int(FALSE, dlen)
val[ind1] <- TRUE
}
else if (is(dm, "lMatrix")) {
val <- rep.int(FALSE, dlen)
val[ind1] <- as.logical(dm@x)
}
else {
val <- rep.int(0, dlen)
## cMatrix not yet active but for future expansion
if (is(dm, "cMatrix")) val <- as.complex(val)
val[ind1] <- dm@x
}
val
})
setMethod("writeMM", "CsparseMatrix",
function(obj, file, ...)
.Call(Csparse_MatrixMarket, obj, as.character(file)))
setMethod("Cholesky", signature(A = "CsparseMatrix"),
function(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, ...)
Cholesky(as(A, "symmetricMatrix"),
perm=perm, LDL=LDL, super=super, Imult=Imult, ...))
## TODO (in ../TODO for quite a while .....):
setMethod("Cholesky", signature(A = "nsparseMatrix"),
function(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, ...)
stop("Cholesky(<nsparse...>) -> *symbolic* factorization -- not yet implemented"))
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...