Python Numpy methods correspond to C++ Eigen make crash

I have a NumPy arithmetic needs to translate C++ with Eigen.

# D is a 20001x13 matrix stacked from Dva and Dvb, then multiply by w_f.
# dtype=complex<double>
D = np.column_stack((Dva, Dvb)) * w_f.reshape((20001, 1)) * np.ones((1, 13))
R = np.dot(D.conj().T, D)

Here is my C++ code (Minimal test):

#include <Eigen/Core>
#include <Eigen/Dense>
#include <vector>
#include <complex>

using namespace std;

typedef complex<double> dcomplex;

void foo()
{
    vector<dcomplex> wf;
    wf.resize(20001);
    Eigen::Matrix<dcomplex, 20001, 13> *tmp_1 = new Eigen::Matrix<dcomplex, 20001, 13>;
    Eigen::Matrix<dcomplex, 20001, 13> *tmp_2 = new Eigen::Matrix<dcomplex, 20001, 13>;
    Eigen::Matrix<dcomplex, 20001, 7> *Dva = new Eigen::Matrix<dcomplex, 20001, 7>;
    Eigen::Matrix<dcomplex, 20001, 6> *Dvb = new Eigen::Matrix<dcomplex, 20001, 6>;

    for (int i = 0; i < 20001; i++){
        for (int j = 0; j < 7; j++)
            (*Dva)(i, j) = 0;
        for (int j = 0; j < 6; j++)
            (*Dvb)(i, j) = 0;
        for (int j = 0; j < 13; j++)
            (*tmp_2)(i, j) = wf[i];
    }
    *tmp_1 << *Dva, *Dvb;
    auto *D = &tmp_1->cwiseProduct(*tmp_2);

    auto R = (D->transpose() * (*D));
    R(0,0);
}

The shape of matrix R is 13x13 in Eigen, which is same as NumPy. But the variable R can't be indicated in C++.

R.rows() == 13;  // true
R.cols() == 13;  // true
R(0, 0);  // or what ever makes it crash

Raised the exception "0xC00000FD: Stack overflow".

First of all, you should barely work with new in C++ code. Use local objects (or std::vector) most of the time, and if necessary use smart pointers like std::unique_ptr or std::shared_ptr.

Regarding the Eigen part of your question, avoid fixed sized matrices which are very large (more than a few KiB). You can have one dimension fixed and the other Dynamic. And finally, avoid auto in combination with Eigen, unless you know what you are doing!

The following should work. I replaced all your loops by corresponding Eigen functionality, and avoided the temporaries by directly doing a product with a diagonal matrix. Alternatively, you could do a cwiseProduct with a replicate()d matrix.

typedef Eigen::Matrix<dcomplex, Eigen::Dynamic,1> VectorXcd;
typedef Eigen::Matrix<dcomplex, Eigen::Dynamic,13> MatrixX13cd;
typedef Eigen::Matrix<dcomplex, Eigen::Dynamic,7> MatrixX7cd;
typedef Eigen::Matrix<dcomplex, Eigen::Dynamic,6> MatrixX6cd;
typedef Eigen::Matrix<dcomplex, 13,13> Matrix13cd;


MatrixX7cd  Dva(20001,  7);
MatrixX6cd  Dvb(20001,  6);

Dva.setZero(); Dvb.setZero();
MatrixX13cd D(20001, 13);
D.leftCols(7).noalias()  = VectorXcd::Map(wf.data(), wf.size()).asDiagonal() * Dva;
D.rightCols(6).noalias() = VectorXcd::Map(wf.data(), wf.size()).asDiagonal() * Dvb;

Matrix13cd R = D.transpose() * D;