WIP fixes/cleanup for MUSIC RNG

src/plugins/random_music_wnoise_generator.cc

@@ -7,6 +7,157 @@
 #include "random.hh"
 #include "random_music_wnoise_generator.hh"
 
+
+template <typename T>
+void rapid_proto_ngenic_rng(size_t res, long baseseed, music_wnoise_generator<T> &R)
+{
+	LOGUSER("Invoking the N-GenIC random number generator");
+
+	unsigned *seedtable = new unsigned[res * res];
+
+	gsl_rng *random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);
+
+	gsl_rng_set(random_generator, baseseed);
+
+	for (size_t i = 0; i < res / 2; i++)
+	{
+		size_t j;
+		for (j = 0; j < i; j++)
+			seedtable[i * res + j] = 0x7fffffff * gsl_rng_uniform(random_generator);
+		for (j = 0; j < i + 1; j++)
+			seedtable[j * res + i] = 0x7fffffff * gsl_rng_uniform(random_generator);
+		for (j = 0; j < i; j++)
+			seedtable[(res - 1 - i) * res + j] = 0x7fffffff * gsl_rng_uniform(random_generator);
+		for (j = 0; j < i + 1; j++)
+			seedtable[(res - 1 - j) * res + i] = 0x7fffffff * gsl_rng_uniform(random_generator);
+		for (j = 0; j < i; j++)
+			seedtable[i * res + (res - 1 - j)] = 0x7fffffff * gsl_rng_uniform(random_generator);
+		for (j = 0; j < i + 1; j++)
+			seedtable[j * res + (res - 1 - i)] = 0x7fffffff * gsl_rng_uniform(random_generator);
+		for (j = 0; j < i; j++)
+			seedtable[(res - 1 - i) * res + (res - 1 - j)] = 0x7fffffff * gsl_rng_uniform(random_generator);
+		for (j = 0; j < i + 1; j++)
+			seedtable[(res - 1 - j) * res + (res - 1 - i)] = 0x7fffffff * gsl_rng_uniform(random_generator);
+	}
+
+	fftw_real *rnoise = new fftw_real[res * res * (res + 2)];
+	fftw_complex *knoise = reinterpret_cast<fftw_complex *>(rnoise);
+
+	double fnorm = 1. / sqrt(res * res * res);
+
+// #warning need to check for race conditions below
+	//#pragma omp parallel for
+	for (size_t i = 0; i < res; i++)
+	{
+		int ii = (int)res - (int)i;
+		if (ii == (int)res)
+			ii = 0;
+
+		for (size_t j = 0; j < res; j++)
+		{
+			gsl_rng_set(random_generator, seedtable[i * res + j]);
+
+			for (size_t k = 0; k < res / 2; k++)
+			{
+				double phase = gsl_rng_uniform(random_generator) * 2 * M_PI;
+				double ampl;
+				do
+					ampl = gsl_rng_uniform(random_generator);
+				while (ampl == 0);
+
+				if (i == res / 2 || j == res / 2 || k == res / 2)
+					continue;
+				if (i == 0 && j == 0 && k == 0)
+					continue;
+
+				T rp = -sqrt(-log(ampl)) * cos(phase) * fnorm;
+				T ip = -sqrt(-log(ampl)) * sin(phase) * fnorm;
+
+				if (k > 0)
+				{
+					RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp;
+					IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip;
+				}
+				else /* k=0 plane needs special treatment */
+				{
+					if (i == 0)
+					{
+						if (j >= res / 2)
+							continue;
+						else
+						{
+							int jj = (int)res - (int)j; /* note: j!=0 surely holds at this point */
+
+							RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp;
+							IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip;
+
+							RE(knoise[(i * res + jj) * (res / 2 + 1) + k]) = rp;
+							IM(knoise[(i * res + jj) * (res / 2 + 1) + k]) = -ip;
+						}
+					}
+					else
+					{
+						if (i >= res / 2)
+							continue;
+						else
+						{
+							int ii = (int)res - (int)i;
+							if (ii == (int)res)
+								ii = 0;
+							int jj = (int)res - (int)j;
+							if (jj == (int)res)
+								jj = 0;
+
+							RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp;
+							IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip;
+
+							if (ii >= 0 && ii < (int)res)
+							{
+								RE(knoise[(ii * res + jj) * (res / 2 + 1) + k]) = rp;
+								IM(knoise[(ii * res + jj) * (res / 2 + 1) + k]) = -ip;
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+
+	delete[] seedtable;
+
+	//... perform FT to real space
+
+#ifdef FFTW3
+#ifdef SINGLE_PRECISION
+	fftwf_plan plan = fftwf_plan_dft_c2r_3d(res, res, res, knoise, rnoise, FFTW_ESTIMATE);
+	fftwf_execute(plan);
+	fftwf_destroy_plan(plan);
+#else
+	fftw_plan plan = fftw_plan_dft_c2r_3d(res, res, res, knoise, rnoise, FFTW_ESTIMATE);
+	fftw_execute(plan);
+	fftw_destroy_plan(plan);
+#endif
+#else
+	rfftwnd_plan plan = rfftw3d_create_plan(res, res, res, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_IN_PLACE);
+#ifndef SINGLETHREAD_FFTW
+	rfftwnd_threads_one_complex_to_real(omp_get_max_threads(), plan, knoise, NULL);
+#else
+	rfftwnd_one_complex_to_real(plan, knoise, NULL);
+#endif
+	rfftwnd_destroy_plan(plan);
+#endif
+
+	// copy to array that holds the random numbers
+
+#pragma omp parallel for
+	for (int i = 0; i < (int)res; ++i)
+		for (size_t j = 0; j < res; ++j)
+			for (size_t k = 0; k < res; ++k)
+				R(i, j, k) = rnoise[((size_t)i * res + j) * res + k];
+
+	delete[] rnoise;
+}
+
 template <typename T>
 music_wnoise_generator<T>::music_wnoise_generator(unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx)
     : res_(res), cubesize_(cubesize), ncubes_(1), baseseed_(baseseed)
@@ -35,18 +186,6 @@ music_wnoise_generator<T>::music_wnoise_generator(unsigned res, unsigned cubesiz
 
   initialize();
 
-  if (musicnoise)
-    mean = fill_all();
-  else
-  {
-    rnums_.push_back(new Meshvar<T>(res, 0, 0, 0));
-    cubemap_[0] = 0; // create dummy map index
-    register_cube(0, 0, 0);
-    // rapid_proto_ngenic_rng( res_, baseseed_, *this );
-  }
-
-  /*
-
   if (musicnoise)
     mean = fill_all();
   else
@@ -57,8 +196,6 @@ music_wnoise_generator<T>::music_wnoise_generator(unsigned res, unsigned cubesiz
     rapid_proto_ngenic_rng( res_, baseseed_, *this );
   }
 
-  */
-
   if (zeromean)
   {
     mean = 0.0;
@@ -847,8 +984,7 @@ double music_wnoise_generator<T>::fill_all(void)
         register_cube(ii, jj, kk);
       }
 
-#pragma omp parallel for reduction(+ \
+  #pragma omp parallel for reduction(+ : sum)
-                                   : sum)
   for (int i = 0; i < (int)ncubes_; ++i)
     for (int j = 0; j < (int)ncubes_; ++j)
       for (int k = 0; k < (int)ncubes_; ++k)
@@ -863,8 +999,7 @@ double music_wnoise_generator<T>::fill_all(void)
       }
 
   //... subtract mean
-#pragma omp parallel for reduction(+ \
+  #pragma omp parallel for reduction(+ : sum)
-                                   : sum)
   for (int i = 0; i < (int)ncubes_; ++i)
     for (int j = 0; j < (int)ncubes_; ++j)
       for (int k = 0; k < (int)ncubes_; ++k)

src/plugins/random_music_wnoise_generator.hh

@@ -1,5 +1,4 @@
-#ifndef __RANDOM_MUSIC_WNOISE_GENERATOR_HH
+#pragma once
-#define __RANDOM_MUSIC_WNOISE_GENERATOR_HH
 
 #include <vector>
 #include <map>
@@ -31,7 +30,6 @@ protected:
   typedef std::map<size_t, size_t>::iterator cubemap_iterator;
 
 protected:
-  
   //! register a cube with the hash map
   void register_cube(int i, int j, int k);
 
@@ -98,7 +96,8 @@ protected:
           register_cube(ii, jj, kk);
         }
 
-    #pragma omp parallel for reduction(+:sum)
+#pragma omp parallel for reduction(+ \
+                                   : sum)
     for (int i = 0; i < (int)ncubes_; ++i)
       for (int j = 0; j < (int)ncubes_; ++j)
         for (int k = 0; k < (int)ncubes_; ++k)
@@ -121,7 +120,6 @@ protected:
   void print_allocated(void);
 
 public:
-  
   //! constructor
   music_wnoise_generator(unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx);
 
@@ -131,11 +129,9 @@ public:
   //! constructor
   music_wnoise_generator(unsigned res, unsigned cubesize, long baseseed, bool zeromean = true);
 
-  
   //! constructor to read white noise from file
   music_wnoise_generator(unsigned res, std::string randfname, bool rndsign);
 
-  
   //! copy constructor for averaged field (not copying) hence explicit!
   explicit music_wnoise_generator(/*const*/ music_wnoise_generator<T> &rc);
 
@@ -199,7 +195,3 @@ public:
   }
 };
 
-
-
-
-#endif //__MUSIC_WNOISE_GENERATOR_HH