mirror/MUSIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/MUSIC.git synced 2024-09-19 17:03:46 +02:00
Code Releases Activity

updated N-GenIC RNG to multithreading

Browse source
This commit is contained in:
Oliver Hahn 2024-02-27 01:59:49 +08:00
parent 75587bbf26
commit b6cef8dd01
1 changed files with 59 additions and 52 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

111
src/plugins/random_music_wnoise_generator.cc
View file

@ -48,77 +48,84 @@ void music_wnoise_generator<T>:: gen_topgrid_NGenIC(size_t res, long baseseed) {
for (j = 0; j < i + 1; j++) for (j = 0; j < i + 1; j++)
seedtable[(res - 1 - j) * res + (res - 1 - i)] = 0x7fffffff * gsl_rng_uniform(random_generator); 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)]; real_t *rnoise = new real_t[res * res * (res + 2)];
complex_t *knoise = reinterpret_cast<complex_t *>(rnoise); complex_t *knoise = reinterpret_cast<complex_t *>(rnoise);
double fnorm = pow((double)res, -1.5); double fnorm = pow((double)res, -1.5);
// /#warning need to check for race conditions below // launch threads with indendent RNGs
//#pragma omp parallel for #pragma omp parallel
for (size_t i = 0; i < res; i++) { {
int ii = (int)res - (int)i; gsl_rng *thread_rng = gsl_rng_alloc(gsl_rng_ranlxd1);
if (ii == (int)res)
ii = 0;
for (size_t j = 0; j < res; j++) { #pragma omp for
gsl_rng_set(random_generator, seedtable[i * res + j]); for (size_t i = 0; i < res; i++) {
for (size_t k = 0; k < res / 2; k++) { int ii = (int)res - (int)i;
double phase = gsl_rng_uniform(random_generator) * 2 * M_PI; if (ii == (int)res)
double ampl; ii = 0;
do {
ampl = gsl_rng_uniform(random_generator);
} while (ampl == 0);
if (i == res / 2 || j == res / 2 || k == res / 2) for (size_t j = 0; j < res; j++) {
continue; gsl_rng_set(thread_rng, seedtable[i * res + j]);
if (i == 0 && j == 0 && k == 0) for (size_t k = 0; k < res / 2; k++) {
continue; 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; if (i == res / 2 || j == res / 2 || k == res / 2)
T ip = -sqrt(-log(ampl)) * sin(phase) * fnorm; continue;
if (i == 0 && j == 0 && k == 0)
continue;
if (k > 0) { T rp = -sqrt(-log(ampl)) * cos(phase) * fnorm;
RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp; T ip = -sqrt(-log(ampl)) * sin(phase) * fnorm;
IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip;
} else /* k=0 plane needs special treatment */ if (k > 0) {
{ RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp;
if (i == 0) { IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip;
if (j >= res / 2) { } else /* k=0 plane needs special treatment */
continue; {
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 { } else {
int jj = if (i >= res / 2) {
(int)res - (int)j; /* note: j!=0 surely holds at this point */ 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; RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp;
IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip; IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip;
RE(knoise[(i * res + jj) * (res / 2 + 1) + k]) = rp; if (ii >= 0 && ii < (int)res) {
IM(knoise[(i * res + jj) * (res / 2 + 1) + k]) = -ip; RE(knoise[(ii * res + jj) * (res / 2 + 1) + k]) = rp;
} IM(knoise[(ii * 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;
} }
} }
} }
} }
} }
} }
gsl_rng_free( thread_rng );
} }
delete[] seedtable; delete[] seedtable;