Template abstraction of arma::mat and arma::sp_mat

[GabrielHoffman]

I am writing code using C++/armadillo, and then writing an interface to R with RcppArmadillo.

In this simple example, I want to multiply two matrices and I can use templates to handle cases where the matrices are any combination of mat and sp_mat. This works fine at the C++/armadillo level:

template <class T1, class T2> arma::mat h_tt(T1 &M1, T2 &M2) {
  return M1 * M2;
}

But I would like to export a function to R that handles any combination of mat and sp_mat. I could to this "manually", and write an exported function h_*() for each combination of mat and sp_mat. This works:

library(Rcpp)
library(Matrix)

# matrix, matrix
cppFunction("
arma::mat h_mm(const arma::mat& M1, const arma::mat& M2){
    return( M1 * M2 );
}
", depends="RcppArmadillo")

# matrix, sparseMatrix
cppFunction("
arma::mat h_ms(const arma::mat& M1, const arma::sp_mat& M2){
    return( M1 * M2 );
}", depends="RcppArmadillo")

A = matrix(1:4, 2,2)
B = matrix(5:8, 2,2)

# matrix, matrix
h_mm( A, B)
#     [,1] [,2]
#[1,]   23   31
#[2,]   34   46

# matrix, sparseMatrix
h_ms( A, as(B, "sparseMatrix"))
#     [,1] [,2]
#[1,]   23   31
#[2,]   34   46

But it quickly becomes unmanageable if I have multiple arguments that can be mat or sp_mat, and of course my real application is more complicated that multiplying two matrices.

I suspect there is a way to do this by having the arguments be of type SEXP or S4, but I can't figure it out.

Any thoughts?

Best, Gabriel

[eddelbuettel]

Quick thought: Yes.

Sparse matrices are S4 objects, dense matrices are straight REAL (ie NumericMatrix in Rcpp). You would need to run-time switch around the SEXP. That is the way it is as only talks in terms of SEXP objects both in and out.

This whole nexus is hard and AFAIK one of the reasons Matrix is in S4 to allow the multitude of dense + sparse matrix and vector ops.

[eddelbuettel]

Did you make any progress? Would it be ok with you to close this as it seems a little out of scope.

[GabrielHoffman]

I couldn't figure this out. Before you close, do you have a suggestion about the right place for this question? But I figure if you didn't know, its unlikely someone else would

Oh well, thanks for your help, Gabriel

[eddelbuettel]

I think it is a loaded / complicated topic and no I don't have a great alternative venue for it either. :crying_cat_face: