Home > Blockchain >  How to convert for loop matrix multiplication using Eigen C : incorrect results
How to convert for loop matrix multiplication using Eigen C : incorrect results

Time:06-21

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
  • Related