Методы Python Numpy соответствуют C ++ Eigen make crash

Question

У меня есть арифметика NumPy для перевода C ++ с 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)

Вот мой код C ++ (минимальный тест):

#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);
}

Форма матрицы R равна 13x13 в Eigen, что соответствует NumPy. Но переменная R не может быть указана в C ++.

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

Возникла исключительная ситуация "0xC00000FD: переполнение стека".

Answer 1

Прежде всего, вам едва нужно работать с new в коде C ++.Используйте локальные объекты (или std::vector) большую часть времени, и, если необходимо, используйте умные указатели, такие как std::unique_ptr или std::shared_ptr.

Что касается собственной части вашего вопроса, избегайте матриц фиксированного размера, которые оченьбольшой (более нескольких КиБ).У вас может быть одно фиксированное измерение, а другое Dynamic.И, наконец, избегайте auto в сочетании с Эйгеном, если только вы не знаете, что делаете!

Следующее должно работать.Я заменил все ваши циклы на соответствующие функции Eigen и избежал временных изменений, выполнив непосредственно продукт с диагональной матрицей.Кроме того, вы можете сделать cwiseProduct с матрицей replicate() d .

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;