Error: MATLAB vs C++ output

[JamieAl]

I am converting a MATLAB code to a C++ and running into small differences b/w my two codes.

For some reason in the C++ code, each term in dVGL is correct when printed out "separately" but somehow when put togther dVGL is incorrect. I get "Inf" and "nan" for some of the values.

I wrote a loop that basically calculates the diagonal matrix given a vector (similar to diag function in MATLAB) and initialized dummy variables that can be used in the final calculation of dVGL.

C++ code:

  
    std::vector<double> ygl(ny+1);  

    for (int i = 0; i< ny+1; i++){
		ygl[i] = -1. * cos(((i) * M_PI )/ny); //This works
		//cout << ygl[i] << endl ; //test
	}

//initialize variables
        double *VGL; 
	VGL = (double*) fftw_malloc(nx*ny*sizeof(double));
	memset(VGL, 42, nx*ny* sizeof(double));

	double *dVGL; 
	dVGL = (double*) fftw_malloc(nx*ny*sizeof(double));
	memset(dVGL, 42, nx*ny* sizeof(double));


	double *dummy1; 
	dummy1 = (double*) fftw_malloc(nx*ny*sizeof(double));
	memset(dummy1, 42, nx*ny* sizeof(double));

	double *dummy2; 
	dummy2 = (double*) fftw_malloc(nx*ny*sizeof(double));
	memset(dummy2, 42, nx*ny* sizeof(double));


	for(int i = 0; i< ny+1; i++){ //ny?
		for(int j = 0; j< ny+1; j++){
			if(i==j){
				dummy1[j + ny*i] =  1./(sqrt(1-pow(ygl[j],2))); 
				dummy2[j + ny*i] = 1. * (i);
			}else {
				dummy1[j + ny*i] = 0;
				dummy2[j + ny*i] = 0;
			}
			VGL[j + ny*i] = cos( acos(ygl[j]) *(i)); //VGL correct
			dVGL[j + ny*i] = dummy1[j + ny*i] * sin(acos(ygl[j]) * (i)) * dummy2[j + ny*i]; // dummy1, dummy2, and sin(acos(ygl[j]) * (i)) are correct separately, but not dVGL[j + ny*i]
	
			cout << dVGL[j + ny*i]<< endl ;
			
		}		
	}

This is my MATLAB code:

Ny = 10; 
ygl = -cos(pi*(0:Ny)/Ny)'; %Gauss-Lobatto chebyshev points 

VGL = cos(acos(ygl(:))*(0:Ny));
dVGL = diag(1./sqrt(1-ygl.^2))*sin(acos(ygl)*(0:Ny))*diag(0:Ny);  

Thanks

[coder777]

By using ny+1/ny+1 in the for loop (line 28/29 and probably 4) the indexes i/j are out of bounds at the end. This causes undefined behavior.

You can check this with an assert:
assert((j + ny*i) < (nx*ny));

[jonnin]

it may be easier to wire your c++ with extra space and use 1-N instead of 0..N-1 indexing so the code is closer to the same logic and lines.

once you get it debugged, do not be alarmed if the answers are still slightly different. Matlab has all sorts of numerical methods that it chooses from when doing things and those algorithms are not always documented for the user. If you get within some reasonable margin of error to the same answer, I would check it but keep in mind that it could just be close enough.

dummy2[j + ny*i] = 1. * (i); //this multiply just casts, use a cast instead.

[JamieAl]

Thanks all for the help.

@coder777 I used this new debugging technique and it helps a lot, now it returns the message:

test.out: Fun.cpp:216: int main(): Assertion `(j + ny*i) < (nx*ny)' failed.
Aborted

So you're right, one thing I was trying to do is have the size of my matrix to be Ny+1 by Ny+1 and wasn't sure how to use memset correctly to do that without returning some error:

malloc(): corrupted top size
Aborted

@jonnin
Not sure what you mean by extra space in the code.

use 1-N instead of 0..N-1 indexing

I am struggling to index properly here for some reason, it isn't as direct as MATLAB. Unfortunately, I am still getting results off with lots of zeros.

[JamieAl]

okay I guess I see the issue, the C++ code is giving me the diagonal matrix with zeros (off-diagonal terms) of the correct answer, which is interesting. So I though then maybe I need to create a separate loop outside to calculate dVGL like:

	for(int i = 0; i< ny+1; i++){ //ny?
		for(int j = 0; j< ny+1; j++){
			if(i==j){
				dummy1[j + ny*i] =  1./(sqrt(1-pow(ygl[j],2))); //dummy1[j*ny+i] =  1./(sqrt(1-pow(ygl[j],2)));
				dummy2[j + ny*i] = 1. * (i);
			}else {
				dummy1[j + ny*i] = 0;
				dummy2[j + ny*i] = 0;
			}
			VGL[j + ny*i] = cos( acos(ygl[j]) *(i)); //VGL correct
			//dVGL[j + ny*i] =  dummy1[j + ny*i] * sin(acos(ygl[j]) * (i)) * dummy2[j + ny*i];
			//dummy1[j + ny*i] * sin(acos(ygl[j]) * (i)) * dummy2[j + ny*i];
			//dummy1[j + ny*i] * sin(acos(ygl[j]) * (i)) * (dummy2[j + ny*i]);
			//	dVGL = diag(1./sqrt(1-ygl.^2))*sin(acos(ygl)*(0:Ny))*diag(0:Ny); 
			//assert((j + ny*i) < (nx*ny)); //this macro is designed to capture programming errors
			//cout << dVGL[j + ny*i]<< endl ;
			
		}	
			
	}
	for(int i = 0; i< ny+1; i++){
		for(int j = 0; j< ny+1; j++){
			dVGL[j + ny*i] =  dummy1[j + ny*i] * sin(acos(ygl[j]) * (i)) * dummy2[j + ny*i];
			cout << dVGL[j + ny*i]<< endl ;

		}

	}

However, this returns:

malloc_consolidate(): invalid chunk size
Aborted

which is odd??!

[jonnin]

the c++ formula for 2d from 1d is
[desired row * number of cols + desired col]

I don't know what random letter variables are, so its hard to point to errors.
some of this stuff may be much easier with a valarray.

for(int i = 0; i< ny+1; i++){ //ny?
that just looks wrong.
most c++ for loops iterate eg
for(int i = 0; i< ny; i++) //happens NY times.

[George P]

@JamieAl, I'm coming into this a bit late, but I noticed in your first code snippet you create a std::vector. Using a C++ container such as a vector removes a lot of the nitty-gritty details of manual memory management and makes it less of a pain.

With a vector you can get/release memory on demand. You can create a sized vector from the start and add/release memory as needed, or create an empty vector and resize as needed.

You've created a 1D vector, creating a true 2D vector (a vector of vectors) isn't that difficult.

Empty:
std::vector<std::vector<int>> a2Dvec;

Easier to deal with is one sized at creation:

int row_size { 3 };
int col_size { 4 };
std::vector<std::vector<int>> aSized2DVec(row_size, std::vector<int>(col_size));

With a sized 2D vector you can access the elements using operator[] as if it were a regular 2D array: aSize2DVec[1][2] accesses the 2nd vector (row) and the 3rd element of that vector (column).

[JamieAl]

@jonnin

Okay, I was able to fix the main "out of band issue" and get the correct results. Thanks everyone!