C++ vs MATLAB for loops syntax

[JamieAl]

I was trying to use the function rand() in C++ to reproduce some for loop results from MATLAB. Are these to for loops in MATLAB and C++ identical? Or I guess how to fix this to be similar to MATLAB syntax. I am not really able to compare easily since I am using rand here.

MATLAB syntax:

Nx = 8;  
Ny = 8;
//meshgrid: [X,Y] = meshgrid(x,y);
double Ly = 128;
pertParam(1) = 1.e-6; 
pertParam(2) = 1.e-5; 
pertParam(3) = 0;     
pertParam(4) = 2*pi;  
pertParam(5) = 64;    
pertParam(6) = 128;    
pertY = 0;
for i = 1:pertParam(6)
            pertY = pertY + ((pertParam(2) - pertParam(1)) * rand + pertParam(1))...
                * cos( 2 * pi / Ly * i * Y2 + ((pertParam(4) -  ...
                pertParam(3)) * rand + pertParam(3)));
end

This is what I have in C++:

double pertParam1 = 1.e-6; 
double pertParam2 = 1.e-5; 
double pertParam3 = 0;     
double pertParam4 = 2*M_PI;  
double pertParam5 = 64;     
double pertParam6 = 128;     

	double *pertY;
	pertY = (double*) fftw_malloc(nx*ny*sizeof(double));
	memset(pertY, 42, nx*ny* sizeof(double)); 
	for (int i = 1; i < nx; i++){
		for (int j = 0; j < ny; j++){
			for (int k = 1; k<pertParam6; k++){
				pertY[j + ny*i] += ((pertParam2 - pertParam1) * ((rand() % 101)/100.) + pertParam1) * cos( (2 * M_PI / Ly) * k * YY[j + ny*i] + ((pertParam4 - pertParam3) * ((rand() % 101)/100.) + pertParam3));


			}

		}
	}

Is the above syntax correct? Thanks.

[seeplus]

Are you using c or C++ as you mention C++ but then show c code...

Why are you initialising pertY to 42? I know it's the answer to everything but in this context...

[Ganado]

seeplus wrote:

as you mention C++ but then show c code...

This is some unnecessary obtuseness/"gatekeeping" (for lack of a better term).

JamieAI,
Neither your Matlab nor your C++ code is complete, so we can't easily run your code to verify anything.

I am not really able to compare easily since I am using rand here.

Factor out the call to rand to something you can control, then. e.g. a generator that just goes {1, 2, 3, 1, 2, 3, ...}

e.g.

int my_rand()
{
#ifdef MATLAB_TEST
    static int i = -1;
    i = (i + 1) % 101;
    return i;
#else 
    return rand();
#endif

}

Or, as a general tip for hypothetical future scenarios where you can't control the data, run multiple trials and statistically compare the outputs to see if they are significantly different from each other (p-test, null hypothesis accepted/rejected).

[jonnin]

for something this aggravating, you should replace the random numbers with a known stream or known value.
eg
double rands[] { 0,1,2,3,4,5...}; //whatever values, have matlab generate this maybe
...
int index{};

... in the loops:
rands[index++] //make sure you generated enough that index won't go out of bounds

and then compare the c++ result to the matlab result for a few different rand sequences. you can expect some tiny differences in double values, moreso if you call any matlab routines inside the code (none here that I saw).

if you get the expected output, its good. This is much safer than trying to eyeball it for some subtle difference or trying to catch the old matlab index 1 c++ index 0 or matlab variable type changed midstream type bugs.

Anything you are not absolutely sure of, you should do a side by side compare of a small isolated test run.

don't forget to put your random numbers back in, and consider using the c++ <random> library because rand() has a LOT of issues (both with randomness and lack of features).
also consider using a vector, which prevents having to use C malloc or C++ pointer nonsense. Dynamic memory is a recipe for problems, esp if you are not used to it, and vectors hide this.
simply
#include <vector>
and this
pertY = (double*) fftw_malloc(nx*ny*sizeof(double));
becomes this:
std::vector<double> pertY(nx*ny);
and you don't even change the code otherwise:
pertY[index];//same as a pointer when accessing!
it will self delete, and you can iterate it with a range based for loop:
for(auto &py:pertY)
cout << py; //py becomes a reference to each item in pertY, it stands for pertyY[i] in a traditional loop.
Important: without the & py is a copy, and you can't change pertY (changes to py are lost). With the &, the changes are kept, as its a reference to the original. To keep it simple, just always use the &, and use const if you don't plan to change it, because copies can be expensive waste of time.

[Ganado]

also consider using a vector, which prevents having to use C malloc or C++ pointer nonsense.

std::vector<double> pertY(nx*ny);

No, do not do this. At least not without a custom allocator(*). The point of using fftw_malloc is that it reserves FFTW the right to change the memory-allocation implementation to better suit its needs. (Edit: Actually, according to their website, it's so "they guarantee that the returned pointer obeys any special alignment restrictions imposed by any algorithm in FFTW".)

But if you wanted to make your own RAII wrapper that uses fftw_malloc/destroy, then that's fine. (*) I think one of the experienced members on this site once set up the correct vector template for working with FFTW, but I can't find immediately find the thread. (Or, this might not be possible with std::vector. I don't remember.)

[helios]

But if you wanted to make your own RAII wrapper that uses fftw_malloc/destroy, then that's fine. (*) I think one of the experienced members on this site once set up the correct vector template for working with FFTW, but I can't find immediately find the thread. (Or, this might not be possible with std::vector. I don't remember.)

Just use std::unique_ptr:

struct FftwRelease{
    void operator()(void *p){
        fftw_free(p);
    }
};

template <typename T>
std::unique_ptr<T, FftwRelease> fftw_new(size_t size){
    //Might be wise to throw on null.
    return {(T *)fftw_malloc(size * sizeof(T)))};
}

[JamieAl]

Thanks everyone! These were all VERY helpful comments/suggestions and helped me a lot to compare my code!

[mbozzi]

I think one of the experienced members on this site once set up the correct vector template for working with FFTW, but I can't immediately find the thread.

That was me, using std::pmr::vector here
https://cplusplus.com/forum/beginner/283931/#msg1229277

[jonnin]

ouch, twice I have looked at that and took it to be some C extension instead of a full library. Apologies again.