Python Numpy 方法对应于 C++ Eigen make crash

Python Numpy methods correspond to C++ Eigen make crash

本文关键字：C++ Eigen make crash Numpy 方法 Python      更新时间：2023-10-16

我有一个 NumPy 算术需要用特征翻译C++。 

# 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，与 NumPy 相同。但是变量 R 不能用 C++ 表示。

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

引发异常"0xC00000FD：堆栈溢出"。

首先，你应该在C++代码中几乎不使用new。大多数时候使用本地对象(或std::vector(，如有必要，请使用智能指针，如std::unique_ptr或std::shared_ptr。

关于问题的特征部分，请避免使用非常大的固定大小矩阵(超过几 KiB(。您可以固定一个维度，而另一个维度Dynamic 。最后，避免与本征结合使用auto，除非您知道自己在做什么！

以下应该有效。我用相应的 Eigen 功能替换了所有循环，并通过直接使用对角矩阵做产品来避免临时循环。或者，您可以使用 replicate() d 矩阵执行cwiseProduct。

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;
相关文章：