MATLAB vs C++ Outputs not matching

[JamieAl]

I have the following function of sum of sines and cosines in both C++ and MATLAB and I am using this function for testing purposes. However, for some reason the output text file of both codes isn't matching. I double checked my parenthesis and floating point precision but with no difference.

C++ code:

#define _USE_MATH_DEFINES // for C++
#include <cmath>
#include<math.h>
#include<stdio.h>
#include "spectralFunctions.h"
#include "arithmeticFunctions.h"
#include "fftw3.h"

static const int nx = 1024; 
static const int ny = 256;

int main(){			
	double Lx = 1.5E6; 
	double Ly = 1.E6;
        double n0 = 1.0e11;
// SET Coeff. of test perturbation function:
  
	 double ax = 0.1;
	 double ay = 0.2;
	 double axy = .6; //sum ax, ay, axy <1
	 double nx1 = 1.; 
	 double ny1 = 100.; 
	 double nxy1 = 3.; 
	 double nxy2 = 10.; 

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

	

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

	// Make the physical grid
	double kmax = makeSpatialMesh2D(dx, dy, XX, YY);
double *ne1;
	ne1 = (double*) fftw_malloc(nx*ny*sizeof(double));
for (int i = 0; i < nx; i++){
		for (int j = 0; j < ny; j++){
		  //ne1[j + ny*i] =  ((rand() % 101)/100. * (Ampl_2 - Ampl_1) + Ampl_1)*n0;
          ne1[j+ny*i] =  (ax * sin((2. *M_PI * nx1 * XX[j + ny*i])/Lx ) + ay * cos((2. *M_PI * ny1 * YY[j + ny*i])/Ly) + axy * cos((2. *M_PI *nxy1 * XX[j + ny*i])/Lx) * sin((2. *M_PI * nxy2* YY[j + ny*i])/Ly) ) *n0;
		 }
}
// free...	   

While the MATLAB code is:

%same initialization and setting constants
  prim_var_pert(:,:,2) = (ax*sin(2 * pi * XX * nx1 / L(1)) + ay * cos(2 * pi * YY * ny1 / L(2)) + ...
            axy*cos(2*pi*XX * nxy1/L(1)) .* sin(2*pi*YY * nxy2/L(2))) *n0;

Where L(1) is Lx and L(2) is Ly here, prim_var_pert(:,:,2) is just ne1 and so on.

The output is a 2D array written in a text file but somehow not all the lines match.

[JamieAl]

Just for you information, I am using this function to write a 2D array to a text file:

void print2DArrf(char filename[], double arr[]){

	FILE *fptr;
	
	fptr = fopen(filename, "w");
	
	
	for (int i = 0; i < nx; i++) {
		for (int j = 0 ; j < ny ; j++) {
			fprintf(fptr, "%+4.2le  ",arr[j + ny*i]);
		}
		fprintf(fptr, "\n");		
	}
	
	fclose(fptr);
}	 

[salem c]

> The output is a 2D array written in a text file but somehow not all the lines match.
Depends on what you mean by "not matching".

If you're talking about random +/-1 in the least significant digit of the results, then there's nothing to worry about, and probably little you can do about it.

Floating point numbers are approximations, so one approximation is as good as any other if the difference is minor.
This is required reading if you're doing any serious work with FP.
https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

[jonnin]

what do you think that memset does?
memset is a byte dump function. it says set every BYTE to a value. Not every double. A double is 8 bytes, and you set all 8 of them to an identical value. If that value is not zero, and you are doing doing something extremely weird, its likely those statements are wrong.

[JamieAl]

@jonnin
Memset has always been a confusing function to me. Here I am initializing my XX and YY to zeros with a specific size so that I can later use them in functions. This has been the only "initializing" way that worked for me with fftw_malloc

@salem c
The type of difference b/w results I am talking about is like this:

MATLAB output:

2.0000000e+10      -8.8139827e+08    3.2185611e+10  ...
2.0061359e+10      -8.2250934e+08    3.2242178e+10  ....
2.0122715e+10      -7.6856172e+08    3.2289160e+10  ....
.
.
.

C++ output:

+2.00e+10      -8.81e+08   +3.22e+10  ...
5.39e+10        -5.68e+10    -4.27e+10  ....
.86e+10         +7.56e+10    +3.14e+10  ....
.
.
.

As you can see, things start okay but suddenly change...

Obviously it's hard to compare by just looking or simply posting this here, since these are large text files. So, instead I found I loaded these files on MATLAB and found the difference between the two, and there's a significant difference.

[jonnin]

you can clear doubles to zero with memset, but just use 0, not 42. You set each individual byte of the doubles to 42. If that somehow makes zero, I do not understand it.

as far as the two result streams drifting apart, the easy way is to check each term on both sets of code for the first 10 or something terms. see which term(s) are not correct, and we can look at that closer.

[helios]

8 consecutive bytes holding 42 is not the double 0.0. According to Ganado it's 1.42603e-105. OP was already told not to use 42 with memset() back in April.

[JamieAl]

I see now...

I just recall that using 0.0 with memset gave same answer as using 42 with memset. That's why I totally forgot about this.

I am testing each term separately as @jonnin recommended.

[jonnin]

hmm... its pretty darn close to zero, though. Interesting trick, if mostly useless. I suddenly feel compelled to learn all 255 of the possible doubles to see if any ARE useful...

you may find it to be an order of operations issue, as well. its aggravating but if you keep comparing the inputs and terms you will find it going wrong somewhere. That is more than roundoff error in your output.

[dhayden]

If you're still having trouble, check the results of makeSpatialMesh2D() too.

[JamieAl]

Hey everyone!

It seems the issue was a huge "round off" error. I didn't think to match the number of significant figures of MATLAB's output to the C++. When I set both to 16 sig figs that seemed to solve the issue. I clearly missed that somehow although it was almost obvious. Thanks!