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.
Aucun commentaire:
Enregistrer un commentaire