I have the following for loop that works:
static const int ny = 10;
std::vector<double> ygl(ny 1);
double *dVGL;
dVGL = (double*) fftw_malloc(((ny 1)*(ny 1))*sizeof(double));
memset(dVGL, 42, ((ny 1)*(ny 1))* sizeof(double));
double *dummy1;
dummy1 = (double*) fftw_malloc(((ny 1)*(ny 1))*sizeof(double));
memset(dummy1, 42, ((ny 1)*(ny 1))* sizeof(double));
double *dummy2;
dummy2 = (double*) fftw_malloc(((ny 1)*(ny 1))*sizeof(double));
memset(dummy2, 42, ((ny 1)*(ny 1))* sizeof(double));
for (int i = 0; i < ny 1; i ){
for (int j = 0; j < ny 1; j ){
ygl[j] = -1. * cos(((j) * EIGEN_PI )/ny);
dummy1[j ny*i] = 1./(sqrt(1-pow(ygl[j],2)));
dummy2[j ny*i] = 1. * (i);
dVGL[j ny*i] = dummy1[j ny*i] * sin(acos(ygl[j]) * (i)) * dummy2[j ny*i];
}
}
The above works fine, now I convert to Eigen and I am getting off results:
Eigen::Matrix< double, 1, ny 1> v1;
v1.setZero();
std::iota(v1.begin(), v1.end(), 0);
Eigen::Matrix< double, ny 1, ny 1> dummy1;
dummy1.setZero();
Eigen::Matrix< double, ny 1, ny 1> dummy2;
dummy2.setZero();
for (int j = 0; j < ny 1; j ){
v[j] = 1./(sqrt(1-pow(ygl[j],2)));
}
dummy1 = v.array().matrix().asDiagonal(); //this is correct
dummy2 = v1.array().matrix().asDiagonal(); //this is correct
Eigen::Matrix< double, ny 1, ny 1> dVGL;
dVGL.setZero();
for (int i = 0; i < ny 1; i ){
for (int j = 0; j < ny 1; j ){
ygl[j] = -1. * cos(((j) * EIGEN_PI )/ny);
dVGL(j ny*i) = sin(acos(ygl[j]) * (i)); //this is correct
}
}
//##################################
dv1 = (dummy1) * (dVGL) * (dummy2); //this is incorrect, dVGL outside for loop is different from inside for loop!!
I have been wracking my brain for a week now over this I cannot seem to fix it. why is dVGL
out side the for loop is different?? (off as if my rows and columns have length ny 1, but When I flatten my array I am using just ny.) why is that?
I feel like this should be a simple issue.
CodePudding user response:
Thanks to @chtz in the comments, I was able to fix this problem simply by changing my indexing:
dVGL(j (ny 1)*i) = sin(acos(ygl[j]) * (i));
dv1 = (dummy1) * (dVGL) * (dummy2); //Then this is correct