From b6cef8dd016f8ead539ebd4fe5fc3161fc09e23c Mon Sep 17 00:00:00 2001
From: Oliver Hahn <oliver.hahn@univie.ac.at>
Date: Tue, 27 Feb 2024 01:59:49 +0800
Subject: [PATCH] updated N-GenIC RNG to multithreading

---
 src/plugins/random_music_wnoise_generator.cc | 111 ++++++++++---------
 1 file changed, 59 insertions(+), 52 deletions(-)

diff --git a/src/plugins/random_music_wnoise_generator.cc b/src/plugins/random_music_wnoise_generator.cc
index 324450f..ae75110 100644
--- a/src/plugins/random_music_wnoise_generator.cc
+++ b/src/plugins/random_music_wnoise_generator.cc
@@ -48,77 +48,84 @@ void music_wnoise_generator<T>:: gen_topgrid_NGenIC(size_t res, long baseseed) {
 		for (j = 0; j < i + 1; j++)
 			seedtable[(res - 1 - j) * res + (res - 1 - i)] = 0x7fffffff * gsl_rng_uniform(random_generator);
 	}
+  gsl_rng_free( random_generator );
 
 	real_t *rnoise = new real_t[res * res * (res + 2)];
 	complex_t *knoise = reinterpret_cast<complex_t *>(rnoise);
 
   double fnorm = pow((double)res, -1.5);
 
-// /#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;
+  // launch threads with indendent RNGs
+  #pragma omp parallel
+  {
+    gsl_rng *thread_rng = gsl_rng_alloc(gsl_rng_ranlxd1);
 
-    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);
+    #pragma omp for
+    for (size_t i = 0; i < res; i++) {
+      int ii = (int)res - (int)i;
+      if (ii == (int)res)
+        ii = 0;
 
-        if (i == res / 2 || j == res / 2 || k == res / 2)
-          continue;
-        if (i == 0 && j == 0 && k == 0)
-          continue;
+      for (size_t j = 0; j < res; j++) {
+        gsl_rng_set(thread_rng, seedtable[i * res + j]);
+        for (size_t k = 0; k < res / 2; k++) {
+          double phase = gsl_rng_uniform(thread_rng) * 2 * M_PI;
+          double ampl;
+          do {
+            ampl = gsl_rng_uniform(thread_rng);
+          } while (ampl == 0);
 
-        T rp = -sqrt(-log(ampl)) * cos(phase) * fnorm;
-        T ip = -sqrt(-log(ampl)) * sin(phase) * fnorm;
+          if (i == res / 2 || j == res / 2 || k == res / 2)
+            continue;
+          if (i == 0 && j == 0 && k == 0)
+            continue;
 
-        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;
+          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 {
-              int jj =
-                  (int)res - (int)j; /* note: j!=0 surely holds at this point */
+              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;
+                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;
+                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;
+                }
               }
             }
           }
         }
       }
     }
+    gsl_rng_free( thread_rng );
   }
 
   delete[] seedtable;