diff --git a/cpp11test/R/cpp11.R b/cpp11test/R/cpp11.R index 038e7b76..29775d0d 100644 --- a/cpp11test/R/cpp11.R +++ b/cpp11test/R/cpp11.R @@ -108,6 +108,14 @@ row_sums <- function(x) { .Call(`_cpp11test_row_sums`, x) } +row_sums_2 <- function(x) { + .Call(`_cpp11test_row_sums_2`, x) +} + +row_sums_3 <- function(x) { + .Call(`_cpp11test_row_sums_3`, x) +} + col_sums <- function(x) { .Call(`_cpp11test_col_sums`, x) } diff --git a/cpp11test/src/cpp11.cpp b/cpp11test/src/cpp11.cpp index 421de637..c9ef3e1b 100644 --- a/cpp11test/src/cpp11.cpp +++ b/cpp11test/src/cpp11.cpp @@ -209,6 +209,20 @@ extern "C" SEXP _cpp11test_row_sums(SEXP x) { END_CPP11 } // matrix.cpp +cpp11::doubles row_sums_2(cpp11::dge_matrix<> x); +extern "C" SEXP _cpp11test_row_sums_2(SEXP x) { + BEGIN_CPP11 + return cpp11::as_sexp(row_sums_2(cpp11::as_cpp>>(x))); + END_CPP11 +} +// matrix.cpp +cpp11::doubles row_sums_3(cpp11::dsy_matrix<> x); +extern "C" SEXP _cpp11test_row_sums_3(SEXP x) { + BEGIN_CPP11 + return cpp11::as_sexp(row_sums_3(cpp11::as_cpp>>(x))); + END_CPP11 +} +// matrix.cpp cpp11::doubles col_sums(cpp11::doubles_matrix x); extern "C" SEXP _cpp11test_col_sums(SEXP x) { BEGIN_CPP11 @@ -518,6 +532,8 @@ static const R_CallMethodDef CallEntries[] = { {"_cpp11test_rcpp_sum_int_for_", (DL_FUNC) &_cpp11test_rcpp_sum_int_for_, 1}, {"_cpp11test_remove_altrep", (DL_FUNC) &_cpp11test_remove_altrep, 1}, {"_cpp11test_row_sums", (DL_FUNC) &_cpp11test_row_sums, 1}, + {"_cpp11test_row_sums_2", (DL_FUNC) &_cpp11test_row_sums_2, 1}, + {"_cpp11test_row_sums_3", (DL_FUNC) &_cpp11test_row_sums_3, 1}, {"_cpp11test_string_proxy_assignment_", (DL_FUNC) &_cpp11test_string_proxy_assignment_, 0}, {"_cpp11test_string_push_back_", (DL_FUNC) &_cpp11test_string_push_back_, 0}, {"_cpp11test_sum_dbl_accumulate2_", (DL_FUNC) &_cpp11test_sum_dbl_accumulate2_, 1}, diff --git a/cpp11test/src/matrix.cpp b/cpp11test/src/matrix.cpp index 10348945..6ff631e6 100644 --- a/cpp11test/src/matrix.cpp +++ b/cpp11test/src/matrix.cpp @@ -1,5 +1,6 @@ #include "cpp11/matrix.hpp" #include "Rmath.h" +#include "cpp11/dmatrix.hpp" #include "cpp11/doubles.hpp" using namespace cpp11; @@ -82,10 +83,42 @@ using namespace Rcpp; } ++i; } - + // 4 return sums; } +[[cpp11::register]] cpp11::doubles row_sums_2(cpp11::dge_matrix<> x) { + int nrow = x.nrow(); + int ncol = x.ncol(); + cpp11::writable::doubles result(nrow); + + for (int i = 0; i < nrow; ++i) { + double sum = 0; + for (int j = 0; j < ncol; ++j) { + sum += x(i, j); + } + result[i] = sum; + } + + return result; +} + +[[cpp11::register]] cpp11::doubles row_sums_3(cpp11::dsy_matrix<> x) { + int nrow = x.nrow(); + int ncol = x.ncol(); + cpp11::writable::doubles result(nrow); + + for (int i = 0; i < nrow; ++i) { + double sum = 0; + for (int j = 0; j < ncol; ++j) { + sum += x(i, j); + } + result[i] = sum; + } + + return result; +} + [[cpp11::register]] cpp11::doubles col_sums(cpp11::doubles_matrix x) { cpp11::writable::doubles sums(x.ncol()); diff --git a/cpp11test/tests/testthat/test-matrix.R b/cpp11test/tests/testthat/test-matrix.R index a43f59b0..0dbf09f6 100644 --- a/cpp11test/tests/testthat/test-matrix.R +++ b/cpp11test/tests/testthat/test-matrix.R @@ -9,6 +9,16 @@ test_that("row_sums gives same result as rowSums", { y <- cbind(x1 = 3, x2 = c(4:1, 2:5)) y[3, ] <- NA; expect_equal(row_sums(y), rowSums(y)) + + # Pacha: test dgeMatrix and other d**Matrix (Matrix package) + + z <- cbind(x1 = 3, x2 = c(4, 4)) + z2 <- Matrix::Matrix(z) # dgeMatrix + expect_equal(rowSums(z), row_sums_2(z2)) + + z3 <- matrix(c(1, 2, 2, 1), nrow = 2) + z4 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 1)) + expect_equal(rowSums(z3), row_sums_3(z4)) }) test_that("col_sums gives same result as colSums", { diff --git a/inst/include/cpp11.hpp b/inst/include/cpp11.hpp index 71e1cf1d..a41a0098 100644 --- a/inst/include/cpp11.hpp +++ b/inst/include/cpp11.hpp @@ -5,6 +5,7 @@ #include "cpp11/as.hpp" #include "cpp11/attribute_proxy.hpp" #include "cpp11/data_frame.hpp" +#include "cpp11/dmatrix.hpp" #include "cpp11/doubles.hpp" #include "cpp11/environment.hpp" #include "cpp11/external_pointer.hpp" diff --git a/inst/include/cpp11/dmatrix.hpp b/inst/include/cpp11/dmatrix.hpp new file mode 100644 index 00000000..7534cd07 --- /dev/null +++ b/inst/include/cpp11/dmatrix.hpp @@ -0,0 +1,133 @@ +#pragma once + +#include "cpp11/R.hpp" +#include "cpp11/r_vector.hpp" +#include "cpp11/sexp.hpp" + +// dgeMatrix: Class "dgeMatrix" of Dense Numeric (S4 Class) Matrices + +namespace cpp11 { + +template +class dge_matrix { + private: + writable::r_vector vector_; + int nrow_; + int ncol_; + + public: + dge_matrix(SEXP data) { + SEXP x_slot = R_do_slot(data, Rf_install("x")); + + // get dimensions + + // nrow_ = Rf_nrows(data); + // ncol_ = Rf_ncols(data); + + // use attr instead (Matrix != matrix class) + SEXP dim = R_do_slot(data, Rf_install("Dim")); + nrow_ = INTEGER(dim)[0]; + ncol_ = INTEGER(dim)[1]; + + if (nrow_ <= 0 || ncol_ <= 0) { + throw std::runtime_error("Invalid matrix dimensions"); + } + + vector_ = writable::r_vector(Rf_allocVector(REALSXP, Rf_length(x_slot))); + for (int i = 0; i < Rf_length(x_slot); ++i) { + vector_[i] = REAL(x_slot)[i]; + } + + SEXP dims = PROTECT(Rf_allocVector(INTSXP, 2)); + INTEGER(dims)[0] = nrow_; + INTEGER(dims)[1] = ncol_; + vector_.attr(R_DimSymbol) = dims; + UNPROTECT(1); + } + + int nrow() const { return nrow_; } + int ncol() const { return ncol_; } + + double operator()(int row, int col) const { return vector_[row + (col * nrow_)]; } +}; + +namespace writable { +template +using dge_matrix = cpp11::dge_matrix; +} // namespace writable + +} // namespace cpp11 + +// dsyMatrix: Symmetric Dense (Packed or Unpacked) Numeric Matrices + +namespace cpp11 { + +template +class dsy_matrix { + private: + writable::r_vector vector_; + int nrow_; + int ncol_; + bool is_packed_; + bool is_upper_; + + public: + dsy_matrix(SEXP data) { + SEXP x_slot = R_do_slot(data, Rf_install("x")); + + // use attr instead (Matrix != matrix class) + SEXP dim = R_do_slot(data, Rf_install("Dim")); + nrow_ = INTEGER(dim)[0]; + ncol_ = INTEGER(dim)[1]; + + if (nrow_ <= 0 || ncol_ <= 0) { + throw std::runtime_error("Invalid matrix dimensions"); + } + + // Check if the matrix is packed and if it's upper or lower triangular + SEXP uplo = R_do_slot(data, Rf_install("uplo")); + is_upper_ = (Rf_asChar(uplo) == Rf_mkChar("U")); + is_packed_ = (Rf_length(x_slot) == (nrow_ * (nrow_ + 1)) / 2); + + vector_ = writable::r_vector(Rf_allocVector(REALSXP, Rf_length(x_slot))); + for (int i = 0; i < Rf_length(x_slot); ++i) { + vector_[i] = REAL(x_slot)[i]; + } + + SEXP dims = PROTECT(Rf_allocVector(INTSXP, 2)); + INTEGER(dims)[0] = nrow_; + INTEGER(dims)[1] = ncol_; + vector_.attr(R_DimSymbol) = dims; + UNPROTECT(1); + } + + int nrow() const { return nrow_; } + int ncol() const { return ncol_; } + + double operator()(int row, int col) const { + if (is_packed_) { + if (is_upper_) { + if (row <= col) { + return vector_[col * (col + 1) / 2 + row]; + } else { + return vector_[row * (row + 1) / 2 + col]; + } + } else { + if (row >= col) { + return vector_[row * (row + 1) / 2 + col]; + } else { + return vector_[col * (col + 1) / 2 + row]; + } + } + } else { + return vector_[row + (col * nrow_)]; + } + } +}; + +namespace writable { +template +using dsy_matrix = cpp11::dsy_matrix; +} // namespace writable + +} // namespace cpp11