Aquila-Consortium/libsharp2
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

parallelize random a_lm generation

Browse source
This commit is contained in:
Martin Reinecke 2013-01-16 12:57:33 +01:00
parent 54bc25edcb
commit 842b790b1f
1 changed files with 14 additions and 6 deletions
Download patch file Download diff file Expand all files Collapse all files

20
libsharp/sharp_testsuite.c
View file

@ -46,29 +46,37 @@ typedef complex double dcmplx;
int ntasks, mytask; int ntasks, mytask;
static double drand (double min, double max) static double drand (double min, double max, int *state)
{ return min + (max-min)*rand()/(RAND_MAX+1.0); } {
*state = (((*state) * 1103515245) + 12345) & 0x7fffffff;
return min + (max-min)*(*state)/(0x7fffffff+1.0);
}
static void random_alm (dcmplx *alm, sharp_alm_info *helper, int spin) static void random_alm (dcmplx *alm, sharp_alm_info *helper, int spin)
{ {
static int cnt=0; static int cnt=0;
++cnt; ++cnt;
for (int mi=0;mi<helper->nm; ++mi) #pragma omp parallel
{
int mi;
#pragma omp for schedule (dynamic,100)
for (mi=0;mi<helper->nm; ++mi)
{ {
int m=helper->mval[mi]; int m=helper->mval[mi];
srand(1234567*cnt+8912*m); int state=1234567*cnt+8912*m; // random seed
for (int l=m;l<=helper->lmax; ++l) for (int l=m;l<=helper->lmax; ++l)
{ {
if ((l<spin)&&(m<spin)) if ((l<spin)&&(m<spin))
alm[sharp_alm_index(helper,l,mi)] = 0.; alm[sharp_alm_index(helper,l,mi)] = 0.;
else else
{ {
double rv = drand(-1,1); double rv = drand(-1,1,&state);
double iv = (m==0) ? 0 : drand(-1,1); double iv = (m==0) ? 0 : drand(-1,1,&state);
alm[sharp_alm_index(helper,l,mi)] = rv+_Complex_I*iv; alm[sharp_alm_index(helper,l,mi)] = rv+_Complex_I*iv;
} }
} }
} }
} // end of parallel region
} }
static unsigned long long totalops (unsigned long long val) static unsigned long long totalops (unsigned long long val)