I am trying to switch my code and start using Eigen library in C since I have heard it is really good with matrices. Now my old C code used mostly fftw_malloc to initialize my arrays like the following:
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;
}
}
So I changed it to the following,
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;
}
}
This however returns the error:
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);
| ^
I see the error is pointing at line: eXX[j ny*i] = (i)*2*EIGEN_PI/nx;
I am fairly new to Eigen and still trying to understand my way around so any feedback/suggestions would help. Thanks.
CodePudding user response:
The assertion message THE_BRACKET_OPERATOR_IS_ONLY_FOR_VECTORS__USE_THE_PARENTHESIS_OPERATOR_INSTEAD
already tells you what the error is: Do not use operator[]
but instead use operator()
, like so:
#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;
}
}
}
Also note that your code has undefined behavior because the j
loop goes to <= ny
due to the 1
in the loop condition, meaning that the last iterations cause out-of-bound accesses to the array. Indeed, Eigen issues an assertion (in debug).
Moreover, since Eigen matrices are column major by default, it is more efficient to exchange the row and column loops (but of course this only matters for larger matrices):
#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;
}
}
}
Since the elements of each row have the same value here, we can also use row()
together with setConstant()
to simplify the code:
#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);
}
}