Eigen错误:用Eigen替换fftw_malloc数组后，c++静态断言失败

Question

我正试图切换我的代码，并开始使用c++中的Eigen库，因为我听说它对矩阵非常好。现在，我的旧c++代码主要使用fftw_malloc来初始化数组，如下所示:

static const int nx = 10;
static const int ny = 10; 
double *XX;
XX = (double*) fftw_malloc(nx*ny*sizeof(double));
memset(XX, 42, nx*ny* sizeof(double)); 
for(int i = 0; i< nx; i++){
for(int j = 0; j< ny+1; j++){
XX[j + ny*i] = (i)*2*M_PI/nx;
}       
}

所以我把它改成了

Matrix <double, nx, ny> eXX;
eXX.setZero();
for(int i = 0; i< nx; i++){
for(int j = 0; j< ny+1; j++){
eXX[j + ny*i] = (i)*2*EIGEN_PI/nx;
}       
}

返回错误:

In file included from /mnt/c/Users/J/Documents/eigen-3.4.0/eigen-3.4.0/Eigen/Core:164,
from /mnt/c/Users/J/Documents/eigen-3.4.0/eigen-3.4.0/Eigen/Dense:1,
from Test.cpp:21:
/mnt/c/Users/J/Documents/eigen-3.4.0/eigen-3.4.0/Eigen/src/Core/DenseCoeffsBase.h: In instantiation of ‘Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::operator[](Eigen::Index) [with Derived = Eigen::Matrix<double, 10, 10>; Eigen::DenseCoeffsBase<Derived, 1>::Scalar = double; Eigen::Index = long int]’:
Test.cpp:134:16:   required from here
/mnt/c/Users/J/Documents/eigen-3.4.0/eigen-3.4.0/Eigen/src/Core/DenseCoeffsBase.h:408:36: error: static assertion failed: THE_BRACKET_OPERATOR_IS_ONLY_FOR_VECTORS__USE_THE_PARENTHESIS_OPERATOR_INSTEAD
408 |       EIGEN_STATIC_ASSERT(Derived::IsVectorAtCompileTime,
|                                    ^~~~~~~~~~~~~~~~~~~~~
/mnt/c/Users/J/Documents/eigen-3.4.0/eigen-3.4.0/Eigen/src/Core/util/StaticAssert.h:33:54: note: in definition of macro ‘EIGEN_STATIC_ASSERT’
33 |     #define EIGEN_STATIC_ASSERT(X,MSG) static_assert(X,#MSG);
|                                                      ^

我看到错误指向行:eXX[j + ny*i] = (i)*2*EIGEN_PI/nx;我对Eigen相当陌生，仍然在努力理解我的方式，所以任何反馈/建议都会有所帮助。谢谢。

Answer 1

断言消息THE_BRACKET_OPERATOR_IS_ONLY_FOR_VECTORS__USE_THE_PARENTHESIS_OPERATOR_INSTEAD已经告诉您错误是什么:不要使用operator[]，而是使用operator()，如下所示:

#include <Eigen/Core>
int main()
{
static const int nx = 10;
static const int ny = 100;
Eigen::Matrix<double, nx, ny> eXX;
eXX.setZero();
for (int i = 0; i < nx; i++) {
for (int j = 0; j < ny; j++) {
eXX(i, j) = i * 2 * EIGEN_PI / nx;
}
}
}

还要注意，您的代码具有未定义的行为，因为由于循环条件中的+1,j循环转到<= ny，这意味着最后一次迭代导致对数组的越界访问。实际上，Eigen发出了一个断言(在调试中)。

此外，由于缺省情况下特征矩阵是列为主的，因此交换行和列循环更有效(但当然这只对较大的矩阵有效):

#include <Eigen/Core>
int main()
{
static const int nx = 10;
static const int ny = 100;
Eigen::Matrix<double, nx, ny> eXX;
eXX.setZero();
for (int j = 0; j < ny; j++) {
for (int i = 0; i < nx; i++) {
eXX(i, j) = i * 2 * EIGEN_PI / nx;
}
}
}

由于每行的元素在这里具有相同的值，我们也可以将row()与setConstant()一起使用以简化代码:

#include <Eigen/Core>
int main()
{
static const int nx = 10;
static const int ny = 100;
Eigen::Matrix<double, nx, ny> eXX;
for (int i = 0; i < nx; i++) {
eXX.row(i).setConstant(i * 2 * EIGEN_PI / nx);
}
}