diff --git a/R/rowCounts.R b/R/rowCounts.R index 40e3cdab..d8830652 100644 --- a/R/rowCounts.R +++ b/R/rowCounts.R @@ -84,7 +84,7 @@ rowCounts <- function(x, rows=NULL, cols=NULL, value=TRUE, na.rm=FALSE, dim.=dim if (is.numeric(x) || is.logical(x)) { na.rm <- as.logical(na.rm) hasNAs <- TRUE - counts <- .Call("rowCounts", x, dim., rows, cols, value, 2L, na.rm, hasNAs, PACKAGE="matrixStats") + counts <- .Call("rowCounts", x, dim., rows, cols, value, 2L, na.rm, hasNAs, mc.cores, PACKAGE="matrixStats") } else { if (is.vector(x)) dim(x) <- dim. @@ -130,7 +130,7 @@ colCounts <- function(x, rows=NULL, cols=NULL, value=TRUE, na.rm=FALSE, dim.=dim if (is.numeric(x) || is.logical(x)) { na.rm <- as.logical(na.rm) hasNAs <- TRUE - counts <- .Call("colCounts", x, dim., rows, cols, value, 2L, na.rm, hasNAs, PACKAGE="matrixStats") + counts <- .Call("colCounts", x, dim., rows, cols, value, 2L, na.rm, hasNAs, mc.cores, PACKAGE="matrixStats") } else { if (is.vector(x)) dim(x) <- dim. @@ -193,7 +193,7 @@ rowAlls <- function(x, rows=NULL, cols=NULL, value=TRUE, na.rm=FALSE, dim.=dim(x if (is.numeric(x) || is.logical(x)) { na.rm <- as.logical(na.rm) hasNAs <- TRUE - counts <- .Call("rowCounts", x, dim., rows, cols, value, 0L, na.rm, hasNAs, PACKAGE="matrixStats") + counts <- .Call("rowCounts", x, dim., rows, cols, value, 0L, na.rm, hasNAs, mc.cores, PACKAGE="matrixStats") as.logical(counts) } else { if (is.vector(x)) dim(x) <- dim. @@ -216,7 +216,7 @@ colAlls <- function(x, rows=NULL, cols=NULL, value=TRUE, na.rm=FALSE, dim.=dim(x if (is.numeric(x) || is.logical(x)) { na.rm <- as.logical(na.rm) hasNAs <- TRUE - counts <- .Call("colCounts", x, dim., rows, cols, value, 0L, na.rm, hasNAs, PACKAGE="matrixStats") + counts <- .Call("colCounts", x, dim., rows, cols, value, 0L, na.rm, hasNAs, mc.cores, PACKAGE="matrixStats") as.logical(counts) } else { if (is.vector(x)) dim(x) <- dim. @@ -260,7 +260,7 @@ rowAnys <- function(x, rows=NULL, cols=NULL, value=TRUE, na.rm=FALSE, dim.=dim(x if (is.numeric(x) || is.logical(x)) { na.rm <- as.logical(na.rm) hasNAs <- TRUE - counts <- .Call("rowCounts", x, dim., rows, cols, value, 1L, na.rm, hasNAs, PACKAGE="matrixStats") + counts <- .Call("rowCounts", x, dim., rows, cols, value, 1L, na.rm, hasNAs, mc.cores, PACKAGE="matrixStats") as.logical(counts) } else { if (is.vector(x)) dim(x) <- dim. @@ -283,7 +283,7 @@ colAnys <- function(x, rows=NULL, cols=NULL, value=TRUE, na.rm=FALSE, dim.=dim(x if (is.numeric(x) || is.logical(x)) { na.rm <- as.logical(na.rm) hasNAs <- TRUE - counts <- .Call("colCounts", x, dim., rows, cols, value, 1L, na.rm, hasNAs, PACKAGE="matrixStats") + counts <- .Call("colCounts", x, dim., rows, cols, value, 1L, na.rm, hasNAs, mc.cores, PACKAGE="matrixStats") as.logical(counts) } else { if (is.vector(x)) dim(x) <- dim. diff --git a/R/rowCumsums.R b/R/rowCumsums.R index b45b3364..26472700 100644 --- a/R/rowCumsums.R +++ b/R/rowCumsums.R @@ -56,44 +56,44 @@ #*/########################################################################### rowCumsums <- function(x, rows=NULL, cols=NULL, dim.=dim(x), mc.cores=1L, ...) { dim <- as.integer(dim.); - .Call("rowCumsums", x, dim, rows, cols, TRUE, PACKAGE="matrixStats") + .Call("rowCumsums", x, dim, rows, cols, TRUE, mc.cores, PACKAGE="matrixStats") } colCumsums <- function(x, rows=NULL, cols=NULL, dim.=dim(x), mc.cores=1L, ...) { dim <- as.integer(dim.); - .Call("rowCumsums", x, dim, rows, cols, FALSE, PACKAGE="matrixStats") + .Call("rowCumsums", x, dim, rows, cols, FALSE, mc.cores, PACKAGE="matrixStats") } rowCumprods <- function(x, rows=NULL, cols=NULL, dim.=dim(x), mc.cores=1L, ...) { dim <- as.integer(dim.); - .Call("rowCumprods", x, dim, rows, cols, TRUE, PACKAGE="matrixStats") + .Call("rowCumprods", x, dim, rows, cols, TRUE, mc.cores, PACKAGE="matrixStats") } colCumprods <- function(x, rows=NULL, cols=NULL, dim.=dim(x), mc.cores=1L, ...) { dim <- as.integer(dim.); - .Call("rowCumprods", x, dim, rows, cols, FALSE, PACKAGE="matrixStats") + .Call("rowCumprods", x, dim, rows, cols, FALSE, mc.cores, PACKAGE="matrixStats") } rowCummins <- function(x, rows=NULL, cols=NULL, dim.=dim(x), mc.cores=1L, ...) { dim <- as.integer(dim.); - .Call("rowCummins", x, dim, rows, cols, TRUE, PACKAGE="matrixStats") + .Call("rowCummins", x, dim, rows, cols, TRUE, mc.cores, PACKAGE="matrixStats") } colCummins <- function(x, rows=NULL, cols=NULL, dim.=dim(x), mc.cores=1L, ...) { dim <- as.integer(dim.); - .Call("rowCummins", x, dim, rows, cols, FALSE, PACKAGE="matrixStats") + .Call("rowCummins", x, dim, rows, cols, FALSE, mc.cores, PACKAGE="matrixStats") } rowCummaxs <- function(x, rows=NULL, cols=NULL, dim.=dim(x), mc.cores=1L, ...) { dim <- as.integer(dim.); - .Call("rowCummaxs", x, dim, rows, cols, TRUE, PACKAGE="matrixStats") + .Call("rowCummaxs", x, dim, rows, cols, TRUE, mc.cores, PACKAGE="matrixStats") } colCummaxs <- function(x, rows=NULL, cols=NULL, dim.=dim(x), mc.cores=1L, ...) { dim <- as.integer(dim.); - .Call("rowCummaxs", x, dim, rows, cols, FALSE, PACKAGE="matrixStats") + .Call("rowCummaxs", x, dim, rows, cols, FALSE, mc.cores, PACKAGE="matrixStats") } diff --git a/R/rowDiffs.R b/R/rowDiffs.R index 595096ae..b01a0ac7 100644 --- a/R/rowDiffs.R +++ b/R/rowDiffs.R @@ -42,11 +42,11 @@ # @keyword univar #*/########################################################################### rowDiffs <- function(x, rows=NULL, cols=NULL, lag=1L, differences=1L, mc.cores=1L, ...) { - .Call("rowDiffs", x, dim(x), rows, cols, as.integer(lag), as.integer(differences), TRUE, PACKAGE="matrixStats") + .Call("rowDiffs", x, dim(x), rows, cols, as.integer(lag), as.integer(differences), TRUE, mc.cores, PACKAGE="matrixStats") } colDiffs <- function(x, rows=NULL, cols=NULL, lag=1L, differences=1L, mc.cores=1L, ...) { - .Call("rowDiffs", x, dim(x), rows, cols, as.integer(lag), as.integer(differences), FALSE, PACKAGE="matrixStats") + .Call("rowDiffs", x, dim(x), rows, cols, as.integer(lag), as.integer(differences), FALSE, mc.cores, PACKAGE="matrixStats") } diff --git a/R/rowLogSumExps.R b/R/rowLogSumExps.R index f93e37ce..ec9258d2 100644 --- a/R/rowLogSumExps.R +++ b/R/rowLogSumExps.R @@ -53,7 +53,7 @@ rowLogSumExps <- function(lx, rows=NULL, cols=NULL, na.rm=FALSE, dim.=dim(lx), m hasNA <- TRUE; res <- .Call("rowLogSumExps", lx, dim., rows, cols, - as.logical(na.rm), as.logical(hasNA), TRUE, + as.logical(na.rm), as.logical(hasNA), TRUE, mc.cores, PACKAGE="matrixStats"); # Preserve names @@ -71,7 +71,7 @@ colLogSumExps <- function(lx, rows=NULL, cols=NULL, na.rm=FALSE, dim.=dim(lx), m hasNA <- TRUE; res <- .Call("rowLogSumExps", lx, dim., rows, cols, - as.logical(na.rm), as.logical(hasNA), FALSE, + as.logical(na.rm), as.logical(hasNA), FALSE, mc.cores, PACKAGE="matrixStats"); # Preserve names diff --git a/R/rowMads.R b/R/rowMads.R index 040fcf6b..531223c8 100644 --- a/R/rowMads.R +++ b/R/rowMads.R @@ -12,7 +12,7 @@ rowMads <- function(x, rows=NULL, cols=NULL, center=NULL, constant=1.4826, na.rm na.rm <- as.logical(na.rm) constant = as.numeric(constant) hasNAs <- TRUE - x <- .Call("rowMads", x, dim., rows, cols, constant, na.rm, hasNAs, TRUE, PACKAGE="matrixStats") + x <- .Call("rowMads", x, dim., rows, cols, constant, na.rm, hasNAs, TRUE, mc.cores, PACKAGE="matrixStats") } else { # Apply subset on 'x' if (is.vector(x)) dim(x) <- dim. @@ -27,7 +27,7 @@ rowMads <- function(x, rows=NULL, cols=NULL, center=NULL, constant=1.4826, na.rm x <- x - center if (is.null(dim(x))) dim(x) <- dim. # prevent from dim dropping x <- abs(x) - x <- rowMedians(x, na.rm=na.rm, ...) + x <- rowMedians(x, na.rm=na.rm, mc.cores=mc.cores, ...) x <- constant*x } x @@ -48,7 +48,7 @@ colMads <- function(x, rows=NULL, cols=NULL, center=NULL, constant=1.4826, na.rm na.rm <- as.logical(na.rm) constant = as.numeric(constant) hasNAs <- TRUE - x <- .Call("rowMads", x, dim., rows, cols, constant, na.rm, hasNAs, FALSE, PACKAGE="matrixStats") + x <- .Call("rowMads", x, dim., rows, cols, constant, na.rm, hasNAs, FALSE, mc.cores, PACKAGE="matrixStats") } else { # Apply subset on 'x' if (is.vector(x)) dim(x) <- dim. @@ -67,7 +67,7 @@ colMads <- function(x, rows=NULL, cols=NULL, center=NULL, constant=1.4826, na.rm ## FAST: x <- t_tx_OP_y(x, center, OP="-", na.rm=FALSE) x <- abs(x) - x <- colMedians(x, na.rm=na.rm, ...) + x <- colMedians(x, na.rm=na.rm, mc.cores=mc.cores, ...) x <- constant*x } x diff --git a/R/rowMedians.S4.R b/R/rowMedians.S4.R index 561de8fe..1cefae0a 100644 --- a/R/rowMedians.S4.R +++ b/R/rowMedians.S4.R @@ -63,7 +63,7 @@ setMethod("rowMedians", signature(x="matrix"), function(x, rows=NULL, cols=NULL, na.rm <- as.logical(na.rm); hasNAs <- TRUE; # Add as an argument? /2007-08-24 - .Call("rowMedians", x, dim., rows, cols, na.rm, hasNAs, TRUE, PACKAGE="matrixStats"); + .Call("rowMedians", x, dim., rows, cols, na.rm, hasNAs, TRUE, mc.cores, PACKAGE="matrixStats"); }) @@ -76,7 +76,7 @@ setMethod("colMedians", signature(x="matrix"), function(x, rows=NULL, cols=NULL, na.rm <- as.logical(na.rm); hasNAs <- TRUE; # Add as an argument? /2007-08-24 - .Call("rowMedians", x, dim., rows, cols, na.rm, hasNAs, FALSE, PACKAGE="matrixStats"); + .Call("rowMedians", x, dim., rows, cols, na.rm, hasNAs, FALSE, mc.cores, PACKAGE="matrixStats"); }) diff --git a/R/rowOrderStats.R b/R/rowOrderStats.R index 37a33feb..6d4e8a40 100644 --- a/R/rowOrderStats.R +++ b/R/rowOrderStats.R @@ -67,7 +67,7 @@ rowOrderStats <- function(x, rows=NULL, cols=NULL, which, dim.=dim(x), mc.cores= } which <- as.integer(which) - .Call("rowOrderStats", x, dim., rows, cols, which, PACKAGE="matrixStats"); + .Call("rowOrderStats", x, dim., rows, cols, which, mc.cores, PACKAGE="matrixStats"); } @@ -80,7 +80,7 @@ colOrderStats <- function(x, rows=NULL, cols=NULL, which, dim.=dim(x), mc.cores= } which <- as.integer(which) - .Call("colOrderStats", x, dim., rows, cols, which, PACKAGE="matrixStats"); + .Call("colOrderStats", x, dim., rows, cols, which, mc.cores, PACKAGE="matrixStats"); } diff --git a/R/rowRanges.R b/R/rowRanges.R index 6bec18cd..1bbebdb5 100644 --- a/R/rowRanges.R +++ b/R/rowRanges.R @@ -57,38 +57,38 @@ rowRanges <- function(x, rows=NULL, cols=NULL, na.rm=FALSE, dim.=dim(x), mc.cores=1L, ...) { dim. <- as.integer(dim.) na.rm <- as.logical(na.rm) - .Call("rowRanges", x, dim., rows, cols, 2L, na.rm, TRUE, PACKAGE="matrixStats") + .Call("rowRanges", x, dim., rows, cols, 2L, na.rm, TRUE, mc.cores, PACKAGE="matrixStats") } rowMins <- function(x, rows=NULL, cols=NULL, na.rm=FALSE, dim.=dim(x), mc.cores=1L, ...) { dim. <- as.integer(dim.) na.rm <- as.logical(na.rm) - .Call("rowRanges", x, dim., rows, cols, 0L, na.rm, TRUE, PACKAGE="matrixStats") + .Call("rowRanges", x, dim., rows, cols, 0L, na.rm, TRUE, mc.cores, PACKAGE="matrixStats") } rowMaxs <- function(x, rows=NULL, cols=NULL, na.rm=FALSE, dim.=dim(x), mc.cores=1L, ...) { dim. <- as.integer(dim.) na.rm <- as.logical(na.rm) - .Call("rowRanges", x, dim., rows, cols, 1L, na.rm, TRUE, PACKAGE="matrixStats") + .Call("rowRanges", x, dim., rows, cols, 1L, na.rm, TRUE, mc.cores, PACKAGE="matrixStats") } colRanges <- function(x, rows=NULL, cols=NULL, na.rm=FALSE, dim.=dim(x), mc.cores=1L, ...) { dim. <- as.integer(dim.) na.rm <- as.logical(na.rm) - .Call("colRanges", x, dim., rows, cols, 2L, na.rm, TRUE, PACKAGE="matrixStats") + .Call("colRanges", x, dim., rows, cols, 2L, na.rm, TRUE, mc.cores, PACKAGE="matrixStats") } colMins <- function(x, rows=NULL, cols=NULL, na.rm=FALSE, dim.=dim(x), mc.cores=1L, ...) { dim. <- as.integer(dim.) na.rm <- as.logical(na.rm) - .Call("colRanges", x, dim., rows, cols, 0L, na.rm, TRUE, PACKAGE="matrixStats") + .Call("colRanges", x, dim., rows, cols, 0L, na.rm, TRUE, mc.cores, PACKAGE="matrixStats") } colMaxs <- function(x, rows=NULL, cols=NULL, na.rm=FALSE, dim.=dim(x), mc.cores=1L, ...) { dim. <- as.integer(dim.) na.rm <- as.logical(na.rm) - .Call("colRanges", x, dim., rows, cols, 1L, na.rm, TRUE, PACKAGE="matrixStats") + .Call("colRanges", x, dim., rows, cols, 1L, na.rm, TRUE, mc.cores, PACKAGE="matrixStats") } diff --git a/R/rowRanks.R b/R/rowRanks.R index fd0798d8..da9bb300 100644 --- a/R/rowRanks.R +++ b/R/rowRanks.R @@ -111,7 +111,7 @@ rowRanks <- function(x, rows=NULL, cols=NULL, ties.method=c("max", "average", "m dim. <- as.integer(dim.) # byrow=TRUE - .Call("rowRanksWithTies", x, dim., rows, cols, tiesMethod, TRUE, PACKAGE="matrixStats") + .Call("rowRanksWithTies", x, dim., rows, cols, tiesMethod, TRUE, mc.cores, PACKAGE="matrixStats") } @@ -133,7 +133,7 @@ colRanks <- function(x, rows=NULL, cols=NULL, ties.method=c("max", "average", "m dim. <- as.integer(dim.) # byrow=FALSE - y <- .Call("rowRanksWithTies", x, dim., rows, cols, tiesMethod, FALSE, PACKAGE="matrixStats") + y <- .Call("rowRanksWithTies", x, dim., rows, cols, tiesMethod, FALSE, mc.cores, PACKAGE="matrixStats") if (!preserveShape) y <- t(y) y } diff --git a/R/rowVars.R b/R/rowVars.R index c14c8a11..78d3cc1e 100644 --- a/R/rowVars.R +++ b/R/rowVars.R @@ -53,7 +53,7 @@ rowVars <- function(x, rows=NULL, cols=NULL, na.rm=FALSE, center=NULL, dim.=dim( if (is.null(center)) { na.rm <- as.logical(na.rm) hasNAs <- TRUE - sigma2 <- .Call("rowVars", x, dim., rows, cols, na.rm, hasNAs, TRUE, PACKAGE="matrixStats"); + sigma2 <- .Call("rowVars", x, dim., rows, cols, na.rm, hasNAs, TRUE, mc.cores, PACKAGE="matrixStats"); return(sigma2) } @@ -112,7 +112,7 @@ colVars <- function(x, rows=NULL, cols=NULL, na.rm=FALSE, center=NULL, dim.=dim( dim. <- as.integer(dim.) na.rm <- as.logical(na.rm) hasNAs <- TRUE - sigma2 <- .Call("rowVars", x, dim., rows, cols, na.rm, hasNAs, FALSE, PACKAGE="matrixStats"); + sigma2 <- .Call("rowVars", x, dim., rows, cols, na.rm, hasNAs, FALSE, mc.cores, PACKAGE="matrixStats"); return(sigma2) } diff --git a/inst/benchmarking/colRowCummins_parallel.md.rsp b/inst/benchmarking/colRowCummins_parallel.md.rsp index 8941fc5f..4f234c36 100644 --- a/inst/benchmarking/colRowCummins_parallel.md.rsp +++ b/inst/benchmarking/colRowCummins_parallel.md.rsp @@ -23,7 +23,7 @@ This report benchmark the performance of <%=colname%>() and <%=rowname%>() on mu ```r <%=withCapture({ <%@include file="R/random-matrices.R"%> -data <- rmatrices(mode=mode) +data <- rmatrices(mode=mode, scale=20) })%> ``` diff --git a/inst/benchmarking/colRowCumprods_parallel.md.rsp b/inst/benchmarking/colRowCumprods_parallel.md.rsp index 1aaf5de0..82b0f228 100644 --- a/inst/benchmarking/colRowCumprods_parallel.md.rsp +++ b/inst/benchmarking/colRowCumprods_parallel.md.rsp @@ -23,7 +23,7 @@ This report benchmark the performance of <%=colname%>() and <%=rowname%>() on mu ```r <%=withCapture({ <%@include file="R/random-matrices.R"%> -data <- rmatrices(mode=mode, range=c(-1,1)) +data <- rmatrices(mode=mode, range=c(-1,1), scale=20) })%> ``` diff --git a/inst/benchmarking/colRowCumsums_parallel.md.rsp b/inst/benchmarking/colRowCumsums_parallel.md.rsp index fa27afc1..ad59191d 100644 --- a/inst/benchmarking/colRowCumsums_parallel.md.rsp +++ b/inst/benchmarking/colRowCumsums_parallel.md.rsp @@ -23,7 +23,7 @@ This report benchmark the performance of <%=colname%>() and <%=rowname%>() on mu ```r <%=withCapture({ <%@include file="R/random-matrices.R"%> -data <- rmatrices(mode=mode) +data <- rmatrices(mode=mode, scale=20) })%> ``` diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 00000000..2f138d97 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,2 @@ +PKG_CFLAGS=-pthread +PKG_LIBS=-pthread diff --git a/src/colCounts.c b/src/colCounts.c index e072b696..f5540316 100644 --- a/src/colCounts.c +++ b/src/colCounts.c @@ -11,7 +11,7 @@ #define METHOD colCounts #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, X_C_TYPE value, int what, int narm, int hasna, int *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, X_C_TYPE value, int what, int narm, int hasna, int *ans, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -21,9 +21,9 @@ #include "templates-gen-matrix.h" -SEXP colCounts(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP value, SEXP what, SEXP naRm, SEXP hasNA) { +SEXP colCounts(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP value, SEXP what, SEXP naRm, SEXP hasNA, SEXP cores) { SEXP ans; - int narm, hasna, what2; + int narm, hasna, what2, cores2; R_xlen_t nrow, ncol; /* Argument 'x' and 'dim': */ @@ -40,6 +40,8 @@ SEXP colCounts(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP value, SEXP what, SE /* Argument 'what': */ what2 = asInteger(what); + if (what2 < 0 || what2 > 2) + error("INTERNAL ERROR: Unknown value of 'what' for rowCounts: %d", what2); /* Argument 'naRm': */ narm = asLogicalNoNA(naRm, "na.rm"); @@ -53,15 +55,24 @@ SEXP colCounts(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP value, SEXP what, SE void *crows = validateIndices(rows, nrow, 0, &nrows, &rowsType); void *ccols = validateIndices(cols, ncol, 0, &ncols, &colsType); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + /* R allocate an integer vector of length 'ncol' */ PROTECT(ans = allocVector(INTSXP, ncols)); if (isReal(x)) { - colCounts_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, asReal(value), what2, narm, hasna, INTEGER(ans)); + colCounts_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, asReal(value), what2, narm, hasna, INTEGER(ans), cores2); } else if (isInteger(x)) { - colCounts_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, asInteger(value), what2, narm, hasna, INTEGER(ans)); + colCounts_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, asInteger(value), what2, narm, hasna, INTEGER(ans), cores2); } else if (isLogical(x)) { - colCounts_Logical[rowsType][colsType](LOGICAL(x), nrow, ncol, crows, nrows, ccols, ncols, asLogical(value), what2, narm, hasna, INTEGER(ans)); + colCounts_Logical[rowsType][colsType](LOGICAL(x), nrow, ncol, crows, nrows, ccols, ncols, asLogical(value), what2, narm, hasna, INTEGER(ans), cores2); } UNPROTECT(1); @@ -105,11 +116,11 @@ SEXP count(SEXP x, SEXP idxs, SEXP value, SEXP what, SEXP naRm, SEXP hasNA) { PROTECT(ans = allocVector(INTSXP, 1)); if (isReal(x)) { - colCounts_Real[rowsType][colsType](REAL(x), nx, 1, crows, nrows, ccols, ncols, asReal(value), what2, narm, hasna, INTEGER(ans)); + colCounts_Real[rowsType][colsType](REAL(x), nx, 1, crows, nrows, ccols, ncols, asReal(value), what2, narm, hasna, INTEGER(ans), 1); } else if (isInteger(x)) { - colCounts_Integer[rowsType][colsType](INTEGER(x), nx, 1, crows, nrows, ccols, ncols, asInteger(value), what2, narm, hasna, INTEGER(ans)); + colCounts_Integer[rowsType][colsType](INTEGER(x), nx, 1, crows, nrows, ccols, ncols, asInteger(value), what2, narm, hasna, INTEGER(ans), 1); } else if (isLogical(x)) { - colCounts_Logical[rowsType][colsType](LOGICAL(x), nx, 1, crows, nrows, ccols, ncols, asLogical(value), what2, narm, hasna, INTEGER(ans)); + colCounts_Logical[rowsType][colsType](LOGICAL(x), nx, 1, crows, nrows, ccols, ncols, asLogical(value), what2, narm, hasna, INTEGER(ans), 1); } UNPROTECT(1); @@ -120,6 +131,8 @@ SEXP count(SEXP x, SEXP idxs, SEXP value, SEXP what, SEXP naRm, SEXP hasNA) { /*************************************************************************** HISTORY: + 2015-07-31 [DJ] + o Pthread processing. 2015-04-21 [DJ] o Supported subsetted computation. 2014-11-14 [HB] diff --git a/src/colCounts_TYPE-template.h b/src/colCounts_TYPE-template.h index 5e30237f..aae374b9 100644 --- a/src/colCounts_TYPE-template.h +++ b/src/colCounts_TYPE-template.h @@ -3,7 +3,7 @@ void colCounts_[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, X_C_TYPE value, int what, int narm, int hasna, int *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, X_C_TYPE value, int what, int narm, int hasna, int *ans, int cores Arguments: The following macros ("arguments") should be defined for the @@ -14,6 +14,8 @@ Copyright: Henrik Bengtsson, 2014 ***********************************************************************/ +#include +#include #include "types.h" /* Expand arguments: @@ -22,7 +24,102 @@ #include "templates-types.h" +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, X_C_TYPE, value); + POP_ARGUMENT(buffer0, ii, int, what); + POP_ARGUMENT(buffer0, ii, int, narm); + POP_ARGUMENT(buffer0, ii, int, hasna); + POP_ARGUMENT(buffer0, ii, int*, ans); + + ans += begin; + + ncols = end - begin; + int colsType; +#ifdef COLS_TYPE + cols = (COLS_C_TYPE*) cols + begin; + colsType = COLS_TYPE_CODE; +#else + cols = indicesFromRange(begin, end, &colsType); +#endif + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + METHOD_NAME[ROWS_TYPE_CODE][colsType](x, nrow, ncol, rows, nrows, cols, ncols, value, what, narm, hasna, ans, 0); + +#ifndef COLS_TYPE + Free(cols); +#endif + return NULL; +} + + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + if (cores > 1 && ncols > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(value) + sizeof(what) + sizeof(narm) + sizeof(hasna) + sizeof(ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, ncols); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, value); + PUSH_ARGUMENT(buffer, ii, what); + PUSH_ARGUMENT(buffer, ii, narm); + PUSH_ARGUMENT(buffer, ii, hasna); + PUSH_ARGUMENT(buffer, ii, ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (ncols + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < ncols) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, ncols); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj; R_xlen_t colBegin, idx; int count; @@ -150,14 +247,14 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { ans[jj] = count; } /* for (jj ...) */ } /* if (X_ISNAN(value)) */ - } else { - error("INTERNAL ERROR: Unknown value of 'what' for colCounts: %d", what); } /* if (what) */ } /*************************************************************************** HISTORY: + 2015-07-31 [DJ] + o Pthread processing. 2015-04-18 [DJ] o Supported subsetted computation. 2014-11-14 [HB] diff --git a/src/colOrderStats.c b/src/colOrderStats.c index 6f7bbfa8..2bcdefce 100644 --- a/src/colOrderStats.c +++ b/src/colOrderStats.c @@ -1,6 +1,6 @@ /*************************************************************************** Public methods: - SEXP colOrderStats(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP which) + SEXP colOrderStats(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP which, SEXP cores) Authors: Henrik Bengtsson @@ -14,7 +14,7 @@ #define METHOD colOrderStats #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, R_xlen_t qq, X_C_TYPE *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, R_xlen_t qq, X_C_TYPE *ans, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -22,9 +22,10 @@ #include "templates-gen-matrix.h" -SEXP colOrderStats(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP which) { +SEXP colOrderStats(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP which, SEXP cores) { SEXP ans = NILSXP; R_xlen_t nrow, ncol, qq; + int cores2; /* Argument 'x' and 'dim': */ assertArgMatrix(x, dim, (R_TYPE_INT | R_TYPE_REAL), "x"); @@ -41,8 +42,18 @@ SEXP colOrderStats(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP which) { /* Argument 'rows' and 'cols': */ R_xlen_t nrows, ncols; int rowsType, colsType; - void *crows = validateIndices(rows, nrow, 0, &nrows, &rowsType); - void *ccols = validateIndices(cols, ncol, 0, &ncols, &colsType); + int rowsHasna, colsHasna; + void *crows = validateIndicesCheckNA(rows, nrow, 0, &nrows, &rowsType, &rowsHasna); + void *ccols = validateIndicesCheckNA(cols, ncol, 0, &ncols, &colsType, &colsHasna); + + // Check missing rows + if (rowsHasna && ncols > 0) { + error("Argument 'rows' must not contain missing value"); + } + // Check missing cols + if (colsHasna && nrows > 0) { + error("Argument 'cols' must not contain missing value"); + } /* Subtract one here, since rPsort does zero based addressing */ qq = asInteger(which) - 1; @@ -52,14 +63,23 @@ SEXP colOrderStats(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP which) { error("Argument 'which' is out of range."); } +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + /* Double matrices are more common to use. */ if (isReal(x)) { PROTECT(ans = allocVector(REALSXP, ncols)); - colOrderStats_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, qq, REAL(ans)); + colOrderStats_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, qq, REAL(ans), cores2); UNPROTECT(1); } else if (isInteger(x)) { PROTECT(ans = allocVector(INTSXP, ncols)); - colOrderStats_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, qq, INTEGER(ans)); + colOrderStats_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, qq, INTEGER(ans), cores2); UNPROTECT(1); } @@ -69,6 +89,8 @@ SEXP colOrderStats(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP which) { /*************************************************************************** HISTORY: + 2015-08-12 [DJ] + o Pthread processing. 2015-07-08 [DJ] o Supported subsetted computation. 2014-11-16 [HB] diff --git a/src/colOrderStats_TYPE-template.h b/src/colOrderStats_TYPE-template.h index 7ccd7da5..7bac38ba 100644 --- a/src/colOrderStats_TYPE-template.h +++ b/src/colOrderStats_TYPE-template.h @@ -3,7 +3,7 @@ void colOrderStats_[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, R_xlen_t qq, X_C_TYPE *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, R_xlen_t qq, X_C_TYPE *ans, int cores Arguments: The following macros ("arguments") should be defined for the @@ -20,6 +20,8 @@ ***********************************************************************/ #include #include +#include +#include #include "types.h" /* Expand arguments: @@ -29,35 +31,109 @@ #include "templates-types.h" +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, qq); + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, ans); + + ans += begin; + + ncols = end - begin; + int colsType; +#ifdef COLS_TYPE + cols = (COLS_C_TYPE*) cols + begin; + colsType = COLS_TYPE_CODE; +#else + cols = indicesFromRange(begin, end, &colsType); +#endif + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + METHOD_NAME[ROWS_TYPE_CODE][colsType](x, nrow, ncol, rows, nrows, cols, ncols, qq, ans, 0); + +#ifndef COLS_TYPE + Free(cols); +#endif + return NULL; +} + + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + if (cores > 1 && ncols > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(qq) + sizeof(ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, ncols); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, qq); + PUSH_ARGUMENT(buffer, ii, ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (ncols + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < ncols) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, ncols); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj; R_xlen_t offset; X_C_TYPE *values; #ifdef ROWS_TYPE ROWS_C_TYPE *crows = (ROWS_C_TYPE*) rows; - // Check missing rows - for (ii=0; ii < nrows; ++ii) { - if (ROW_INDEX(crows,ii) == NA_R_XLEN_T) break; - } - if (ii < nrows && ncols > 0) { - error("Argument 'rows' must not contain missing value"); - } #endif #ifdef COLS_TYPE COLS_C_TYPE *ccols = (COLS_C_TYPE*) cols; - // Check missing cols - for (jj=0; jj < ncols; ++jj) { - if (COL_INDEX(ccols,jj) == NA_R_XLEN_T) break; - } - if (jj < ncols && nrows > 0) { - error("Argument 'cols' must not contain missing value"); - } #endif - /* R allocate memory for the 'values'. This will be - taken care of by the R garbage collector later on. */ - values = (X_C_TYPE *) R_alloc(nrows, sizeof(X_C_TYPE)); + /* C allocate memory for the 'values'. Should be freed manually. */ + values = Calloc(nrows, X_C_TYPE); for (jj=0; jj < ncols; jj++) { offset = COL_INDEX_NONA(ccols,jj) * nrow; @@ -72,11 +148,15 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { ans[jj] = values[qq]; } + + Free(values); } /*************************************************************************** HISTORY: + 2015-08-12 [DJ] + o Pthread processing. 2015-07-08 [DJ] o Supported subsetted computation. 2014-11-16 [HB] diff --git a/src/colRanges.c b/src/colRanges.c index bcc45ec8..ca26650d 100644 --- a/src/colRanges.c +++ b/src/colRanges.c @@ -12,7 +12,7 @@ #define METHOD colRanges #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int what, int narm, int hasna, X_C_TYPE *ans, int *is_counted +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int what, int narm, int hasna, X_C_TYPE *ans, R_xlen_t nans, int *is_counted, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -20,12 +20,12 @@ #include "templates-gen-matrix.h" -SEXP colRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEXP hasNA) { +SEXP colRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEXP hasNA, SEXP cores) { SEXP ans = NILSXP, ans2 = NILSXP; int *mins, *maxs; double *mins2, *maxs2; int *is_counted, all_counted = 0; - int what2, narm, hasna; + int what2, narm, hasna, cores2; R_xlen_t nrow, ncol, jj; /* Argument 'x' and 'dim': */ @@ -54,6 +54,15 @@ SEXP colRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEX void *crows = validateIndices(rows, nrow, 0, &nrows, &rowsType); void *ccols = validateIndices(cols, ncol, 0, &ncols, &colsType); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + is_counted = (int *) R_alloc(ncols, sizeof(int)); if (isReal(x)) { @@ -62,7 +71,7 @@ SEXP colRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEX } else { PROTECT(ans = allocVector(REALSXP, ncols)); } - colRanges_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, what2, narm, hasna, REAL(ans), is_counted); + colRanges_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, what2, narm, hasna, REAL(ans), ncols, is_counted, cores2); UNPROTECT(1); } else if (isInteger(x)) { if (what2 == 2) { @@ -70,7 +79,7 @@ SEXP colRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEX } else { PROTECT(ans = allocVector(INTSXP, ncols)); } - colRanges_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, what2, narm, hasna, INTEGER(ans), is_counted); + colRanges_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, what2, narm, hasna, INTEGER(ans), ncols, is_counted, cores2); /* Any entries with zero non-missing values? */ all_counted = 1; @@ -138,6 +147,8 @@ SEXP colRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEX /*************************************************************************** HISTORY: + 2015-08-12 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-16 [HB] diff --git a/src/colRanges_TYPE-template.h b/src/colRanges_TYPE-template.h index 72e97887..27721e23 100644 --- a/src/colRanges_TYPE-template.h +++ b/src/colRanges_TYPE-template.h @@ -3,10 +3,10 @@ void colRanges_[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int what, int narm, int hasna, X_C_TYPE *ans, int *is_counted + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int what, int narm, int hasna, X_C_TYPE *ans, R_xlen_t nans, int *is_counted, int cores Arguments: - The following macros ("arguments") should be defined for the + The following macros ("arguments") should be defined for the template to work as intended. - METHOD_NAME: the name of the resulting function @@ -17,18 +17,118 @@ Henrik Bengtsson. Copyright: Henrik Bengtsson, 2014 - ***********************************************************************/ + ***********************************************************************/ #include +#include +#include #include "types.h" /* Expand arguments: X_TYPE => (X_C_TYPE, X_IN_C, [METHOD_NAME]) ANS_TYPE => (ANS_SXP, ANS_NA, ANS_C_TYPE, ANS_IN_C) */ -#include "templates-types.h" +#include "templates-types.h" + + +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, int, what); + POP_ARGUMENT(buffer0, ii, int, narm); + POP_ARGUMENT(buffer0, ii, int, hasna); + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, ans); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nans); + POP_ARGUMENT(buffer0, ii, int*, is_counted); + + ans += begin; + is_counted += begin; + + ncols = end - begin; + int colsType; +#ifdef COLS_TYPE + cols = (COLS_C_TYPE*) cols + begin; + colsType = COLS_TYPE_CODE; +#else + cols = indicesFromRange(begin, end, &colsType); +#endif + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + METHOD_NAME[ROWS_TYPE_CODE][colsType](x, nrow, ncol, rows, nrows, cols, ncols, what, narm, hasna, ans, nans, is_counted, 0); + +#ifndef COLS_TYPE + Free(cols); +#endif + return NULL; +} RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + if (cores > 1 && ncols > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(what) + sizeof(narm) + sizeof(hasna) + sizeof(ans) + sizeof(nans) + sizeof(is_counted); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, ncols); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, what); + PUSH_ARGUMENT(buffer, ii, narm); + PUSH_ARGUMENT(buffer, ii, hasna); + PUSH_ARGUMENT(buffer, ii, ans); + PUSH_ARGUMENT(buffer, ii, nans); + PUSH_ARGUMENT(buffer, ii, is_counted); + + pthread_t threads[cores]; + + R_xlen_t gap = (ncols + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < ncols) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, ncols); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj; R_xlen_t colBegin, idx; X_C_TYPE value, *mins = NULL, *maxs = NULL; @@ -92,7 +192,7 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } else if (what == 1) { /* colMaxs() */ maxs = ans; - + for (jj=0; jj < ncols; jj++) { colBegin = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); @@ -131,8 +231,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } else if (what == 2) { /* colRanges() */ mins = ans; - maxs = &ans[ncols]; - + maxs = ans + nans; + for (jj=0; jj < ncols; jj++) { colBegin = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); @@ -179,12 +279,12 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { if (what == 0) { /* colMins() */ mins = ans; - + /* Initiate results */ for (jj=0; jj < ncols; jj++) { mins[jj] = x[jj]; } - + for (jj=1; jj < ncols; jj++) { colBegin = COL_INDEX_NONA(ccols,jj) * nrow; for (ii=0; ii < nrows; ii++) { @@ -195,12 +295,12 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } else if (what == 1) { /* colMax() */ maxs = ans; - + /* Initiate results */ for (jj=0; jj < ncols; jj++) { maxs[jj] = x[jj]; } - + for (jj=1; jj < ncols; jj++) { colBegin = COL_INDEX_NONA(ccols,jj) * nrow; for (ii=0; ii < nrows; ii++) { @@ -211,14 +311,14 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } else if (what == 2) { /* colRanges()*/ mins = ans; - maxs = &ans[ncols]; - + maxs = ans + nans; + /* Initiate results */ for (jj=0; jj < ncols; jj++) { mins[jj] = x[jj]; maxs[jj] = x[jj]; } - + for (jj=1; jj < ncols; jj++) { colBegin = COL_INDEX_NONA(ccols,jj) * nrow; for (ii=0; ii < nrows; ii++) { @@ -237,6 +337,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { /*************************************************************************** HISTORY: + 2015-08-12 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-16 [HB] diff --git a/src/logSumExp_TYPE-template.h b/src/logSumExp_TYPE-template.h index c45dd085..ee17f810 100644 --- a/src/logSumExp_TYPE-template.h +++ b/src/logSumExp_TYPE-template.h @@ -98,7 +98,8 @@ RETURN_TYPE METHOD_NAME_IDXS(ARGUMENTS_LIST) { xMaxIsNA = ISNAN(xMax); } - R_CHECK_USER_INTERRUPT(ii); + // TODO: interrupt under pthread + // R_CHECK_USER_INTERRUPT(ii); } /* for (ii ...) */ } else { for (ii=1; ii < nidxs; ii++) { @@ -120,7 +121,8 @@ RETURN_TYPE METHOD_NAME_IDXS(ARGUMENTS_LIST) { xMaxIsNA = ISNAN(xMax); } - R_CHECK_USER_INTERRUPT(ii); + // TODO: interrupt under pthread + // R_CHECK_USER_INTERRUPT(ii); } /* for (ii ...) */ } /* by */ @@ -152,7 +154,8 @@ RETURN_TYPE METHOD_NAME_IDXS(ARGUMENTS_LIST) { /* Early LDOUBLE stopping on -Inf/+Inf and user interrupt? */ if (ii % 1048576 == 0) { if (!R_FINITE(sum)) break; - R_CheckUserInterrupt(); + // TODO: interrupt under pthread + // R_CheckUserInterrupt(); } } /* for (ii ...) */ } else { @@ -171,7 +174,8 @@ RETURN_TYPE METHOD_NAME_IDXS(ARGUMENTS_LIST) { /* Early LDOUBLE stopping on -Inf/+Inf and user interrupt? */ if (ii % 1048576 == 0) { if (!R_FINITE(sum)) break; - R_CheckUserInterrupt(); + // TODO: interrupt under pthread + // R_CheckUserInterrupt(); } } /* for (ii ...) */ } /* if (by) */ diff --git a/src/rowCounts.c b/src/rowCounts.c index fd015b2c..645c1924 100644 --- a/src/rowCounts.c +++ b/src/rowCounts.c @@ -11,7 +11,7 @@ #define METHOD rowCounts #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, X_C_TYPE value, int what, int narm, int hasna, int *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, X_C_TYPE value, int what, int narm, int hasna, int *ans, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -21,9 +21,9 @@ #include "templates-gen-matrix.h" -SEXP rowCounts(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP value, SEXP what, SEXP naRm, SEXP hasNA) { +SEXP rowCounts(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP value, SEXP what, SEXP naRm, SEXP hasNA, SEXP cores) { SEXP ans; - int narm, hasna, what2; + int narm, hasna, what2, cores2; R_xlen_t nrow, ncol; /* Argument 'x' & 'dim': */ @@ -40,6 +40,8 @@ SEXP rowCounts(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP value, SEXP what, SE /* Argument 'what': */ what2 = asInteger(what); + if (what2 < 0 || what2 > 2) + error("INTERNAL ERROR: Unknown value of 'what' for rowCounts: %d", what2); /* Argument 'naRm': */ narm = asLogicalNoNA(naRm, "na.rm"); @@ -53,16 +55,25 @@ SEXP rowCounts(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP value, SEXP what, SE void *crows = validateIndices(rows, nrow, 0, &nrows, &rowsType); void *ccols = validateIndices(cols, ncol, 0, &ncols, &colsType); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + /* R allocate a double vector of length 'nrow' */ PROTECT(ans = allocVector(INTSXP, nrows)); /* Double matrices are more common to use. */ if (isReal(x)) { - rowCounts_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, asReal(value), what2, narm, hasna, INTEGER(ans)); + rowCounts_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, asReal(value), what2, narm, hasna, INTEGER(ans), cores2); } else if (isInteger(x)) { - rowCounts_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, asInteger(value), what2, narm, hasna, INTEGER(ans)); + rowCounts_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, asInteger(value), what2, narm, hasna, INTEGER(ans), cores2); } else if (isLogical(x)) { - rowCounts_Logical[rowsType][colsType](LOGICAL(x), nrow, ncol, crows, nrows, ccols, ncols, asLogical(value), what2, narm, hasna, INTEGER(ans)); + rowCounts_Logical[rowsType][colsType](LOGICAL(x), nrow, ncol, crows, nrows, ccols, ncols, asLogical(value), what2, narm, hasna, INTEGER(ans), cores2); } UNPROTECT(1); @@ -73,6 +84,8 @@ SEXP rowCounts(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP value, SEXP what, SE /*************************************************************************** HISTORY: + 2015-07-30 [DJ] + o Pthread processing. 2015-04-13 [DJ] o Supported subsetted computation. 2014-06-02 [HB] diff --git a/src/rowCounts_TYPE-template.h b/src/rowCounts_TYPE-template.h index c7035630..c6cf44ef 100644 --- a/src/rowCounts_TYPE-template.h +++ b/src/rowCounts_TYPE-template.h @@ -3,26 +3,123 @@ void rowCounts_[ROWS_TYPE][COLS_TYPE](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, X_C_TYPE value, int what, int narm, int hasna, int *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, X_C_TYPE value, int what, int narm, int hasna, int *ans, int cores Arguments: - The following macros ("arguments") should be defined for the + The following macros ("arguments") should be defined for the template to work as intended. - METHOD_NAME: the name of the resulting function - X_TYPE: 'i', 'r', or 'l' Copyright: Henrik Bengtsson, 2014 - ***********************************************************************/ + ***********************************************************************/ +#include +#include #include "types.h" /* Expand arguments: X_TYPE => (X_C_TYPE, X_IN_C, [METHOD_NAME]) */ -#include "templates-types.h" +#include "templates-types.h" + + +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, X_C_TYPE, value); + POP_ARGUMENT(buffer0, ii, int, what); + POP_ARGUMENT(buffer0, ii, int, narm); + POP_ARGUMENT(buffer0, ii, int, hasna); + POP_ARGUMENT(buffer0, ii, int*, ans); + + ans += begin; + + nrows = end - begin; + int rowsType; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; + rowsType = ROWS_TYPE_CODE; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + METHOD_NAME[rowsType][COLS_TYPE_CODE](x, nrow, ncol, rows, nrows, cols, ncols, value, what, narm, hasna, ans, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + return NULL; +} RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + if (cores > 1 && nrows > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(value) + sizeof(what) + sizeof(narm) + sizeof(hasna) + sizeof(ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nrows); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, value); + PUSH_ARGUMENT(buffer, ii, what); + PUSH_ARGUMENT(buffer, ii, narm); + PUSH_ARGUMENT(buffer, ii, hasna); + PUSH_ARGUMENT(buffer, ii, ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (nrows + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nrows) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nrows); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj; R_xlen_t colBegin, idx; int count; @@ -144,7 +241,7 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { count = ans[ii]; /* Nothing more to do on this row? */ if (count == NA_INTEGER) continue; - + idx = R_INDEX_OP(colBegin, +, ROW_INDEX(crows,ii)); xvalue = R_INDEX_GET(x, idx, X_NA); if (xvalue == value) { @@ -158,14 +255,14 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } /* for (ii ...) */ } /* for (jj ...) */ } - } else { - error("INTERNAL ERROR: Unknown value of 'what' for colCounts: %d", what); - } /* if (what ...) */ + } /* if (what) */ } /*************************************************************************** HISTORY: + 2015-07-31 [DJ] + o Pthread processing. 2015-04-13 [DJ] o Supported subsetted computation. 2014-11-06 [HB] diff --git a/src/rowCumMinMaxs.c b/src/rowCumMinMaxs.c index f7756390..2fcc4d49 100644 --- a/src/rowCumMinMaxs.c +++ b/src/rowCumMinMaxs.c @@ -1,7 +1,6 @@ /*************************************************************************** Public methods: SEXP rowCummins(SEXP x, ...) - SEXP colCummins(SEXP x, ...) Authors: Henrik Bengtsson @@ -12,7 +11,7 @@ #include "utils.h" #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans, R_xlen_t nrow_ans, int *oks, int cores #define METHOD rowCummins #define COMP '<' @@ -23,8 +22,8 @@ #include "templates-gen-matrix.h" -SEXP rowCummins(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { - int byrow; +SEXP rowCummins(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow, SEXP cores) { + int byrow, cores2; SEXP ans = NILSXP; R_xlen_t nrow, ncol; @@ -42,14 +41,25 @@ SEXP rowCummins(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + + int *oks = NULL; + if (byrow) oks = (int *) R_alloc(nrows, sizeof(int)); /* Double matrices are more common to use. */ if (isReal(x)) { PROTECT(ans = allocMatrix(REALSXP, nrows, ncols)); - rowCummins_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, REAL(ans)); + rowCummins_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, REAL(ans), nrows, oks, cores2); UNPROTECT(1); } else if (isInteger(x)) { PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - rowCummins_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, INTEGER(ans)); + rowCummins_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, INTEGER(ans), nrows, oks, cores2); UNPROTECT(1); } @@ -61,9 +71,6 @@ SEXP rowCummins(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { -#define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans - #define METHOD rowCummaxs #define COMP '>' #define X_TYPE 'i' @@ -72,8 +79,8 @@ SEXP rowCummins(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { #include "templates-gen-matrix.h" -SEXP rowCummaxs(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { - int byrow; +SEXP rowCummaxs(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow, SEXP cores) { + int byrow, cores2; SEXP ans = NILSXP; R_xlen_t nrow, ncol; @@ -91,14 +98,25 @@ SEXP rowCummaxs(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + + int *oks = NULL; + if (byrow) oks = (int *) R_alloc(nrows, sizeof(int)); /* Double matrices are more common to use. */ if (isReal(x)) { PROTECT(ans = allocMatrix(REALSXP, nrows, ncols)); - rowCummaxs_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, REAL(ans)); + rowCummaxs_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, REAL(ans), nrows, oks, cores2); UNPROTECT(1); } else if (isInteger(x)) { PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - rowCummaxs_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, INTEGER(ans)); + rowCummaxs_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, INTEGER(ans), nrows, oks, cores2); UNPROTECT(1); } @@ -110,6 +128,8 @@ SEXP rowCummaxs(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { /*************************************************************************** HISTORY: + 2015-08-01 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-26 [HB] diff --git a/src/rowCumMinMaxs_TYPE-template.h b/src/rowCumMinMaxs_TYPE-template.h index 4f047f4c..1e682028 100644 --- a/src/rowCumMinMaxs_TYPE-template.h +++ b/src/rowCumMinMaxs_TYPE-template.h @@ -1,9 +1,9 @@ /*********************************************************************** TEMPLATE: - void rowCummins_[rowsType][colsType](ARGUMENTS_LIST) + void rowCum_[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans, R_xlen_t nrow_ans, int *oks, int cores Arguments: The following macros ("arguments") should be defined for the @@ -19,6 +19,8 @@ ***********************************************************************/ #include #include +#include +#include #include "types.h" /* Expand arguments: @@ -32,12 +34,122 @@ #define OP > #endif + +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, int, byrow); + POP_ARGUMENT(buffer0, ii, ANS_C_TYPE*, ans); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow_ans); + POP_ARGUMENT(buffer0, ii, int*, oks); + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + int rowsType = ROWS_TYPE_CODE, colsType = COLS_TYPE_CODE; + if (byrow) { + ans += begin; + nrows = end - begin; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + oks += begin; + + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, byrow, ans, nrow_ans, oks, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + + } else { + ans += begin * nrow_ans; + ncols = end - begin; +#ifdef COLS_TYPE + cols = (COLS_C_TYPE*) cols + begin; +#else + cols = indicesFromRange(begin, end, &colsType); +#endif + + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, byrow, ans, nrow_ans, oks, 0); + +#ifndef COLS_TYPE + Free(cols); +#endif + } + return NULL; +} + + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + R_xlen_t nv = byrow ? nrows : ncols; + if (cores > 1 && nv > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(byrow) + sizeof(ans) + sizeof(nrow_ans) + sizeof(oks); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nv); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, byrow); + PUSH_ARGUMENT(buffer, ii, ans); + PUSH_ARGUMENT(buffer, ii, nrow_ans); + PUSH_ARGUMENT(buffer, ii, oks); + + pthread_t threads[cores]; + + R_xlen_t gap = (nv + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nv) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nv); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj, kk, kk_prev, idx; R_xlen_t colBegin; ANS_C_TYPE value; int ok; - int *oks = NULL; #ifdef ROWS_TYPE ROWS_C_TYPE *crows = (ROWS_C_TYPE*) rows; @@ -49,8 +161,6 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { if (ncols == 0 || nrows == 0) return; if (byrow) { - oks = (int *) R_alloc(nrows, sizeof(int)); - colBegin = R_INDEX_OP(COL_INDEX(ccols,0), *, nrow); for (kk=0; kk < nrows; kk++) { idx = R_INDEX_OP(colBegin, +, ROW_INDEX(crows,kk)); @@ -65,9 +175,10 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } } - kk_prev = 0; for (jj=1; jj < ncols; jj++) { colBegin = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); + kk_prev = kk - nrows; + kk = kk_prev + nrow_ans; for (ii=0; ii < nrows; ii++) { idx = R_INDEX_OP(colBegin, +, ROW_INDEX(crows,ii)); value = (ANS_C_TYPE) R_INDEX_GET(x, idx, X_NA); @@ -88,7 +199,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { kk++; kk_prev++; - R_CHECK_USER_INTERRUPT(kk); + // TODO: interrupt subthreads + if (cores) R_CHECK_USER_INTERRUPT(kk); } /* for (ii ...) */ } /* for (jj ...) */ } else { @@ -130,7 +242,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { kk++; } - R_CHECK_USER_INTERRUPT(kk); + // TODO: interrupt subthreads + if (cores) R_CHECK_USER_INTERRUPT(kk); } /* for (ii ...) */ } /* for (jj ...) */ } /* if (byrow) */ @@ -141,6 +254,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { /*************************************************************************** HISTORY: + 2015-08-01 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-26 [HB] diff --git a/src/rowCumprods.c b/src/rowCumprods.c index 98d97b25..4976d3f3 100644 --- a/src/rowCumprods.c +++ b/src/rowCumprods.c @@ -1,7 +1,6 @@ /*************************************************************************** Public methods: SEXP rowCumprods(SEXP x, ...) - SEXP colCumprods(SEXP x, ...) Authors: Henrik Bengtsson @@ -13,7 +12,7 @@ #define METHOD rowCumprods #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans, R_xlen_t nrow_ans, int *oks, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -21,8 +20,8 @@ #include "templates-gen-matrix.h" -SEXP rowCumprods(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { - int byrow; +SEXP rowCumprods(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow, SEXP cores) { + int byrow, cores2; SEXP ans = NILSXP; R_xlen_t nrow, ncol; @@ -40,14 +39,25 @@ SEXP rowCumprods(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + + int *oks = NULL; /* Double matrices are more common to use. */ if (isReal(x)) { PROTECT(ans = allocMatrix(REALSXP, nrows, ncols)); - rowCumprods_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, REAL(ans)); + rowCumprods_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, REAL(ans), nrows, oks, cores2); UNPROTECT(1); } else if (isInteger(x)) { PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - rowCumprods_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, INTEGER(ans)); + if (byrow) oks = (int*)R_alloc(nrows, sizeof(int)); + rowCumprods_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, INTEGER(ans), nrows, oks, cores2); UNPROTECT(1); } @@ -57,6 +67,8 @@ SEXP rowCumprods(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { /*************************************************************************** HISTORY: + 2015-08-01 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-26 [HB] diff --git a/src/rowCumprods_TYPE-template.h b/src/rowCumprods_TYPE-template.h index 47962691..efd04da0 100644 --- a/src/rowCumprods_TYPE-template.h +++ b/src/rowCumprods_TYPE-template.h @@ -3,7 +3,7 @@ void rowCumprods_[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans, R_xlen_t nrow_ans, int *oks, int cores Arguments: The following macros ("arguments") should be defined for the @@ -19,6 +19,8 @@ ***********************************************************************/ #include #include +#include +#include #include "types.h" /* Expand arguments: @@ -27,7 +29,119 @@ #include "templates-types.h" +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, int, byrow); + POP_ARGUMENT(buffer0, ii, ANS_C_TYPE*, ans); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow_ans); + POP_ARGUMENT(buffer0, ii, int*, oks); + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + int rowsType = ROWS_TYPE_CODE, colsType = COLS_TYPE_CODE; + if (byrow) { + ans += begin; + nrows = end - begin; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + +#if ANS_TYPE == 'i' + oks += begin; +#endif + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, byrow, ans, nrow_ans, oks, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + + } else { + ans += begin * nrow_ans; + ncols = end - begin; +#ifdef COLS_TYPE + cols = (COLS_C_TYPE*) cols + begin; +#else + cols = indicesFromRange(begin, end, &colsType); +#endif + + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, byrow, ans, nrow_ans, oks, 0); + +#ifndef COLS_TYPE + Free(cols); +#endif + } + return NULL; +} + + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + R_xlen_t nv = byrow ? nrows : ncols; + if (cores > 1 && nv > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(byrow) + sizeof(ans) + sizeof(nrow_ans) + sizeof(oks); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nv); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, byrow); + PUSH_ARGUMENT(buffer, ii, ans); + PUSH_ARGUMENT(buffer, ii, nrow_ans); + PUSH_ARGUMENT(buffer, ii, oks); + + pthread_t threads[cores]; + + R_xlen_t gap = (nv + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nv) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nv); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj, kk, kk_prev, idx; R_xlen_t colBegin; X_C_TYPE xvalue; @@ -41,19 +155,15 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { #endif #if ANS_TYPE == 'i' - double R_INT_MIN_d = (double)R_INT_MIN, + double R_INT_MIN_d = (double)R_INT_MIN, R_INT_MAX_d = (double)R_INT_MAX; /* OK, i.e. no integer overflow yet? */ - int warn = 0, ok, *oks = NULL; + int warn = 0, ok; #endif if (ncols == 0 || nrows == 0) return; if (byrow) { -#if ANS_TYPE == 'i' - oks = (int *) R_alloc(nrows, sizeof(int)); -#endif - colBegin = R_INDEX_OP(COL_INDEX(ccols,0), *, nrow); for (kk=0; kk < nrows; kk++) { idx = R_INDEX_OP(colBegin, +, ROW_INDEX(crows,kk)); @@ -64,9 +174,10 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { #endif } - kk_prev = 0; for (jj=1; jj < ncols; jj++) { colBegin = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); + kk_prev = kk - nrows; + kk = kk_prev + nrow_ans; for (ii=0; ii < nrows; ii++) { idx = R_INDEX_OP(colBegin, +, ROW_INDEX(crows,ii)); xvalue = R_INDEX_GET(x, idx, X_NA); @@ -93,10 +204,11 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { ans[kk] = (ANS_C_TYPE) ((LDOUBLE) ans[kk_prev] * (LDOUBLE) xvalue); #endif - kk++; + kk++; kk_prev++; - R_CHECK_USER_INTERRUPT(kk); + // TODO: interrupt subthreads + if (cores) R_CHECK_USER_INTERRUPT(kk); } /* for (ii ...) */ } /* for (jj ...) */ } else { @@ -133,9 +245,10 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { value *= xvalue; ans[kk] = (ANS_C_TYPE) value; #endif - kk++; + kk++; - R_CHECK_USER_INTERRUPT(kk); + // TODO: interrupt subthreads + if (cores) R_CHECK_USER_INTERRUPT(kk); } /* for (ii ...) */ } /* for (jj ...) */ } /* if (byrow) */ @@ -143,14 +256,17 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { #if ANS_TYPE == 'i' /* Warn on integer overflow? */ if (warn) { - warning("Integer overflow. Detected one or more elements whose absolute values were out of the range [%d,%d] that can be used to for integers. Such values are set to NA_integer_.", R_INT_MIN, R_INT_MAX); + // TODO: warning in subthreads + if (cores) warning("Integer overflow. Detected one or more elements whose absolute values were out of the range [%d,%d] that can be used to for integers. Such values are set to NA_integer_.", R_INT_MIN, R_INT_MAX); } -#endif +#endif } /*************************************************************************** HISTORY: + 2015-08-01 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-26 [HB] diff --git a/src/rowCumsums.c b/src/rowCumsums.c index ea7fc3ae..c87b4340 100644 --- a/src/rowCumsums.c +++ b/src/rowCumsums.c @@ -1,7 +1,6 @@ /*************************************************************************** Public methods: SEXP rowCumsums(SEXP x, ...) - SEXP colCumsums(SEXP x, ...) Authors: Henrik Bengtsson @@ -13,7 +12,7 @@ #define METHOD rowCumsums #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans, R_xlen_t nrow_ans, int *oks, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -21,8 +20,8 @@ #include "templates-gen-matrix.h" -SEXP rowCumsums(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { - int byrow; +SEXP rowCumsums(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow, SEXP cores) { + int byrow, cores2; SEXP ans = NILSXP; R_xlen_t nrow, ncol; @@ -40,14 +39,25 @@ SEXP rowCumsums(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + + int *oks = NULL; /* Double matrices are more common to use. */ if (isReal(x)) { PROTECT(ans = allocMatrix(REALSXP, nrows, ncols)); - rowCumsums_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, REAL(ans)); + rowCumsums_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, REAL(ans), nrows, oks, cores2); UNPROTECT(1); } else if (isInteger(x)) { PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - rowCumsums_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, INTEGER(ans)); + if (byrow) oks = (int*)R_alloc(nrows, sizeof(int)); + rowCumsums_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, INTEGER(ans), nrows, oks, cores2); UNPROTECT(1); } @@ -57,6 +67,8 @@ SEXP rowCumsums(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP byRow) { /*************************************************************************** HISTORY: + 2015-08-01 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-26 [HB] diff --git a/src/rowCumsums_TYPE-template.h b/src/rowCumsums_TYPE-template.h index b435bc14..c7cd8633 100644 --- a/src/rowCumsums_TYPE-template.h +++ b/src/rowCumsums_TYPE-template.h @@ -3,7 +3,7 @@ void rowCumsums_[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, ANS_C_TYPE *ans, R_xlen_t nrow_ans, int *oks, int cores Arguments: The following macros ("arguments") should be defined for the @@ -19,6 +19,8 @@ ***********************************************************************/ #include #include +#include +#include #include "types.h" /* Expand arguments: @@ -27,7 +29,120 @@ #include "templates-types.h" +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, int, byrow); + POP_ARGUMENT(buffer0, ii, ANS_C_TYPE*, ans); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow_ans); + POP_ARGUMENT(buffer0, ii, int*, oks); + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + int rowsType = ROWS_TYPE_CODE, colsType = COLS_TYPE_CODE; + if (byrow) { + ans += begin; + nrows = end - begin; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + +#if ANS_TYPE == 'i' + oks += begin; +#endif + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, byrow, ans, nrow_ans, oks, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + + } else { + ans += begin * nrow_ans; + ncols = end - begin; +#ifdef COLS_TYPE + cols = (COLS_C_TYPE*) cols + begin; +#else + cols = indicesFromRange(begin, end, &colsType); +#endif + + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, byrow, ans, nrow_ans, oks, 0); + +#ifndef COLS_TYPE + Free(cols); +#endif + } + return NULL; +} + + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + R_xlen_t nv = byrow ? nrows : ncols; + if (cores > 1 && nv > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(byrow) + sizeof(ans) + sizeof(nrow_ans) + sizeof(oks); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nv); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, byrow); + PUSH_ARGUMENT(buffer, ii, ans); + PUSH_ARGUMENT(buffer, ii, nrow_ans); + PUSH_ARGUMENT(buffer, ii, oks); + + pthread_t threads[cores]; + + R_xlen_t gap = (nv + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nv) { + uint8_t *args = buffer + ii; + buffer0 = buffer; + end = MIN(begin + gap, nv); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj, kk, kk_prev, idx; R_xlen_t colBegin; X_C_TYPE xvalue; @@ -41,19 +156,15 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { #endif #if ANS_TYPE == 'i' - double R_INT_MIN_d = (double)R_INT_MIN, + double R_INT_MIN_d = (double)R_INT_MIN, R_INT_MAX_d = (double)R_INT_MAX; /* OK, i.e. no integer overflow yet? */ - int warn = 0, ok, *oks = NULL; + int warn = 0, ok; #endif if (ncols == 0 || nrows == 0) return; if (byrow) { -#if ANS_TYPE == 'i' - oks = (int *) R_alloc(nrows, sizeof(int)); -#endif - colBegin = R_INDEX_OP(COL_INDEX(ccols,0), *, nrow); for (kk=0; kk < nrows; kk++) { idx = R_INDEX_OP(colBegin, +, ROW_INDEX(crows,kk)); @@ -64,9 +175,10 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { #endif } - kk_prev = 0; for (jj=1; jj < ncols; jj++) { colBegin = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); + kk_prev = kk - nrows; + kk = kk_prev + nrow_ans; for (ii=0; ii < nrows; ii++) { idx = R_INDEX_OP(colBegin, +, ROW_INDEX(crows,ii)); xvalue = R_INDEX_GET(x, idx, X_NA); @@ -94,10 +206,11 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { ans[kk] = (ANS_C_TYPE) ((LDOUBLE) ans[kk_prev] + (LDOUBLE) xvalue); #endif - kk++; - kk_prev++; + kk++; + kk_prev++; - R_CHECK_USER_INTERRUPT(kk); + // TODO: interrupt subthreads + if (cores) R_CHECK_USER_INTERRUPT(kk); } /* for (ii ...) */ } /* for (jj ...) */ } else { @@ -136,9 +249,10 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { ans[kk] = (ANS_C_TYPE) value; #endif - kk++; + kk++; - R_CHECK_USER_INTERRUPT(kk); + // TODO: interrupt subthreads + if (cores) R_CHECK_USER_INTERRUPT(kk); } /* for (ii ...) */ } /* for (jj ...) */ } /* if (byrow) */ @@ -146,14 +260,17 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { #if ANS_TYPE == 'i' /* Warn on integer overflow? */ if (warn) { - warning("Integer overflow. Detected one or more elements whose absolute values were out of the range [%d,%d] that can be used to for integers. Such values are set to NA_integer_.", R_INT_MIN, R_INT_MAX); + // TODO: warning in subthreads + if (cores) warning("Integer overflow. Detected one or more elements whose absolute values were out of the range [%d,%d] that can be used to for integers. Such values are set to NA_integer_.", R_INT_MIN, R_INT_MAX); } -#endif +#endif } /*************************************************************************** HISTORY: + 2015-08-01 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-26 [HB] diff --git a/src/rowDiffs.c b/src/rowDiffs.c index 7b456560..af4f8e39 100644 --- a/src/rowDiffs.c +++ b/src/rowDiffs.c @@ -1,7 +1,6 @@ /*************************************************************************** Public methods: SEXP rowDiffs(SEXP x, ...) - SEXP colDiffs(SEXP x, ...) Authors: Henrik Bengtsson @@ -13,7 +12,7 @@ #define METHOD rowDiffs #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, R_xlen_t lag, R_xlen_t differences, X_C_TYPE *ans, R_xlen_t nrow_ans, R_xlen_t ncol_ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, R_xlen_t lag, R_xlen_t differences, X_C_TYPE *ans, R_xlen_t nrow_ans, R_xlen_t ncol_ans, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -21,8 +20,8 @@ #include "templates-gen-matrix.h" -SEXP rowDiffs(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP lag, SEXP differences, SEXP byRow) { - int byrow; +SEXP rowDiffs(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP lag, SEXP differences, SEXP byRow, SEXP cores) { + int byrow, cores2; SEXP ans = NILSXP; R_xlen_t lagg, diff; R_xlen_t nrow, ncol; @@ -54,25 +53,33 @@ SEXP rowDiffs(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP lag, SEXP differences /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif /* Dimension of result matrix */ if (byrow) { nrow_ans = nrows; ncol_ans = (R_xlen_t)((double)ncols - ((double)diff*(double)lagg)); - if (ncol_ans < 0) ncol_ans = 0; + if (ncol_ans < 0) ncol_ans = 0; } else { nrow_ans = (R_xlen_t)((double)nrows - ((double)diff*(double)lagg)); - if (nrow_ans < 0) nrow_ans = 0; + if (nrow_ans < 0) nrow_ans = 0; ncol_ans = ncols; } if (isReal(x)) { PROTECT(ans = allocMatrix(REALSXP, nrow_ans, ncol_ans)); - rowDiffs_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, lagg, diff, REAL(ans), nrow_ans, ncol_ans); + rowDiffs_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, lagg, diff, REAL(ans), nrow_ans, ncol_ans, cores2); UNPROTECT(1); } else if (isInteger(x)) { PROTECT(ans = allocMatrix(INTSXP, nrow_ans, ncol_ans)); - rowDiffs_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, lagg, diff, INTEGER(ans), nrow_ans, ncol_ans); + rowDiffs_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, byrow, lagg, diff, INTEGER(ans), nrow_ans, ncol_ans, cores2); UNPROTECT(1); } @@ -82,8 +89,10 @@ SEXP rowDiffs(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP lag, SEXP differences /*************************************************************************** HISTORY: + 2015-08-02 [DJ] + o Pthread processing. 2015-06-13 [DJ] o Supported subsetted computation. 2014-12-29 [HB] - o Created. + o Created. **************************************************************************/ diff --git a/src/rowDiffs_TYPE-template.h b/src/rowDiffs_TYPE-template.h index 61c68ce5..b242fba5 100644 --- a/src/rowDiffs_TYPE-template.h +++ b/src/rowDiffs_TYPE-template.h @@ -3,7 +3,7 @@ void rowDiffs_(ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, R_xlen_t lag, R_xlen_t differences, X_C_TYPE *ans, R_xlen_t nrow_ans, R_xlen_t ncol_ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int byrow, R_xlen_t lag, R_xlen_t differences, X_C_TYPE *ans, R_xlen_t nrow_ans, R_xlen_t ncol_ans, int cores Arguments: The following macros ("arguments") should be defined for the @@ -14,6 +14,8 @@ Copyright: Henrik Bengtsson, 2014 ***********************************************************************/ +#include +#include #include "types.h" /* Expand arguments: @@ -49,15 +51,16 @@ tt = 0; ss = 0; for (jj=0; jj < ncol_y; jj++) { - for (ii=0; ii < nrow_y; ii++) { + for (ii=0; ii < nrow_x; ii++) { y[ss++] = X_DIFF(x[uu++], x[tt++]); } + ss += nrow_y - nrow_x; } } else { uu = lag; tt = 0; ss = 0; - for (jj=0; jj < ncol_y; jj++) { + for (jj=0; jj < ncol_x; jj++) { for (ii=0; ii < nrow_y; ii++) { /* Rprintf("y[%d] = x[%d] - x[%d] = %g - %g = %g\n", ss, uu, tt, (double)x[uu], (double)x[tt], (double)X_DIFF(x[uu], x[tt])); */ y[ss++] = X_DIFF(x[uu++], x[tt++]); @@ -117,7 +120,7 @@ static R_INLINE void DIFF_X_MATRIX_ROWS_COLS(X_C_TYPE *x, R_xlen_t nrow, void *r colBegin1 = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); colBegin2 = R_INDEX_OP(COL_INDEX(ccols,(jj+lag)), *, nrow); - for (ii=0; ii < nrow_ans; ii++) { + for (ii=0; ii < nrows; ii++) { idx = R_INDEX_OP(colBegin1, +, ROW_INDEX(crows,ii)); xvalue1 = R_INDEX_GET(x, idx, X_NA); idx = R_INDEX_OP(colBegin2, +, ROW_INDEX(crows,ii)); @@ -125,9 +128,10 @@ static R_INLINE void DIFF_X_MATRIX_ROWS_COLS(X_C_TYPE *x, R_xlen_t nrow, void *r ans[ss++] = X_DIFF(xvalue2, xvalue1); } + ss += nrow_ans - nrows; } } else { - for (jj=0; jj < ncol_ans; jj++) { + for (jj=0; jj < ncols; jj++) { colBegin1 = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); for (ii=0; ii < nrow_ans; ii++) { @@ -143,8 +147,120 @@ static R_INLINE void DIFF_X_MATRIX_ROWS_COLS(X_C_TYPE *x, R_xlen_t nrow, void *r } +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, int, byrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, lag); + POP_ARGUMENT(buffer0, ii, R_xlen_t, differences); + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, ans); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow_ans); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol_ans); + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + int rowsType = ROWS_TYPE_CODE, colsType = COLS_TYPE_CODE; + if (byrow) { + ans += begin; + nrows = end - begin; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, byrow, lag, differences, ans, nrow_ans, ncol_ans, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + + } else { + ans += begin * nrow_ans; + ncols = end - begin; +#ifdef COLS_TYPE + cols = (COLS_C_TYPE*) cols + begin; +#else + cols = indicesFromRange(begin, end, &colsType); +#endif + + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, byrow, lag, differences, ans, nrow_ans, ncol_ans, 0); + +#ifndef COLS_TYPE + Free(cols); +#endif + } + return NULL; +} + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + R_xlen_t nv = byrow ? nrows : ncols; + if (cores > 1 && nv > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(byrow) + sizeof(lag) + sizeof(differences) + sizeof(ans) + sizeof(nrow_ans) + sizeof(ncol_ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nv); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, byrow); + PUSH_ARGUMENT(buffer, ii, lag); + PUSH_ARGUMENT(buffer, ii, differences); + PUSH_ARGUMENT(buffer, ii, ans); + PUSH_ARGUMENT(buffer, ii, nrow_ans); + PUSH_ARGUMENT(buffer, ii, ncol_ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (nv + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nv) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nv); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t nrow_tmp, ncol_tmp; X_C_TYPE *tmp = NULL; @@ -194,8 +310,10 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { /*************************************************************************** HISTORY: + 2015-08-02 [DJ] + o Pthread processing. 2015-06-13 [DJ] o Supported subsetted computation. 2014-12-29 [HB] - o Created. + o Created. **************************************************************************/ diff --git a/src/rowLogSumExp.c b/src/rowLogSumExp.c index 69e812e2..907d314a 100644 --- a/src/rowLogSumExp.c +++ b/src/rowLogSumExp.c @@ -1,6 +1,6 @@ /*************************************************************************** Public methods: - SEXP rowLogSumExps(SEXP lx, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP byRow) + SEXP rowLogSumExps(SEXP lx, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP byRow, SEXP cores) Authors: Henrik Bengtsson @@ -14,14 +14,14 @@ #define METHOD rowLogSumExp #define METHOD_NAME rowLogSumExps_double #define RETURN_TYPE void -#define ARGUMENTS_LIST double *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, int rowsType, void *cols, R_xlen_t ncols, int colsType, int narm, int hasna, R_xlen_t byrow, double *ans +#define ARGUMENTS_LIST double *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, int rowsType, void *cols, R_xlen_t ncols, int colsType, int narm, int hasna, int byrow, double *ans, int cores #include "templates-gen-vector.h" -SEXP rowLogSumExps(SEXP lx, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP byRow) { +SEXP rowLogSumExps(SEXP lx, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP byRow, SEXP cores) { SEXP ans; - int narm, hasna, byrow; + int narm, hasna, byrow, cores2; R_xlen_t nrow, ncol; /* Argument 'lx' and 'dim': */ @@ -44,12 +44,21 @@ SEXP rowLogSumExps(SEXP lx, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasN /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + if (byrow) { ans = PROTECT(allocVector(REALSXP, nrows)); - rowLogSumExps_double[rowsType](REAL(lx), nrow, ncol, crows, nrows, rowsType, ccols, ncols, colsType, narm, hasna, 1, REAL(ans)); + rowLogSumExps_double[rowsType](REAL(lx), nrow, ncol, crows, nrows, rowsType, ccols, ncols, colsType, narm, hasna, 1, REAL(ans), cores2); } else { ans = PROTECT(allocVector(REALSXP, ncols)); - rowLogSumExps_double[colsType](REAL(lx), nrow, ncol, crows, nrows, rowsType, ccols, ncols, colsType, narm, hasna, 0, REAL(ans)); + rowLogSumExps_double[colsType](REAL(lx), nrow, ncol, crows, nrows, rowsType, ccols, ncols, colsType, narm, hasna, 0, REAL(ans), cores2); } UNPROTECT(1); /* ans = PROTECT(...) */ @@ -60,6 +69,8 @@ SEXP rowLogSumExps(SEXP lx, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasN /*************************************************************************** HISTORY: + 2015-08-01 [DJ] + o Pthread processing. 2015-06-12 [DJ] o Supported subsetted computation. 2013-05-02 [HB] diff --git a/src/rowLogSumExp_TYPE-template.h b/src/rowLogSumExp_TYPE-template.h index 62810c25..0a957b4b 100644 --- a/src/rowLogSumExp_TYPE-template.h +++ b/src/rowLogSumExp_TYPE-template.h @@ -3,21 +3,131 @@ double rowLogSumExp_double[idxsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - double *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, int rowsType, void *cols, R_xlen_t ncols, int colsType, int narm, int hasna, R_xlen_t byrow, double *ans + double *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, int rowsType, void *cols, R_xlen_t ncols, int colsType, int narm, int hasna, int byrow, double *ans, int cores ***********************************************************************/ +#include +#include #include "types.h" #include "templates-types.h" -/* extern 1-D function 'logSumExp' */ -extern double (*logSumExp_double[3])(double *x, void *idxs, R_xlen_t nidxs, int narm, int hasna, int by, double *xx); +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_IDXS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, double*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, int, rowsType); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, int, colsType); + POP_ARGUMENT(buffer0, ii, int, narm); + POP_ARGUMENT(buffer0, ii, int, hasna); + POP_ARGUMENT(buffer0, ii, int, byrow); + POP_ARGUMENT(buffer0, ii, double*, ans); + + ans += begin; + extern RETURN_TYPE (*METHOD_NAME[3])(ARGUMENTS_LIST); + if (byrow) { + nrows = end - begin; +#ifdef IDXS_TYPE + rows = (IDXS_C_TYPE*) rows + begin; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + + METHOD_NAME[rowsType](x, nrow, ncol, rows, nrows, rowsType, cols, ncols, colsType, narm, hasna, byrow, ans, 0); + +#ifndef IDXS_TYPE + Free(rows); +#endif + + } else { + ncols = end - begin; +#ifdef IDXS_TYPE + cols = (IDXS_C_TYPE*) cols + begin; +#else + cols = indicesFromRange(begin, end, &colsType); +#endif + + METHOD_NAME[colsType](x, nrow, ncol, rows, nrows, rowsType, cols, ncols, colsType, narm, hasna, byrow, ans, 0); + +#ifndef IDXS_TYPE + Free(cols); +#endif + } + return NULL; +} RETURN_TYPE METHOD_NAME_IDXS(ARGUMENTS_LIST) { + // Apply pthread + R_xlen_t nv = byrow ? nrows : ncols; + if (cores > 1 && nv > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(rowsType) + sizeof(cols) + sizeof(ncols) + sizeof(colsType) + + sizeof(narm) + sizeof(hasna) + sizeof(byrow) + sizeof(ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nv); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, rowsType); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, colsType); + PUSH_ARGUMENT(buffer, ii, narm); + PUSH_ARGUMENT(buffer, ii, hasna); + PUSH_ARGUMENT(buffer, ii, byrow); + PUSH_ARGUMENT(buffer, ii, ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (nv + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nv) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nv); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_IDXS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, idx; double navalue; - double (*logsumexp)(double *x, void *idxs, R_xlen_t nidxs, int narm, int hasna, int by, double *xx); + /* extern 1-D function 'logSumExp' */ + extern double (*logSumExp_double[3])(double *x, void *idxs, R_xlen_t nidxs, int narm, int hasna, int by, double *xx); #ifdef IDXS_TYPE IDXS_C_TYPE *crows = (IDXS_C_TYPE*) rows; @@ -25,31 +135,30 @@ RETURN_TYPE METHOD_NAME_IDXS(ARGUMENTS_LIST) { #endif if (byrow) { - /* R allocate memory for row-vector 'xx' of length 'ncol'. - This will be taken care of by the R garbage collector later on. */ - double *xx = (double *) R_alloc(ncols, sizeof(double)); - + /* C allocate memory for row-vector 'xx' of length 'ncol'. + Should be freed manually. */ + double *xx = Calloc(ncols, double); navalue = (narm || ncols == 0) ? R_NegInf : NA_REAL; - logsumexp = logSumExp_double[colsType]; for (ii=0; ii < nrows; ++ii) { idx = IDX_INDEX(crows,ii); if (idx == NA_R_XLEN_T) { ans[ii] = navalue; } else { - ans[ii] = logsumexp(x+idx, cols, ncols, narm, hasna, nrow, xx); + ans[ii] = logSumExp_double[colsType](x+idx, cols, ncols, narm, hasna, nrow, xx); } } + Free(xx); + } else { navalue = (narm || nrows == 0) ? R_NegInf : NA_REAL; - logsumexp = logSumExp_double[rowsType]; for (ii=0; ii < ncols; ++ii) { idx = R_INDEX_OP(IDX_INDEX(ccols,ii), *, nrow); if (idx == NA_R_XLEN_T) { ans[ii] = navalue; } else { - ans[ii] = logsumexp(x+idx, rows, nrows, narm, hasna, 0, NULL); + ans[ii] = logSumExp_double[rowsType](x+idx, rows, nrows, narm, hasna, 0, NULL); } } } /* if (byrow) */ @@ -58,6 +167,8 @@ RETURN_TYPE METHOD_NAME_IDXS(ARGUMENTS_LIST) { /*************************************************************************** HISTORY: + 2015-08-01 [DJ] + o Pthread processing. 2013-06-12 [DH] o Created. **************************************************************************/ diff --git a/src/rowMads.c b/src/rowMads.c index a289f3a8..04956efc 100644 --- a/src/rowMads.c +++ b/src/rowMads.c @@ -1,7 +1,6 @@ /*************************************************************************** Public methods: SEXP rowMads(SEXP x, ...) - SEXP colMads(SEXP x, ...) Authors: Henrik Bengtsson @@ -14,7 +13,7 @@ #define METHOD rowMads #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, double scale, int narm, int hasna, int byrow, double *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, double scale, int narm, int hasna, int byrow, double *ans, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -22,8 +21,8 @@ #include "templates-gen-matrix.h" -SEXP rowMads(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP constant, SEXP naRm, SEXP hasNA, SEXP byRow) { - int narm, hasna, byrow; +SEXP rowMads(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP constant, SEXP naRm, SEXP hasNA, SEXP byRow, SEXP cores) { + int narm, hasna, byrow, cores2; SEXP ans; R_xlen_t nrow, ncol; double scale; @@ -53,6 +52,15 @@ SEXP rowMads(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP constant, SEXP naRm, S /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + if (!byrow) { SWAP(R_xlen_t, nrow, ncol); SWAP(void*, crows, ccols); @@ -66,9 +74,9 @@ SEXP rowMads(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP constant, SEXP naRm, S /* Double matrices are more common to use. */ if (isReal(x)) { - rowMads_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, scale, narm, hasna, byrow, REAL(ans)); + rowMads_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, scale, narm, hasna, byrow, REAL(ans), cores2); } else if (isInteger(x)) { - rowMads_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, scale, narm, hasna, byrow, REAL(ans)); + rowMads_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, scale, narm, hasna, byrow, REAL(ans), cores2); } UNPROTECT(1); @@ -79,6 +87,8 @@ SEXP rowMads(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP constant, SEXP naRm, S /*************************************************************************** HISTORY: + 2015-08-13 [DJ] + o Pthread processing. 2015-06-13 [DJ] o Supported subsetted computation. 2014-11-17 [HB] diff --git a/src/rowMads_TYPE-template.h b/src/rowMads_TYPE-template.h index a1cff71e..db20675a 100644 --- a/src/rowMads_TYPE-template.h +++ b/src/rowMads_TYPE-template.h @@ -3,7 +3,7 @@ void rowMads_(ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, double scale, int narm, int hasna, int byrow, double *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, double scale, int narm, int hasna, int byrow, double *ans, int cores Arguments: The following macros ("arguments") should be defined for the @@ -20,6 +20,8 @@ ***********************************************************************/ #include #include +#include +#include #include "types.h" #include /* abs() and fabs() */ @@ -29,7 +31,100 @@ #include "templates-types.h" +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE *, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, double, scale); + POP_ARGUMENT(buffer0, ii, int, narm); + POP_ARGUMENT(buffer0, ii, int, hasna); + POP_ARGUMENT(buffer0, ii, int, byrow); + POP_ARGUMENT(buffer0, ii, double*, ans); + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + int rowsType = ROWS_TYPE_CODE; + + ans += begin; + nrows = end - begin; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + + METHOD_NAME[rowsType][COLS_TYPE_CODE](x, nrow, ncol, rows, nrows, cols, ncols, scale, narm, hasna, byrow, ans, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + return NULL; +} + + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + if (cores > 1 && nrows > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(scale) + sizeof(narm) + sizeof(hasna) + sizeof(byrow) + sizeof(ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nrows); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, scale); + PUSH_ARGUMENT(buffer, ii, narm); + PUSH_ARGUMENT(buffer, ii, hasna); + PUSH_ARGUMENT(buffer, ii, byrow); + PUSH_ARGUMENT(buffer, ii, ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (nrows + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nrows) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nrows); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + int isOdd; R_xlen_t ii, jj, kk, qq, idx; R_xlen_t *colOffset; @@ -43,10 +138,9 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { COLS_C_TYPE *ccols = (COLS_C_TYPE*) cols; #endif - /* R allocate memory for the 'values'. This will be - taken care of by the R garbage collector later on. */ - values = (X_C_TYPE *) R_alloc(ncols, sizeof(X_C_TYPE)); - values_d = (double *) R_alloc(ncols, sizeof(double)); + /* C allocate memory for the 'values'. Should be freed manually. */ + values = Calloc(ncols, X_C_TYPE); + values_d = Calloc(ncols, double); /* If there are no missing values, don't try to remove them. */ if (hasna == FALSE) @@ -64,7 +158,7 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { value = 0; /* Pre-calculate the column offsets */ - colOffset = (R_xlen_t *) R_alloc(ncols, sizeof(R_xlen_t)); + colOffset = Calloc(ncols, R_xlen_t); // HJ begin if (byrow) { @@ -192,8 +286,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { ans[ii] = scale * (values_d[qq] + values_d[qq+1])/2; } } /* if (kk == 0) */ - - R_CHECK_USER_INTERRUPT(ii); + // TODO: interrupt under pthread + // R_CHECK_USER_INTERRUPT(ii); } /* for (ii ...) */ } else { for (ii=0; ii < nrows; ii++) { @@ -215,15 +309,21 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { X_PSORT(values, qq+1, qq); ans[ii] = ((double)values[qq] + value)/2; } - - R_CHECK_USER_INTERRUPT(ii); + // TODO: interrupt under pthread + // R_CHECK_USER_INTERRUPT(ii); } /* for (ii ...) */ } /* if (hasna ...) */ + + Free(values); + Free(values_d); + Free(colOffset); } /*************************************************************************** HISTORY: + 2015-08-13 [DJ] + o Pthread processing. 2015-06-13 [DJ] o Supported subsetted computation. 2014-11-17 [HB] diff --git a/src/rowMedians.c b/src/rowMedians.c index aef76562..fa889822 100644 --- a/src/rowMedians.c +++ b/src/rowMedians.c @@ -12,7 +12,7 @@ #define METHOD rowMedians #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int narm, int hasna, int byrow, double *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int narm, int hasna, int byrow, double *ans, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -20,8 +20,8 @@ #include "templates-gen-matrix.h" -SEXP rowMedians(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP byRow) { - int narm, hasna, byrow; +SEXP rowMedians(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP byRow, SEXP cores) { + int narm, hasna, byrow, cores2; SEXP ans; R_xlen_t nrow, ncol; @@ -46,6 +46,15 @@ SEXP rowMedians(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, S /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + if (!byrow) { SWAP(R_xlen_t, nrow, ncol); SWAP(void*, crows, ccols); @@ -59,9 +68,9 @@ SEXP rowMedians(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, S /* Double matrices are more common to use. */ if (isReal(x)) { - rowMedians_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, narm, hasna, byrow, REAL(ans)); + rowMedians_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, narm, hasna, byrow, REAL(ans), cores2); } else if (isInteger(x)) { - rowMedians_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, narm, hasna, byrow, REAL(ans)); + rowMedians_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, narm, hasna, byrow, REAL(ans), cores2); } UNPROTECT(1); @@ -72,6 +81,8 @@ SEXP rowMedians(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, S /*************************************************************************** HISTORY: + 2015-08-13 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2013-01-13 [HB] diff --git a/src/rowMedians_TYPE-template.h b/src/rowMedians_TYPE-template.h index 9ef82b2e..12bfd7b2 100644 --- a/src/rowMedians_TYPE-template.h +++ b/src/rowMedians_TYPE-template.h @@ -3,7 +3,7 @@ void rowMedians_[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int narm, int hasna, int byrow, double *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int narm, int hasna, int byrow, double *ans, int cores Arguments: The following macros ("arguments") should be defined for the @@ -20,6 +20,8 @@ ***********************************************************************/ #include #include +#include +#include #include "types.h" /* Expand arguments: @@ -28,7 +30,97 @@ #include "templates-types.h" +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE *, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, int, narm); + POP_ARGUMENT(buffer0, ii, int, hasna); + POP_ARGUMENT(buffer0, ii, int, byrow); + POP_ARGUMENT(buffer0, ii, double*, ans); + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + int rowsType = ROWS_TYPE_CODE; + + ans += begin; + nrows = end - begin; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + + METHOD_NAME[rowsType][COLS_TYPE_CODE](x, nrow, ncol, rows, nrows, cols, ncols, narm, hasna, byrow, ans, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + return NULL; +} + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + if (cores > 1 && nrows > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(narm) + sizeof(hasna) + sizeof(byrow) + sizeof(ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nrows); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, narm); + PUSH_ARGUMENT(buffer, ii, hasna); + PUSH_ARGUMENT(buffer, ii, byrow); + PUSH_ARGUMENT(buffer, ii, ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (nrows + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nrows) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nrows); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + int isOdd; R_xlen_t ii, jj, kk, qq, idx; R_xlen_t *colOffset; @@ -41,9 +133,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { COLS_C_TYPE *ccols = (COLS_C_TYPE*) cols; #endif - /* R allocate memory for the 'values'. This will be - taken care of by the R garbage collector later on. */ - values = (X_C_TYPE *) R_alloc(ncols, sizeof(X_C_TYPE)); + /* C allocate memory for the 'values'. Should be freed manually. */ + values = Calloc(ncols, X_C_TYPE); /* If there are no missing values, don't try to remove them. */ if (hasna == FALSE) @@ -61,7 +152,7 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { value = 0; /* Pre-calculate the column offsets */ - colOffset = (R_xlen_t *) R_alloc(ncols, sizeof(R_xlen_t)); + colOffset = Calloc(ncols, R_xlen_t); // HJ begin if (byrow) { @@ -120,8 +211,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { ans[ii] = ((double)values[qq] + (double)value)/2; } } - - R_CHECK_USER_INTERRUPT(ii); + // TODO: interrupt under pthread + // R_CHECK_USER_INTERRUPT(ii); } /* for (ii ...) */ } else { for (ii=0; ii < nrows; ii++) { @@ -143,15 +234,20 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { X_PSORT(values, qq+1, qq); ans[ii] = ((double)values[qq] + (double)value)/2; } - - R_CHECK_USER_INTERRUPT(ii); + // TODO: interrupt under pthread + // R_CHECK_USER_INTERRUPT(ii); } /* for (ii ...) */ } /* if (hasna ...) */ + + Free(values); + Free(colOffset); } /*************************************************************************** HISTORY: + 2015-08-13 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-06 [HB] diff --git a/src/rowOrderStats.c b/src/rowOrderStats.c index de298879..19e72d97 100644 --- a/src/rowOrderStats.c +++ b/src/rowOrderStats.c @@ -1,6 +1,6 @@ /*************************************************************************** Public methods: - SEXP rowOrderStats(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP which) + SEXP rowOrderStats(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP which, SEXP cores) Authors: Henrik Bengtsson. Adopted from rowQ() by R. Gentleman. @@ -14,7 +14,7 @@ #define METHOD rowOrderStats #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, R_xlen_t qq, X_C_TYPE *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, R_xlen_t qq, X_C_TYPE *ans, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -22,9 +22,10 @@ #include "templates-gen-matrix.h" -SEXP rowOrderStats(SEXP x, SEXP dim,SEXP rows, SEXP cols, SEXP which) { +SEXP rowOrderStats(SEXP x, SEXP dim,SEXP rows, SEXP cols, SEXP which, SEXP cores) { SEXP ans = NILSXP; R_xlen_t nrow, ncol, qq; + int cores2; /* Argument 'x' and 'dim': */ assertArgMatrix(x, dim, (R_TYPE_INT | R_TYPE_REAL), "x"); @@ -41,8 +42,18 @@ SEXP rowOrderStats(SEXP x, SEXP dim,SEXP rows, SEXP cols, SEXP which) { /* Argument 'rows' and 'cols': */ R_xlen_t nrows, ncols; int rowsType, colsType; - void *crows = validateIndices(rows, nrow, 0, &nrows, &rowsType); - void *ccols = validateIndices(cols, ncol, 0, &ncols, &colsType); + int rowsHasna, colsHasna; + void *crows = validateIndicesCheckNA(rows, nrow, 0, &nrows, &rowsType, &rowsHasna); + void *ccols = validateIndicesCheckNA(cols, ncol, 0, &ncols, &colsType, &colsHasna); + + // Check missing rows + if (rowsHasna && ncols > 0) { + error("Argument 'rows' must not contain missing value"); + } + // Check missing cols + if (colsHasna && nrows > 0) { + error("Argument 'cols' must not contain missing value"); + } /* Subtract one here, since rPsort does zero based addressing */ qq = asInteger(which) - 1; @@ -52,14 +63,23 @@ SEXP rowOrderStats(SEXP x, SEXP dim,SEXP rows, SEXP cols, SEXP which) { error("Argument 'which' is out of range."); } +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + /* Double matrices are more common to use. */ if (isReal(x)) { PROTECT(ans = allocVector(REALSXP, nrows)); - rowOrderStats_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, qq, REAL(ans)); + rowOrderStats_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, qq, REAL(ans), cores2); UNPROTECT(1); } else if (isInteger(x)) { PROTECT(ans = allocVector(INTSXP, nrows)); - rowOrderStats_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, qq, INTEGER(ans)); + rowOrderStats_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, qq, INTEGER(ans), cores2); UNPROTECT(1); } @@ -69,6 +89,8 @@ SEXP rowOrderStats(SEXP x, SEXP dim,SEXP rows, SEXP cols, SEXP which) { /*************************************************************************** HISTORY: + 2015-08-09 [DJ] + o Pthread processing. 2015-07-11 [DJ] o Supported subsetted computation. 2009-02-04 [HB] diff --git a/src/rowOrderStats_TYPE-template.h b/src/rowOrderStats_TYPE-template.h index d57ace64..615d36dd 100644 --- a/src/rowOrderStats_TYPE-template.h +++ b/src/rowOrderStats_TYPE-template.h @@ -3,7 +3,7 @@ void rowOrderStats_[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, R_xlen_t qq, X_C_TYPE *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, R_xlen_t qq, X_C_TYPE *ans, int cores Arguments: The following macros ("arguments") should be defined for the @@ -21,6 +21,8 @@ ***********************************************************************/ #include #include +#include +#include #include "types.h" /* Expand arguments: @@ -30,38 +32,112 @@ #include "templates-types.h" +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, qq); + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, ans); + + ans += begin; + + nrows = end - begin; + int rowsType; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; + rowsType = ROWS_TYPE_CODE; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + METHOD_NAME[rowsType][COLS_TYPE_CODE](x, nrow, ncol, rows, nrows, cols, ncols, qq, ans, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + return NULL; +} + + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + if (cores > 1 && nrows > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(qq) + sizeof(ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nrows); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, qq); + PUSH_ARGUMENT(buffer, ii, ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (nrows + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nrows) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nrows); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj; R_xlen_t *colOffset, rowIdx; X_C_TYPE *values; #ifdef ROWS_TYPE ROWS_C_TYPE *crows = (ROWS_C_TYPE*) rows; - // Check missing rows - for (ii=0; ii < nrows; ++ii) { - if (ROW_INDEX(crows,ii) == NA_R_XLEN_T) break; - } - if (ii < nrows && ncols > 0) { - error("Argument 'rows' must not contain missing value"); - } #endif #ifdef COLS_TYPE COLS_C_TYPE *ccols = (COLS_C_TYPE*) cols; - // Check missing cols - for (jj=0; jj < ncols; ++jj) { - if (COL_INDEX(ccols,jj) == NA_R_XLEN_T) break; - } - if (jj < ncols && nrows > 0) { - error("Argument 'cols' must not contain missing value"); - } #endif - /* R allocate memory for the 'values'. This will be - taken care of by the R garbage collector later on. */ - values = (X_C_TYPE *) R_alloc(ncols, sizeof(X_C_TYPE)); + /* C allocate memory for the 'values'. Should be freed manually. */ + values = Calloc(ncols, X_C_TYPE); /* Pre-calculate the column offsets */ - colOffset = (R_xlen_t *) R_alloc(ncols, sizeof(R_xlen_t)); + colOffset = Calloc(ncols, R_xlen_t); for (jj=0; jj < ncols; jj++) colOffset[jj] = COL_INDEX_NONA(ccols,jj) * nrow; @@ -78,11 +154,16 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { ans[ii] = values[qq]; } + + Free(values); + Free(colOffset); } /*************************************************************************** HISTORY: + 2015-08-09 [DJ] + o Pthread processing. 2015-07-08 [DJ] o Supported subsetted computation. 2014-11-06 [HB] diff --git a/src/rowRanges.c b/src/rowRanges.c index f0130ab9..f8415f9b 100644 --- a/src/rowRanges.c +++ b/src/rowRanges.c @@ -12,7 +12,7 @@ #define METHOD rowRanges #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int what, int narm, int hasna, X_C_TYPE *ans, int *is_counted +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int what, int narm, int hasna, X_C_TYPE *ans, R_xlen_t nans, int *is_counted, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -20,12 +20,12 @@ #include "templates-gen-matrix.h" -SEXP rowRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEXP hasNA) { +SEXP rowRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEXP hasNA, SEXP cores) { SEXP ans = NILSXP, ans2 = NILSXP; int *mins, *maxs; double *mins2, *maxs2; int *is_counted, all_counted = 0; - int what2, narm, hasna; + int what2, narm, hasna, cores2; R_xlen_t nrow, ncol, ii; /* Argument 'x' and 'dim': */ @@ -54,6 +54,15 @@ SEXP rowRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEX void *crows = validateIndices(rows, nrow, 0, &nrows, &rowsType); void *ccols = validateIndices(cols, ncol, 0, &ncols, &colsType); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + is_counted = (int *) R_alloc(nrows, sizeof(int)); if (isReal(x)) { @@ -62,7 +71,7 @@ SEXP rowRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEX } else { PROTECT(ans = allocVector(REALSXP, nrows)); } - rowRanges_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, what2, narm, hasna, REAL(ans), is_counted); + rowRanges_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, what2, narm, hasna, REAL(ans), nrows, is_counted, cores2); UNPROTECT(1); } else if (isInteger(x)) { if (what2 == 2) { @@ -70,7 +79,7 @@ SEXP rowRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEX } else { PROTECT(ans = allocVector(INTSXP, nrows)); } - rowRanges_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, what2, narm, hasna, INTEGER(ans), is_counted); + rowRanges_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, what2, narm, hasna, INTEGER(ans), nrows, is_counted, cores2); /* Any entries with zero non-missing values? */ all_counted = 1; @@ -138,6 +147,8 @@ SEXP rowRanges(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP what, SEXP naRm, SEX /*************************************************************************** HISTORY: + 2015-08-10 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-16 [HB] diff --git a/src/rowRanges_TYPE-template.h b/src/rowRanges_TYPE-template.h index d2740983..65cc9f99 100644 --- a/src/rowRanges_TYPE-template.h +++ b/src/rowRanges_TYPE-template.h @@ -3,10 +3,10 @@ void rowRanges_[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int what, int narm, int hasna, X_C_TYPE *ans, int *is_counted + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int what, int narm, int hasna, X_C_TYPE *ans, R_xlen_t nans, int *is_counted, int cores Arguments: - The following macros ("arguments") should be defined for the + The following macros ("arguments") should be defined for the template to work as intended. - METHOD_NAME: the name of the resulting function @@ -17,18 +17,118 @@ Henrik Bengtsson. Copyright: Henrik Bengtsson, 2014 - ***********************************************************************/ + ***********************************************************************/ #include +#include +#include #include "types.h" /* Expand arguments: X_TYPE => (X_C_TYPE, X_IN_C, [METHOD_NAME]) ANS_TYPE => (ANS_SXP, ANS_NA, ANS_C_TYPE, ANS_IN_C) */ -#include "templates-types.h" +#include "templates-types.h" + + +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, int, what); + POP_ARGUMENT(buffer0, ii, int, narm); + POP_ARGUMENT(buffer0, ii, int, hasna); + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, ans); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nans); + POP_ARGUMENT(buffer0, ii, int*, is_counted); + + ans += begin; + is_counted += begin; + + nrows = end - begin; + int rowsType; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; + rowsType = ROWS_TYPE_CODE; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + METHOD_NAME[rowsType][COLS_TYPE_CODE](x, nrow, ncol, rows, nrows, cols, ncols, what, narm, hasna, ans, nans, is_counted, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + return NULL; +} RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + if (cores > 1 && nrows > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(what) + sizeof(narm) + sizeof(hasna) + sizeof(ans) + sizeof(nans) + sizeof(is_counted); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nrows); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, what); + PUSH_ARGUMENT(buffer, ii, narm); + PUSH_ARGUMENT(buffer, ii, hasna); + PUSH_ARGUMENT(buffer, ii, ans); + PUSH_ARGUMENT(buffer, ii, nans); + PUSH_ARGUMENT(buffer, ii, is_counted); + + pthread_t threads[cores]; + + R_xlen_t gap = (nrows + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nrows) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nrows); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj; R_xlen_t colBegin, idx; X_C_TYPE value, *mins = NULL, *maxs = NULL; @@ -48,7 +148,7 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { narm = FALSE; if (hasna) { - skip = (int *) R_alloc(nrows, sizeof(int)); + skip = Calloc(nrows, int); for (ii=0; ii < nrows; ii++) { is_counted[ii] = 0; skip[ii] = 0; @@ -99,7 +199,7 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } else if (what == 1) { /* rowMaxs() */ maxs = ans; - + for (jj=0; jj < ncols; jj++) { colBegin = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); @@ -140,8 +240,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } else if (what == 2) { /* rowRanges() */ mins = ans; - maxs = &ans[nrows]; - + maxs = ans + nans; + for (jj=0; jj < ncols; jj++) { colBegin = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); @@ -185,17 +285,20 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } #endif } /* if (what ...) */ + + Free(skip); + } else { /* No missing values */ if (what == 0) { /* rowMins() */ mins = ans; - + /* Initiate results */ for (ii=0; ii < nrows; ii++) { mins[ii] = x[ii]; } - + for (jj=1; jj < ncols; jj++) { colBegin = COL_INDEX_NONA(ccols,jj) * nrow; for (ii=0; ii < nrows; ii++) { @@ -206,12 +309,12 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } else if (what == 1) { /* rowMax() */ maxs = ans; - + /* Initiate results */ for (ii=0; ii < nrows; ii++) { maxs[ii] = x[ii]; } - + for (jj=1; jj < ncols; jj++) { colBegin = COL_INDEX_NONA(ccols,jj) * nrow; for (ii=0; ii < nrows; ii++) { @@ -222,14 +325,14 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { } else if (what == 2) { /* rowRanges()*/ mins = ans; - maxs = &ans[nrows]; - + maxs = ans + nans; + /* Initiate results */ for (ii=0; ii < nrows; ii++) { mins[ii] = x[ii]; maxs[ii] = x[ii]; } - + for (jj=1; jj < ncols; jj++) { colBegin = COL_INDEX_NONA(ccols,jj) * nrow; for (ii=0; ii < nrows; ii++) { @@ -248,6 +351,8 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { /*************************************************************************** HISTORY: + 2015-08-10 [DJ] + o Pthread processing. 2015-06-07 [DJ] o Supported subsetted computation. 2014-11-16 [HB] diff --git a/src/rowRanksWithTies.c b/src/rowRanksWithTies.c index 4449d850..cf647e13 100644 --- a/src/rowRanksWithTies.c +++ b/src/rowRanksWithTies.c @@ -1,6 +1,6 @@ /*************************************************************************** Public methods: - SEXP rowRanksWithTies(SEXP x, SEXP rows, SEXP cols, SEXP tiesMethod, SEXP byRow) + SEXP rowRanksWithTies(SEXP x, SEXP rows, SEXP cols, SEXP tiesMethod, SEXP byRow, SEXP cores) Authors: Hector Corrada Bravo, Peter Langfelder and Henrik Bengtsson @@ -13,7 +13,7 @@ #define METHOD_TEMPLATE_H "rowRanksWithTies_TYPE_TIES-template.h" #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, int nrow, int ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, ANS_C_TYPE *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, ANS_C_TYPE *ans, R_xlen_t nrow_ans, int cores /***************************************************************** @@ -100,10 +100,10 @@ * tiesMethod: 1: maximum, 2: average, 3:minimum * The returned rank is a REAL matrix to accomodate average ranks */ -SEXP rowRanksWithTies(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP tiesMethod, SEXP byRow) { - int tiesmethod, byrow; +SEXP rowRanksWithTies(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP tiesMethod, SEXP byRow, SEXP cores) { + int tiesmethod, byrow, cores2; SEXP ans = NILSXP; - int nrow, ncol; + R_xlen_t nrow, ncol; /* Argument 'x' and 'dim': */ assertArgMatrix(x, dim, (R_TYPE_INT | R_TYPE_REAL), "x"); @@ -125,23 +125,32 @@ SEXP rowRanksWithTies(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP tiesMethod, S /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + /* Double matrices are more common to use. */ if (isReal(x)) { if (byrow) { switch (tiesmethod) { case 1: PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - rowRanks_tiesMax_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans)); + rowRanks_tiesMax_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans), nrows, cores2); UNPROTECT(1); break; case 2: PROTECT(ans = allocMatrix(REALSXP, nrows, ncols)); - rowRanks_tiesAverage_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, REAL(ans)); + rowRanks_tiesAverage_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, REAL(ans), nrows, cores2); UNPROTECT(1); break; case 3: PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - rowRanks_tiesMin_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans)); + rowRanks_tiesMin_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans), nrows, cores2); UNPROTECT(1); break; } /* switch */ @@ -149,17 +158,17 @@ SEXP rowRanksWithTies(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP tiesMethod, S switch (tiesmethod) { case 1: PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - colRanks_tiesMax_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans)); + colRanks_tiesMax_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans), nrows, cores2); UNPROTECT(1); break; case 2: PROTECT(ans = allocMatrix(REALSXP, nrows, ncols)); - colRanks_tiesAverage_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, REAL(ans)); + colRanks_tiesAverage_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, REAL(ans), nrows, cores2); UNPROTECT(1); break; case 3: PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - colRanks_tiesMin_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans)); + colRanks_tiesMin_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans), nrows, cores2); UNPROTECT(1); break; } /* switch */ @@ -169,17 +178,17 @@ SEXP rowRanksWithTies(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP tiesMethod, S switch (tiesmethod) { case 1: PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - rowRanks_tiesMax_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans)); + rowRanks_tiesMax_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans), nrows, cores2); UNPROTECT(1); break; case 2: PROTECT(ans = allocMatrix(REALSXP, nrows, ncols)); - rowRanks_tiesAverage_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, REAL(ans)); + rowRanks_tiesAverage_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, REAL(ans), nrows, cores2); UNPROTECT(1); break; case 3: PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - rowRanks_tiesMin_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans)); + rowRanks_tiesMin_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans), nrows, cores2); UNPROTECT(1); break; } /* switch */ @@ -187,17 +196,17 @@ SEXP rowRanksWithTies(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP tiesMethod, S switch (tiesmethod) { case 1: PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - colRanks_tiesMax_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans)); + colRanks_tiesMax_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans), nrows, cores2); UNPROTECT(1); break; case 2: PROTECT(ans = allocMatrix(REALSXP, nrows, ncols)); - colRanks_tiesAverage_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, REAL(ans)); + colRanks_tiesAverage_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, REAL(ans), nrows, cores2); UNPROTECT(1); break; case 3: PROTECT(ans = allocMatrix(INTSXP, nrows, ncols)); - colRanks_tiesMin_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans)); + colRanks_tiesMin_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, INTEGER(ans), nrows, cores2); UNPROTECT(1); break; } /* switch */ @@ -210,8 +219,10 @@ SEXP rowRanksWithTies(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP tiesMethod, S /*************************************************************************** HISTORY: + 2015-08-03 [DJ] + o Pthread processing. 2015-06-12 [DJ] o Supported subsetted computation. 2013-01-13 [HB] - o Added argument 'tiesMethod' to rowRanks(). + o Added argument 'tiesMethod' to rowRanks(). **************************************************************************/ diff --git a/src/rowRanksWithTies_TYPE_TIES-template.h b/src/rowRanksWithTies_TYPE_TIES-template.h index 10ad62cc..413e992d 100644 --- a/src/rowRanksWithTies_TYPE_TIES-template.h +++ b/src/rowRanksWithTies_TYPE_TIES-template.h @@ -3,10 +3,10 @@ Ranks_Real_ties[rowsType][colsType](ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, int nrow, int ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, ANS_C_TYPE *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, ANS_C_TYPE *ans, R_xlen_t nrow_ans, int cores Arguments: - The following macros ("arguments") should be defined for the + The following macros ("arguments") should be defined for the template to work as intended. - METHOD_NAME: the name of the resulting function @@ -21,6 +21,8 @@ Henrik Bengtsson [HB] ***********************************************************************/ #include +#include +#include #undef RANK #if TIESMETHOD == '0' /* min */ @@ -55,7 +57,117 @@ #endif -void METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { +/** Wrapper is used to call `METHOD_NAME`, as pthread worker. **/ +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE*, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, ANS_C_TYPE*, ans); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow_ans); + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + int rowsType = ROWS_TYPE_CODE, colsType = COLS_TYPE_CODE; + #if MARGIN == 'r' + ans += begin; + nrows = end - begin; + #ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; + #else + rows = indicesFromRange(begin, end, &rowsType); + #endif + + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, ans, nrow_ans, 0); + + #ifndef ROWS_TYPE + Free(rows); + #endif + + #elif MARGIN == 'c' + ans += begin * nrow_ans; + ncols = end - begin; + #ifdef COLS_TYPE + cols = (COLS_C_TYPE*) cols + begin; + #else + cols = indicesFromRange(begin, end, &colsType); + #endif + + METHOD_NAME[rowsType][colsType](x, nrow, ncol, rows, nrows, cols, ncols, ans, nrow_ans, 0); + + #ifndef COLS_TYPE + Free(cols); + #endif + #endif // MARGIN + + return NULL; +} + + +RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread +#if MARGIN == 'r' + R_xlen_t nv = nrows; +#elif MARGIN == 'c' + R_xlen_t nv = ncols; +#endif + if (cores > 1 && nv > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(ans) + sizeof(nrow_ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nv); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, ans); + PUSH_ARGUMENT(buffer, ii, nrow_ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (nv + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nv) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nv); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + ANS_C_TYPE rank; X_C_TYPE *values, current, tmp; R_xlen_t *colOffset; @@ -76,7 +188,7 @@ void METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { nVec = nrows; /* Pre-calculate the column offsets */ - colOffset = (R_xlen_t *) R_alloc(ncols, sizeof(R_xlen_t)); + colOffset = Calloc(ncols, R_xlen_t); for (jj=0; jj < ncols; jj++) colOffset[jj] = R_INDEX_OP(COL_INDEX(ccols,jj), *, nrow); @@ -85,13 +197,13 @@ void METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { nVec = ncols; /* Pre-calculate the column offsets */ - colOffset = (R_xlen_t *) R_alloc(nrows, sizeof(R_xlen_t)); + colOffset = Calloc(nrows, R_xlen_t); for (jj=0; jj < nrows; jj++) colOffset[jj] = ROW_INDEX(crows,jj); #endif - values = (X_C_TYPE *) R_alloc(nvalues, sizeof(X_C_TYPE)); - I = (int *) R_alloc(nvalues, sizeof(int)); + values = Calloc(nvalues, X_C_TYPE); + I = Calloc(nvalues, int); for (ii=0; ii < nVec; ii++) { #if MARGIN == 'r' @@ -130,26 +242,26 @@ void METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { /* Rprintf("Swapped vector:\n"); - for (jj=0; jj < nvalues; jj++) - { + for (jj=0; jj < nvalues; jj++) + { Rprintf(" %8.4f,", values[jj]); if (((jj+1) % 5==0) || (jj==nvalues-1)) Rprintf("\n"); } Rprintf("Index vector:\n"); - for (jj=0; jj 0) X_QSORT_I(values, I, 1, lastFinite + 1); - // Calculate the ranks. + // Calculate the ranks. for (jj=0; jj <= lastFinite;) { firstTie = jj; current = values[jj]; @@ -158,23 +270,28 @@ void METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { // Depending on rank method, get maximum, average, or minimum rank rank = RANK(firstTie, aboveTie); for (kk=firstTie; kk < aboveTie; kk++) { - ans[ ANS_INDEX_OF(I[kk], ii, nrows) ] = rank; + ans[ ANS_INDEX_OF(I[kk], ii, nrow_ans) ] = rank; } } // At this point jj = lastFinite + 1, no need to re-initialize again. for (; jj < nvalues; jj++) { - ans[ ANS_INDEX_OF(I[jj], ii, nrows) ] = ANS_NA; + ans[ ANS_INDEX_OF(I[jj], ii, nrow_ans) ] = ANS_NA; } // Rprintf("\n"); } -} + Free(colOffset); + Free(values); + Free(I); +} /*************************************************************************** HISTORY: + 2015-08-03 [DJ] + o Pthread processing. 2015-06-12 [DJ] o Supported subsetted computation. 2014-11-06 [HB] diff --git a/src/rowVars.c b/src/rowVars.c index b38405f4..ad971542 100644 --- a/src/rowVars.c +++ b/src/rowVars.c @@ -13,7 +13,7 @@ #define METHOD rowVars #define RETURN_TYPE void -#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int narm, int hasna, int byrow, double *ans +#define ARGUMENTS_LIST X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int narm, int hasna, int byrow, double *ans, int cores #define X_TYPE 'i' #include "templates-gen-matrix.h" @@ -21,11 +21,11 @@ #define X_TYPE 'r' #include "templates-gen-matrix.h" -#undef METHOD +#undef METHOD -SEXP rowVars(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP byRow) { - int narm, hasna, byrow; +SEXP rowVars(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP byRow, SEXP cores) { + int narm, hasna, byrow, cores2; SEXP ans; R_xlen_t nrow, ncol; @@ -49,6 +49,15 @@ SEXP rowVars(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP /* Argument 'byRow': */ byrow = asLogical(byRow); +#ifdef _USE_PTHREAD_ + /* Argument 'cores': */ + cores2 = asInteger(cores); + if (cores2 <= 0) + error("Argument 'cores' must be a positive value."); +#else + cores2 = 1; +#endif + if (!byrow) { SWAP(R_xlen_t, nrow, ncol); SWAP(void*, crows, ccols); @@ -62,9 +71,9 @@ SEXP rowVars(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP /* Double matrices are more common to use. */ if (isReal(x)) { - rowVars_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, narm, hasna, byrow, REAL(ans)); + rowVars_Real[rowsType][colsType](REAL(x), nrow, ncol, crows, nrows, ccols, ncols, narm, hasna, byrow, REAL(ans), cores2); } else if (isInteger(x)) { - rowVars_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, narm, hasna, byrow, REAL(ans)); + rowVars_Integer[rowsType][colsType](INTEGER(x), nrow, ncol, crows, nrows, ccols, ncols, narm, hasna, byrow, REAL(ans), cores2); } UNPROTECT(1); @@ -75,6 +84,8 @@ SEXP rowVars(SEXP x, SEXP dim, SEXP rows, SEXP cols, SEXP naRm, SEXP hasNA, SEXP /*************************************************************************** HISTORY: + 2015-08-02 [DJ] + o Pthread processing. 2015-06-13 [DJ] o Supported subsetted computation. 2014-11-18 [HB] diff --git a/src/rowVars_TYPE-template.h b/src/rowVars_TYPE-template.h index 3ed4ffce..0ec2df06 100644 --- a/src/rowVars_TYPE-template.h +++ b/src/rowVars_TYPE-template.h @@ -3,7 +3,7 @@ void rowVars_(ARGUMENTS_LIST) ARGUMENTS_LIST: - X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int narm, int hasna, int byrow, double *ans + X_C_TYPE *x, R_xlen_t nrow, R_xlen_t ncol, void *rows, R_xlen_t nrows, void *cols, R_xlen_t ncols, int narm, int hasna, int byrow, double *ans, int cores Arguments: The following macros ("arguments") should be defined for the @@ -20,6 +20,8 @@ ***********************************************************************/ #include #include +#include +#include #include "types.h" /* Expand arguments: @@ -28,7 +30,98 @@ #include "templates-types.h" +static void *WRAPPER_METHOD_NAME_ROWS_COLS(void *args) { + uint8_t *buffer = (uint8_t*) args; + int ii = 0; + POP_ARGUMENT(buffer, ii, uint8_t*, buffer0); + POP_ARGUMENT(buffer, ii, R_xlen_t, begin); + POP_ARGUMENT(buffer, ii, R_xlen_t, end); + + ii = 0; + POP_ARGUMENT(buffer0, ii, X_C_TYPE *, x); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrow); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncol); + POP_ARGUMENT(buffer0, ii, void*, rows); + POP_ARGUMENT(buffer0, ii, R_xlen_t, nrows); + POP_ARGUMENT(buffer0, ii, void*, cols); + POP_ARGUMENT(buffer0, ii, R_xlen_t, ncols); + POP_ARGUMENT(buffer0, ii, int, narm); + POP_ARGUMENT(buffer0, ii, int, hasna); + POP_ARGUMENT(buffer0, ii, int, byrow); + POP_ARGUMENT(buffer0, ii, double*, ans); + + extern RETURN_TYPE (*METHOD_NAME[3][3])(ARGUMENTS_LIST); + int rowsType = ROWS_TYPE_CODE; + + ans += begin; + nrows = end - begin; +#ifdef ROWS_TYPE + rows = (ROWS_C_TYPE*) rows + begin; +#else + rows = indicesFromRange(begin, end, &rowsType); +#endif + + METHOD_NAME[rowsType][COLS_TYPE_CODE](x, nrow, ncol, rows, nrows, cols, ncols, narm, hasna, byrow, ans, 0); + +#ifndef ROWS_TYPE + Free(rows); +#endif + return NULL; +} + + RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { + // Apply pthread + if (cores > 1 && nrows > 1) { + const static int memSize0 = sizeof(x) + sizeof(nrow) + sizeof(ncol) + + sizeof(rows) + sizeof(nrows) + sizeof(cols) + sizeof(ncols) + + sizeof(narm) + sizeof(hasna) + sizeof(byrow) + sizeof(ans); + R_xlen_t begin, end; + uint8_t *buffer0; + const static int memSize1 = sizeof(buffer0) + sizeof(begin) + sizeof(end); + + cores = MIN(cores, nrows); + uint8_t buffer[memSize0 + memSize1 * cores]; + buffer0 = buffer; + + int ii = 0; + PUSH_ARGUMENT(buffer, ii, x); + PUSH_ARGUMENT(buffer, ii, nrow); + PUSH_ARGUMENT(buffer, ii, ncol); + PUSH_ARGUMENT(buffer, ii, rows); + PUSH_ARGUMENT(buffer, ii, nrows); + PUSH_ARGUMENT(buffer, ii, cols); + PUSH_ARGUMENT(buffer, ii, ncols); + PUSH_ARGUMENT(buffer, ii, narm); + PUSH_ARGUMENT(buffer, ii, hasna); + PUSH_ARGUMENT(buffer, ii, byrow); + PUSH_ARGUMENT(buffer, ii, ans); + + pthread_t threads[cores]; + + R_xlen_t gap = (nrows + cores - 1) / cores; + int jj = 0; + begin = 0; + while (begin < nrows) { + uint8_t *args = buffer + ii; + end = MIN(begin + gap, nrows); + + PUSH_ARGUMENT(buffer, ii, buffer0); + PUSH_ARGUMENT(buffer, ii, begin); + PUSH_ARGUMENT(buffer, ii, end); + + pthread_create(threads+(jj++), NULL, &WRAPPER_METHOD_NAME_ROWS_COLS, args); + + begin = end; + } + + for (jj = 0; jj < cores; ++jj) { + pthread_join(threads[jj], NULL); + } + return; + } + + R_xlen_t ii, jj, kk, idx; R_xlen_t *colOffset; X_C_TYPE *values, value; @@ -41,16 +134,15 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { COLS_C_TYPE *ccols = (COLS_C_TYPE*) cols; #endif - /* R allocate memory for the 'values'. This will be - taken care of by the R garbage collector later on. */ - values = (X_C_TYPE *) R_alloc(ncols, sizeof(X_C_TYPE)); + /* C allocate memory for the 'values'. Should be freed manually. */ + values = Calloc(ncols, X_C_TYPE); /* If there are no missing values, don't try to remove them. */ if (hasna == FALSE) narm = FALSE; /* Pre-calculate the column offsets */ - colOffset = (R_xlen_t *) R_alloc(ncols, sizeof(R_xlen_t)); + colOffset = Calloc(ncols, R_xlen_t); if (byrow) { for (jj=0; jj < ncols; jj++) @@ -102,14 +194,19 @@ RETURN_TYPE METHOD_NAME_ROWS_COLS(ARGUMENTS_LIST) { ans[ii] = sigma2_d; } /* if (kk <= 1) */ - - R_CHECK_USER_INTERRUPT(ii); + // TODO: interrupt under pthread + // R_CHECK_USER_INTERRUPT(ii); } /* for (ii ...) */ + + Free(values); + Free(colOffset); } /*************************************************************************** HISTORY: + 2015-08-13 [DJ] + o Pthread processing. 2015-06-13 [DJ] o Supported subsetted computation. 2014-11-18 [HB] diff --git a/src/templates-types.h b/src/templates-types.h index e3fb8ad9..77c96e65 100644 --- a/src/templates-types.h +++ b/src/templates-types.h @@ -1,4 +1,5 @@ #include +#include "types.h" #include "macros.h" @@ -177,21 +178,25 @@ */ #undef ROW_INDEX_NONA #undef ROW_INDEX +#undef ROWS_TYPE_CODE #undef ROWS_C_TYPE #undef METHOD_NAME_ROWS #undef COL_INDEX_NONA #undef COL_INDEX +#undef COLS_TYPE_CODE #undef COLS_C_TYPE #undef METHOD_NAME_ROWS_COLS #ifdef ROWS_TYPE #define ROW_INDEX_NONA(rows, ii) ((R_xlen_t)rows[ii]-1) #if ROWS_TYPE == 'i' + #define ROWS_TYPE_CODE SUBSETTED_INTEGER #define ROWS_C_TYPE int #define ROW_INDEX(rows, ii) (rows[ii] == NA_INTEGER ? NA_R_XLEN_T : (R_xlen_t)rows[ii]-1) #define METHOD_NAME_ROWS CONCAT_MACROS(METHOD_NAME, intRows) #elif ROWS_TYPE == 'r' + #define ROWS_TYPE_CODE SUBSETTED_REAL #define ROWS_C_TYPE double #define ROW_INDEX(rows, ii) (ISNAN(rows[ii]) ? NA_R_XLEN_T : (R_xlen_t)rows[ii]-1) #define METHOD_NAME_ROWS CONCAT_MACROS(METHOD_NAME, realRows) @@ -201,6 +206,7 @@ #else #define ROW_INDEX_NONA(rows, ii) (ii) #define ROW_INDEX(rows, ii) (ii) + #define ROWS_TYPE_CODE SUBSETTED_ALL #define ROWS_C_TYPE void #define METHOD_NAME_ROWS CONCAT_MACROS(METHOD_NAME, noRows) #endif @@ -208,10 +214,12 @@ #ifdef COLS_TYPE #define COL_INDEX_NONA(cols, jj) ((R_xlen_t)cols[jj]-1) #if COLS_TYPE == 'i' + #define COLS_TYPE_CODE SUBSETTED_INTEGER #define COLS_C_TYPE int #define COL_INDEX(cols, jj) (cols[jj] == NA_INTEGER ? NA_R_XLEN_T : (R_xlen_t)cols[jj]-1) #define METHOD_NAME_ROWS_COLS CONCAT_MACROS(METHOD_NAME_ROWS, intCols) #elif COLS_TYPE == 'r' + #define COLS_TYPE_CODE SUBSETTED_REAL #define COLS_C_TYPE double #define COL_INDEX(cols, jj) (ISNAN(cols[jj]) ? NA_R_XLEN_T : (R_xlen_t)cols[jj]-1) #define METHOD_NAME_ROWS_COLS CONCAT_MACROS(METHOD_NAME_ROWS, realCols) @@ -221,6 +229,7 @@ #else #define COL_INDEX_NONA(cols, jj) (jj) #define COL_INDEX(cols, jj) (jj) + #define COLS_TYPE_CODE SUBSETTED_ALL #define COLS_C_TYPE void #define METHOD_NAME_ROWS_COLS CONCAT_MACROS(METHOD_NAME_ROWS, noCols) #endif @@ -252,21 +261,34 @@ #define METHOD_NAME_realRows_realCols CONCAT_MACROS(METHOD_NAME_realRows, realCols) +/* + WRAPPER used for pthread parameter passing + */ +#undef WRAPPER_METHOD_NAME_IDXS +#define WRAPPER_METHOD_NAME_IDXS CONCAT_MACROS(WRAPPER, METHOD_NAME_IDXS) + +#undef WRAPPER_METHOD_NAME_ROWS_COLS +#define WRAPPER_METHOD_NAME_ROWS_COLS CONCAT_MACROS(WRAPPER, METHOD_NAME_ROWS_COLS) + + /* Subsetted indexing: vector */ #undef IDX_INDEX_NONA #undef IDX_INDEX +#undef IDXS_TYPE_CODE #undef IDXS_C_TYPE #undef METHOD_NAME_IDXS #ifdef IDXS_TYPE #define IDX_INDEX_NONA(idxs, ii) ((R_xlen_t)idxs[ii]-1) #if IDXS_TYPE == 'i' + #define IDXS_TYPE_CODE SUBSETTED_INTEGER #define IDXS_C_TYPE int #define IDX_INDEX(idxs, ii) (idxs[ii] == NA_INTEGER ? NA_R_XLEN_T : (R_xlen_t)idxs[ii]-1) #define METHOD_NAME_IDXS CONCAT_MACROS(METHOD_NAME, intIdxs) #elif IDXS_TYPE == 'r' + #define IDXS_TYPE_CODE SUBSETTED_REAL #define IDXS_C_TYPE double #define IDX_INDEX(idxs, ii) (ISNAN(idxs[ii]) ? NA_R_XLEN_T : (R_xlen_t)idxs[ii]-1) #define METHOD_NAME_IDXS CONCAT_MACROS(METHOD_NAME, realIdxs) @@ -276,6 +298,7 @@ #else #define IDX_INDEX_NONA(idxs, ii) (ii) #define IDX_INDEX(idxs, ii) (ii) + #define IDXS_TYPE_CODE SUBSETTED_ALL #define IDXS_C_TYPE void #define METHOD_NAME_IDXS CONCAT_MACROS(METHOD_NAME, noIdxs) #endif diff --git a/src/types.h b/src/types.h index 8a63cde1..d96bb921 100644 --- a/src/types.h +++ b/src/types.h @@ -18,10 +18,9 @@ /* Subsetting index mode */ #ifndef SUBSETTED_MODE_INDEX #define SUBSETTED_MODE_INDEX -#define SUBSETTED_ALL 0 -#define SUBSETTED_INTEGER 1 -#define SUBSETTED_REAL 2 -#define SUBSETTED_NA 3 + #define SUBSETTED_ALL 0 + #define SUBSETTED_INTEGER 1 + #define SUBSETTED_REAL 2 #endif @@ -53,3 +52,12 @@ /* Macro to check for user interrupts every 2^20 iteration */ #define R_CHECK_USER_INTERRUPT(i) if (i % 1048576 == 0) R_CheckUserInterrupt() + +/* Detect OS */ +#if !defined(_WIN32) && (defined(__unix__) || defined(__unix) || (defined(__APPLE__) && defined(__MACH__))) +#define _UNIX_STYLE_OS_ +#endif + +#ifdef _UNIX_STYLE_OS_ +#define _USE_PTHREAD_ +#endif diff --git a/src/utils.h b/src/utils.h index 58fa447d..33068ba9 100644 --- a/src/utils.h +++ b/src/utils.h @@ -110,7 +110,13 @@ static R_INLINE R_xlen_t asR_xlen_t(SEXP x, R_xlen_t i) { /* Specified in validateIndices.c */ -void *validateIndices(SEXP idxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *type); +void *validateIndicesCheckNA(SEXP idxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *type, int *hasna); +void *indicesFromRange(R_xlen_t begin, R_xlen_t end, int *subsettedType); + +static R_INLINE void *validateIndices(SEXP idxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *type) { + int hasna; + return validateIndicesCheckNA(idxs, maxIdx, allowOutOfBound, ansNidxs, type, &hasna); +} static R_INLINE int IntegerFromReal(double x) { @@ -131,3 +137,15 @@ type tmp = x; \ x = y; \ y = tmp; \ } + +#define MIN(x, y) (x < y ? x : y) +#define MAX(x, y) (x < y ? y : x) + +#define PUSH_ARGUMENT(buffer, i, v) \ +memcpy(buffer+i, &v, sizeof(v)); \ +i += sizeof(v) + +#define POP_ARGUMENT(buffer, i, type, v) \ +type v; \ +memcpy(&v, buffer+i, sizeof(v)); \ +i += sizeof(v) diff --git a/src/validateIndices.c b/src/validateIndices.c index 7f5101c2..d5cb7745 100644 --- a/src/validateIndices.c +++ b/src/validateIndices.c @@ -36,10 +36,12 @@ for (ii = 0; ii < n; ++ ii) { \ /** idxs must not be NULL, which should be checked before calling this function. **/ -void* validateIndices_Logical(int *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType) { +void* validateIndices_Logical(int *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType, int *hasna) { R_xlen_t ii, jj, kk; R_xlen_t count1 = 0, count2 = 0; + // set default as no NA. + *hasna = FALSE; // set default type as SUBSETTED_INTEGER *subsettedType = SUBSETTED_INTEGER; if (nidxs == 0) { @@ -51,6 +53,7 @@ void* validateIndices_Logical(int *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int al if (!allowOutOfBound) { error("logical subscript too long"); } + *hasna = TRUE; // out-of-bound index is NA // count how many idx items for (ii = 0; ii < maxIdx; ++ ii) { @@ -114,6 +117,7 @@ void* validateIndices_Logical(int *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int al *subsettedType = SUBSETTED_ALL; return NULL; } + if (naCount) *hasna = TRUE; *ansNidxs = maxIdx / nidxs * count + count1; if (*subsettedType == SUBSETTED_INTEGER) { @@ -154,17 +158,20 @@ void* validateIndices_Logical(int *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int al * `allowOutOfBound` indicates whether to allow positve out of bound indexing. * `ansNidxs` is used for returning the new idxs array's length. * `subsettedType` is used for returning the new idxs array's datatype. + * `hasna` is TRUE, if NA is included in returned result. ************************************************************/ -void *validateIndices(SEXP idxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType) { +void *validateIndicesCheckNA(SEXP idxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType, int *hasna) { R_xlen_t nidxs = xlength(idxs); int mode = TYPEOF(idxs); + // Set no NA as default. + *hasna = FALSE; switch (mode) { case INTSXP: - return validateIndices_Integer(INTEGER(idxs), nidxs, maxIdx, allowOutOfBound, ansNidxs, subsettedType); + return validateIndices_Integer(INTEGER(idxs), nidxs, maxIdx, allowOutOfBound, ansNidxs, subsettedType, hasna); case REALSXP: - return validateIndices_Real(REAL(idxs), nidxs, maxIdx, allowOutOfBound, ansNidxs, subsettedType); + return validateIndices_Real(REAL(idxs), nidxs, maxIdx, allowOutOfBound, ansNidxs, subsettedType, hasna); case LGLSXP: - return validateIndices_Logical(LOGICAL(idxs), nidxs, maxIdx, allowOutOfBound, ansNidxs, subsettedType); + return validateIndices_Logical(LOGICAL(idxs), nidxs, maxIdx, allowOutOfBound, ansNidxs, subsettedType, hasna); case NILSXP: *subsettedType = SUBSETTED_ALL; *ansNidxs = maxIdx; @@ -192,16 +199,18 @@ SEXP validate(SEXP idxs, SEXP maxIdx, SEXP allowOutOfBound) { int callowOutOfBound = asLogicalNoNA(allowOutOfBound, "allowOutOfBound"); void *cidxs; + // Set no NA as default. + int hasna = FALSE; int mode = TYPEOF(idxs); switch (mode) { case INTSXP: - cidxs = validateIndices_Integer(INTEGER(idxs), nidxs, cmaxIdx, callowOutOfBound, &ansNidxs, &subsettedType); + cidxs = validateIndices_Integer(INTEGER(idxs), nidxs, cmaxIdx, callowOutOfBound, &ansNidxs, &subsettedType, &hasna); break; case REALSXP: - cidxs = validateIndices_Real(REAL(idxs), nidxs, cmaxIdx, callowOutOfBound, &ansNidxs, &subsettedType); + cidxs = validateIndices_Real(REAL(idxs), nidxs, cmaxIdx, callowOutOfBound, &ansNidxs, &subsettedType, &hasna); break; case LGLSXP: - cidxs = validateIndices_Logical(LOGICAL(idxs), nidxs, cmaxIdx, callowOutOfBound, &ansNidxs, &subsettedType); + cidxs = validateIndices_Logical(LOGICAL(idxs), nidxs, cmaxIdx, callowOutOfBound, &ansNidxs, &subsettedType, &hasna); break; case NILSXP: return R_NilValue; @@ -224,3 +233,25 @@ SEXP validate(SEXP idxs, SEXP maxIdx, SEXP allowOutOfBound) { UNPROTECT(1); return ans; } + + +/************************************************************* + * Used to generate indices begin+1:end. + * `subsettedType` is used for returning the new idxs array's datatype. + ************************************************************/ +void *indicesFromRange(R_xlen_t begin, R_xlen_t end, int *subsettedType) { + R_xlen_t ii, kk; + + ii = 0; + if (end > R_INT_MAX) { + double *ans = Calloc(end - begin, double); + for (kk = begin + 1; kk <= end; ++ kk) ans[ii ++] = kk; + *subsettedType = SUBSETTED_REAL; + return ans; + } + + int *ans = Calloc(end - begin, int); + for (kk = begin + 1; kk <= end; ++ kk) ans[ii ++] = kk; + *subsettedType = SUBSETTED_INTEGER; + return ans; +} diff --git a/src/validateIndices_TYPE-template.h b/src/validateIndices_TYPE-template.h index 518e2b6e..6aa2522f 100644 --- a/src/validateIndices_TYPE-template.h +++ b/src/validateIndices_TYPE-template.h @@ -1,10 +1,6 @@ /*********************************************************************** TEMPLATE: - void validateIndices_[ROWS_TYPE][COLS_TYPE](X_C_TYPE *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType) - - GENERATES: - void validateIndices_Real[ROWS_TYPE][COLS_TYPE](double *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType) - void validateIndices_Integer[ROWS_TYPE][COLS_TYPE](int *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType) + void validateIndices_(X_C_TYPE *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType, int *hasna) Arguments: The following macros ("arguments") should be defined for the @@ -51,7 +47,9 @@ static R_INLINE int RealFromIndex_TYPE(X_C_TYPE x, R_xlen_t maxIdx) { } /** idxs must not be NULL, which should be checked before calling this function. **/ -void* METHOD_NAME(X_C_TYPE *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType) { +void* METHOD_NAME(X_C_TYPE *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutOfBound, R_xlen_t *ansNidxs, int *subsettedType, int *hasna) { + // set default as no NA. + *hasna = FALSE; // For a un-full positive legal idxs array, we should use SUBSETTED_INTEGER as default. *subsettedType = SUBSETTED_INTEGER; @@ -78,12 +76,15 @@ void* METHOD_NAME(X_C_TYPE *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutO if (!X_ISNAN(idx)) { if (idx > maxIdx) { if (!allowOutOfBound) error("subscript out of bounds"); + *hasna = TRUE; // out-of-bound index is NA needReAlloc = TRUE; } #if X_TYPE == 'r' if (idx > R_INT_MAX) *subsettedType = SUBSETTED_REAL; #endif - }// else error("NA index is not supported"); // NOTE: currently + } else { + *hasna = TRUE; + } state = 1; ++ count; @@ -107,18 +108,10 @@ void* METHOD_NAME(X_C_TYPE *idxs, R_xlen_t nidxs, R_xlen_t maxIdx, int allowOutO if (state >= 0) { if (*subsettedType == SUBSETTED_INTEGER) { // NOTE: braces is needed here, because of macro-defined function -//#if X_TYPE == 'i' RETURN_VALIDATED_ANS(int, nidxs, idxs[ii], IntegerFromIndex_TYPE(idxs[ii],maxIdx),); -//#else // X_TYPE == 'r' -// RETURN_VALIDATED_ANS(int, nidxs, idxs[ii], idxs[ii] > maxIdx || IS_INF(idxs[ii]) ? NA_INTEGER : IntegerFromReal(idxs[ii]),); -//#endif } // *subsettedType == SUBSETTED_REAL -//#if X_TYPE == 'i' RETURN_VALIDATED_ANS(double, nidxs, idxs[ii], RealFromIndex_TYPE(idxs[ii],maxIdx),); -//#else // X_TYPE == 'r' -// RETURN_VALIDATED_ANS(double, nidxs, idxs[ii], idxs[ii] > maxIdx ? NA_REAL : idxs[ii],); -//#endif } // state < 0 diff --git a/tests/rowDiffs_parallel.R b/tests/rowDiffs_parallel.R index 5f6b6306..76c004fc 100644 --- a/tests/rowDiffs_parallel.R +++ b/tests/rowDiffs_parallel.R @@ -21,12 +21,13 @@ source("utils/validateIndicesFramework.R") x <- matrix(runif(6*6, min=-6, max=6), nrow=6, ncol=6) storage.mode(x) <- "integer" lag <- 2L -differences <- 1L for (rows in indexCases) { for (cols in indexCases) { - validateIndicesTestMatrix(x, rows, cols, ftest=rowDiffs, fsure=rowDiffs_R, lag=lag, differences=differences, mc.cores=2L) - validateIndicesTestMatrix(x, rows, cols, ftest=function(x, rows, cols, ...) { - t(colDiffs(t(x), rows=cols, cols=rows, ...)) - }, fsure=rowDiffs_R, lag=lag, differences=differences, mc.cores=2L) + for (differences in 1:3) { + validateIndicesTestMatrix(x, rows, cols, ftest=rowDiffs, fsure=rowDiffs_R, lag=lag, differences=differences, mc.cores=2L) + validateIndicesTestMatrix(x, rows, cols, ftest=function(x, rows, cols, ...) { + t(colDiffs(t(x), rows=cols, cols=rows, ...)) + }, fsure=rowDiffs_R, lag=lag, differences=differences, mc.cores=2L) + } } } diff --git a/tests/rowDiffs_subset.R b/tests/rowDiffs_subset.R index ffafd653..9310d4ef 100644 --- a/tests/rowDiffs_subset.R +++ b/tests/rowDiffs_subset.R @@ -23,7 +23,7 @@ storage.mode(x) <- "integer" for (rows in indexCases) { for (cols in indexCases) { for (lag in 1:2) { - for (differences in 1:2) { + for (differences in 1:3) { validateIndicesTestMatrix(x, rows, cols, ftest=rowDiffs, fsure=rowDiffs_R, lag=lag, differences=differences) validateIndicesTestMatrix(x, rows, cols, ftest=function(x, rows, cols, ...) { t(colDiffs(t(x), rows=cols, cols=rows, ...)) diff --git a/tests/rowMads_parallel.R b/tests/rowMads_parallel.R index cc8ec29d..c9ac7b3d 100644 --- a/tests/rowMads_parallel.R +++ b/tests/rowMads_parallel.R @@ -15,7 +15,7 @@ x <- matrix(runif(6*6, min=-6, max=6), nrow=6, ncol=6) storage.mode(x) <- "integer" for (rows in indexCases) { for (cols in indexCases) { - validateIndicesTestMatrix(x, rows, cols, ftest=rowMads, fsure=rowMads_R) - validateIndicesTestMatrix(x, rows, cols, fcolTest=colMads, fsure=rowMads_R) + validateIndicesTestMatrix(x, rows, cols, ftest=rowMads, fsure=rowMads_R, mc.cores=2L) + validateIndicesTestMatrix(x, rows, cols, fcolTest=colMads, fsure=rowMads_R, mc.cores=2L) } } diff --git a/tests/rowRanks_parallel.R b/tests/rowRanks_parallel.R index cd13afa2..3a3811be 100644 --- a/tests/rowRanks_parallel.R +++ b/tests/rowRanks_parallel.R @@ -15,10 +15,10 @@ x <- matrix(runif(6*6, min=-6, max=6), nrow=6, ncol=6) storage.mode(x) <- "integer" for (rows in indexCases) { for (cols in indexCases) { - validateIndicesTestMatrix(x, rows, cols, ftest=rowRanks, fsure=rowRanks_R, ties.method="average", mc.cores=2L) + validateIndicesTestMatrix(x, rows, cols, ftest=rowRanks, fsure=rowRanks_R, ties.method="average", mc.cores=2L, debug=TRUE) validateIndicesTestMatrix(x, rows, cols, ftest=function(x, rows, cols, ...) { t(colRanks(t(x), rows=cols, cols=rows, preserveShape=TRUE, ...)) - }, fsure=rowRanks_R, ties.method="average", mc.cores=2L) + }, fsure=rowRanks_R, ties.method="average", mc.cores=2L, debug=TRUE) } }