0

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.

Dharman
  • 30,962
  • 25
  • 85
  • 135
mshaffer
  • 959
  • 1
  • 9
  • 19

0 Answers0