vendredi 17 juin 2022

How to convert for loop matrix multiplication using Eigen C++: incorrect results

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