# R/methods_bdMatrix.R In davidruegamer/effortless: efficient operations on row-wise tensor product linked evaluations of special structures

```#' @rdname bdMatrix-class
setMethod("nrow", c("bdMatrix"),
function(x) sum( sapply(x@listOfBlocks, NROW)))

#' @rdname bdMatrix-class
setMethod("length", signature("bdMatrix"),
function(x) length(x@listOfBlocks))

#' @rdname bdMatrix-class
setMethod("ncol", c("bdMatrix"),
function(x) sum( sapply(x@listOfBlocks, NCOL)))

#' @rdname bdMatrix-class
setMethod("dim", c("bdMatrix"),
function(x) c(nrow(x),ncol(x)))

#' @rdname bdMatrix-class
setMethod("abs", c("bdMatrix"),
function(x) bdMatrix(lapply(x@listOfBlocks, abs)))

#' @rdname bdMatrix-class
setMethod("[", c("bdMatrix"),
function(x, i, j, ..., drop=FALSE)
{

indsplit <- getBlockIndices(x)

iMiss <- missing(i)
jMiss <- missing(j)

x <- lapply(1:length(x),
function(k){
ii <- if(!iMiss) which(indsplit\$rows[[k]] %in% i) else 1:length(indsplit\$rows[[k]])
jj <- if(!jMiss) which(indsplit\$cols[[k]] %in% j) else 1:length(indsplit\$cols[[k]])

if(length(ii)==0 | length(jj)==0) return(NULL) else
x[[k]][ii, jj, ..., drop=drop]
}
)

bdMatrix(x[!sapply(x, is.null)])

}
)

#' @rdname bdMatrix-class
setMethod("[[", c("bdMatrix"),
function(x, i, ...)
{

x@listOfBlocks[[i]]

}
)

#' @rdname bdMatrix-class
setMethod("*",
signature(e1="bdMatrix", e2="numeric"),
function(e1, e2){

e1@listOfBlocks <- mclapply(e1@listOfBlocks, function(l)l*e2)

return(e1)
}
)

#' @rdname bdMatrix-class
setMethod("crossprod", signature(x="bdMatrix"),
function(x) {

bdMatrix(mclapply(x@listOfBlocks, function(y) crossprod(y)))

} )

#' @rdname bdMatrix-class
setMethod("%*%", signature(x = "bdMatrix", y = "bdMatrix"),
function(x, y)
{

if( length(unique(c(sapply(x@listOfBlocks, NCOL), sapply(y@listOfBlocks, NROW)))) != 1 )
stop("Multiplication only implemented if all blocks have matching sizes.")
bdMatrix(mclapply(1:length(x), function(i) x[[i]] %*% y[[i]]))

})

#' @rdname bdMatrix-class
setMethod("%*%", signature(x = "bdMatrix", y = "numeric"),
function(x, y)
{

if(ncol(x)!=length(y))
stop("Dimension mismatch.")

blocks <- sapply(x@listOfBlocks, NROW)

sta <- c(1, cumsum(blocks)[-length(blocks)] + 1)
end <- cumsum(blocks)

unlist(mclapply(1:length(x), function(i) x[[i]] %*% y[sta[i]:end[i]]))

})

#' @rdname bdMatrix-class
setMethod("chol", signature(x="bdMatrix"),
function(x) {

bdMatrix(mclapply(x@listOfBlocks, function(y) chol(y)))

} )

#' @export
setGeneric("svd")
# setGeneric("svd", def = function(x, nu, nv) svd(x, nu, nv),
#            signature(x = "matrix", nu = "numeric", nv = "numeric"))

#' svd method for bdMatrix class
#' @export
#'
#' @param x a object of class \code{bdMatrix} whose SVD decomposition is to be computed
#' @param nu the number of left singular vectors to be computed. This must between \code{0} and \code{n = nrow(x)}.
#' @param nv the number of right singular vectors to be computed. This must be between \code{0} and \code{p = ncol(x)}.
#'
#' @export
#'
setMethod("svd", signature(x="bdMatrix"),
function(x, nu = min(nrow(x), p = ncol(x)), nv = min(nrow(x), p = ncol(x)))
{

res <- mclapply(x@listOfBlocks, function(y) svd(y, nu = nu, nv = nv))
list(d = unlist(lapply(res,"[[","d")),
u = if(nu!=0) bdiag(lapply(res,"[[","u")) else NULL,
v = if(nv!=0) bdiag(lapply(res,"[[","v")) else NULL)

}
)

# setMethod("svd", signature(x="bdMatrix", nu="numeric", nv="numeric"), svd.bdMatrix)
#
# svd.bdMatrix <- function(x) {
#
#   res <- mclapply([email protected], function(y) svd(y))
#   list(d = unlist(lapply(res,"[[","d")),
#        u = bdiag(lapply(res,"[[","u")),
#        v = bdiag(lapply(res,"[[","v")))
#
# }

# setMethod("svd", signature(x="bdMatrix"), svd.bdMatrix)

#' @rdname bdMatrix-class
setMethod("forceSymmetric", c("bdMatrix"),
function(x)
{

bdMatrix(mclapply(x, forceSymmetric))

}
)

#' @rdname bdMatrix-class
setMethod("+", signature(e1="bdMatrix", e2="dgCMatrix"),
function(e1, e2) {

stopifnot(all.equal(dim(e1),dim(e2))=="TRUE")

minDim <- min(sapply(e1@listOfBlocks, function(l) min(dim(l))))
if(any(as.numeric(triu(e2, k = minDim))!=0)){

warning("dgCMatrix with non-zero values outside the block range of bdMatrix.
Coercing bdMatrix to dgCMatrix.")

res <- bdiag(e1@listOfBlocks) + e2

}else{

indsplits <- getBlockIndices(e1)
res <- bdMatrix(mclapply(1:length(e1), function(k) e1[[k]] + e2[indsplits\$rows[[k]],
indsplits\$cols[[k]]]))

}

return(res)

} )

#' solve method for class bdMatrix
#'
#' @param a object of class \code{bdMatrix}
#' @details Due to the block diagonal structure of \code{a}, solving can be performed separately on block levels
#' if \code{a} only consists of quadratic blocks.
#' @import parallel
#' @export
setMethod("solve", c("bdMatrix"),
function(a) {

if(any(sapply(a, function(xx) diff(dim(xx))!=0))) stop("solve only implemented for quadratic blocks.")

res <- suppressWarnings(try(mclapply(a@listOfBlocks, solve)))
if(length(cres <- which(sapply(res,class)=="try-error"))>0) stop(res[cres]) else bdMatrix(res)

}
)

#' solve method for class bdMatrix
#'
#' @param a object of class \code{bdMatrix}
#' @param b object of class \code{bdMatrix}
#' @param ... further arguments passed to \code{solve} (see \code{?solve})
#'
#' @details Due to the block diagonal structure of \code{a}, solving can be performed separately on block levels
#' if \code{a} only consists of quadratic blocks. However, \code{solve} does not exploit the second argument \code{b}
#' as it is done in solve with vectors or matrices. Instead, \code{solve(a,b)} simply calculates \code{solve(a)\%*\%b}
#' blockwise for corresponding blocks.
#' @import parallel
#' @export
setMethod("solve", c("bdMatrix", "bdMatrix"),
function(a, b, ...) {

res <- try(bdMatrix(mclapply(1:length(a), function(i) solve(a@listOfBlocks[[i]],
as.vector(b@listOfBlocks[[i]]), ...))))
if(class(res)=="try-error") return(solve(a)%*%b) else return(res)

}
)

#' transpose method for bdMatrix objects
#'
#' @param x object of class \code{bdMatrix}
#'
#' @import parallel
#' @export
setMethod("t", c("bdMatrix"),
function(x) {

bdMatrix(mclapply(x@listOfBlocks, t))

}
)

#' @rdname bdMatrix-class
setMethod("max", c("bdMatrix"),
function(x)
{

max(sapply(x@listOfBlocks, max))

})
```
davidruegamer/effortless documentation built on May 15, 2019, 1:15 a.m.