I am using Eigen
. For a real matrix, things work as expected:
// [[Rcpp::export]]
VectorXd matrix_diagonal(Map<MatrixXd> M)
{
VectorXd RES = M.diagonal();
return RES;
}
For a complex matrix, I am able to get it to work correctly if I pass in the elements Re(m)
and Im(m)
separately:
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]
// should have COMPLEX already ... #include <complex>
using Eigen::Map; // 'maps' rather than copies
using Eigen::MatrixXd; // variable size matrix, double precision
using Eigen::MatrixXcd; // ABOVE, complex
using Eigen::VectorXd; // variable size vector, double precision
using Eigen::VectorXcd; // ABOVE, complex
MatrixXcd toComplex_(const MatrixXd& Re, const MatrixXd& Im)
{
const std::complex<double> I_ {0.0, 1.0};
return Re.cast<std::complex<double>>() +
I_ * Im.cast<std::complex<double>>();
}
// [[Rcpp::export]]
VectorXd matrix_diagonal(Map<MatrixXd> M)
{
VectorXd RES = M.diagonal();
return RES;
}
// [[Rcpp::export]]
VectorXcd matrix_diagonal_complex(const MatrixXd& Re, const MatrixXd& Im)
{
MatrixXcd M = toComplex_(Re, Im);
VectorXcd RES = M.diagonal();
return RES;
}
using a toComplex_
function I found on the forums. This works as expected.
The forums talk about a ComplexMatrix class in Rccp. I am assuming it wraps the MatrixXcd
, but maybe not a Map
?
MAP of complex, ideal solution
// [[Rcpp::export]]
VectorXcd matrix_diagonal_complex2(Map<MatrixXcd> M)
{
// MatrixXcd M = toComplex_(Re, Im);
VectorXcd RES = M.diagonal();
return RES;
}
/*
// ERROR snippet
include/RcppEigen.h:25,
include/Eigen/src/Core/arch/SSE/PacketMath.h:60:39: warning: ignoring attributes on template argument '__m128' [-Wignored-attributes]
60 | template<> struct is_arithmetic<__m128> { enum { value = true }; };
*/
The ideal solution uses the Map
and MatrixXcd
syntax of eigen
. I thought Rcpp
would wrap the result for me. I don't understand if the error is the eigen
MAP or something else.
Question: Using Rcpp
, how to define a eigen
complex map matrix Map<MatrixXcd> M
m = as.matrix(structure(c(1, 0, 4, 0, 3, 0, 2, 0, 5), .Dim = c(3L, 3L)));
mc = m + (3+2i);
# as expected
matrix_diagonal(m);
# [1] 1 3 5
# as expected, throws error regarding type
matrix_diagonal(mc);
# Error in matrix_diagonal(mc) : Wrong R type for mapped matrix
# In addition: Warning message:
# In matrix_diagonal(mc) : imaginary parts discarded in coercion
# as expected with prep work( separating into buckets )
matrix_diagonal_complex(Re(mc), Im(mc));
# [1] 4+2i 6+2i 8+2i
# OBJECTIVE (direct call, preferably as MAP)
matrix_diagonal_complex2(mc);
NOTE: The const MatrixXcd& M
version does work. Not the MAP. Hence the question.