initial headers for dense Matrix (not matrix)

Author: pachadotdev

cpp11test/R/cpp11.R

@@ -108,6 +108,18 @@ 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)
+}
+
+row_sums_4 <- function(x) {
+  .Call(`_cpp11test_row_sums_4`, x)
+}
+
 col_sums <- function(x) {
   .Call(`_cpp11test_col_sums`, x)
 }

cpp11test/src/cpp11.cpp

@@ -209,6 +209,27 @@ 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<cpp11::decay_t<cpp11::dge_matrix<>>>(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<cpp11::decay_t<cpp11::dsy_matrix<>>>(x)));
+  END_CPP11
+}
+// matrix.cpp
+cpp11::writable::doubles row_sums_4(cpp11::dsp_matrix<> x);
+extern "C" SEXP _cpp11test_row_sums_4(SEXP x) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(row_sums_4(cpp11::as_cpp<cpp11::decay_t<cpp11::dsp_matrix<>>>(x)));
+  END_CPP11
+}
+// matrix.cpp
 cpp11::doubles col_sums(cpp11::doubles_matrix<cpp11::by_column> x);
 extern "C" SEXP _cpp11test_col_sums(SEXP x) {
   BEGIN_CPP11
@@ -518,6 +539,9 @@ 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_row_sums_4",               (DL_FUNC) &_cpp11test_row_sums_4,               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},

cpp11test/src/matrix.cpp

@@ -1,4 +1,5 @@
 #include "cpp11/matrix.hpp"
+#include "cpp11/dmatrix.hpp"
 #include "Rmath.h"
 #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<cpp11::by_column> x) {
   cpp11::writable::doubles sums(x.ncol());
 

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", {

inst/include/cpp11.hpp

@@ -14,6 +14,7 @@
 #include "cpp11/list_of.hpp"
 #include "cpp11/logicals.hpp"
 #include "cpp11/matrix.hpp"
+#include "cpp11/dmatrix.hpp"
 #include "cpp11/named_arg.hpp"
 #include "cpp11/protect.hpp"
 #include "cpp11/r_bool.hpp"

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 <typename S = by_column>
+class dge_matrix {
+ private:
+  writable::r_vector<double> 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<double>(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 <typename S = by_column>
+using dge_matrix = cpp11::dge_matrix<S>;
+}  // namespace writable
+
+}  // namespace cpp11
+
+// dsyMatrix: Symmetric Dense (Packed or Unpacked) Numeric Matrices
+
+namespace cpp11 {
+
+template <typename S = by_column>
+class dsy_matrix {
+ private:
+  writable::r_vector<double> 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<double>(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 <typename S = by_column>
+using dsy_matrix = cpp11::dsy_matrix<S>;
+}  // namespace writable
+
+}  // namespace cpp11