Search code examples
c++eigenelementwise-operations

Error using Eigen: Perform element-wise multiplication between a vector and matrix


I am trying to perform an element-wise multiplication of a row vector with matrix. In MATLAB this would be simply done by the "dot" operator or:

  deriv  = 1i * k .* fk;

where k is row vector and fk is a matrix. Now in C++ I have this code:

static const int nx = 10;
static const int ny = 10; 
static const int nyk = ny/2 + 1;
static const int nxk = nx/2 + 1;
static const int ncomp = 2;

Matrix <double, 1, nx> eK; 
eK.setZero();

for(int i = 0; i < nx; i++){
            eK[i] = //some expression
    }
fftw_complex *UOut; 
    UOut= (fftw_complex*) fftw_malloc((((nx)*(ny+1))*nyk)* sizeof(fftw_complex));

for (int i = 0; i < nx; i++){
        for (int j = 0; j < ny+1; j++){ 
            for (int k = 0; k < ncomp; k++){
                UOut[i*(ny+1)+j][k] = //FFT of some expression 
            }
        }
    }

Eigen::Map<Eigen::MatrixXcd, Eigen::Unaligned> U(reinterpret_cast<std::complex<double>*>(UOut),(ny+1),nx); 

Now, I am trying to take the product of eK which is a row vector of 1 x 10 and the matrix U of a 11 x 10. I tried few things, none of which seem to work really:

U = 1i * eKX.array() * euhX.array() ; //ERROR
static assertion failed: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES

( \
      |      ~~~
  176 |         (int(Eigen::internal::size_of_xpr_at_compile_time<TYPE0>::ret)==0 && int(Eigen::internal::size_of_xpr_at_compile_time<TYPE1>::ret)==0) \
      |         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  177 |     || (\
      |     ^~~~~
  178 |           (int(TYPE0::RowsAtCompileTime)==Eigen::Dynamic \

Solution

  • Eigen doesn't do broadcasting the same way Matlab or Numpy do unless you explicitely ask for it, for example with matrix.array().rowwise() * vector.array()

    The IMHO clearer form is to interpret the vector as a diagonal matrix.

    Eigen::VectorXd eK = ...;
    Eigen::Map<Eigen::MatrixXcd, Eigen::Unaligned> U = ...;
    
    Eigen::MatrixXcd result = U * (eK * 1i).asDiagonal();