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

preparation for more flexible arrangement of data in the phase array

Browse source
This commit is contained in:
Martin Reinecke 2013-01-10 10:45:29 +01:00
parent 03d4880cc5
commit 9b1c631ce1
3 changed files with 15 additions and 9 deletions
Download patch file Download diff file Expand all files Collapse all files

15
libsharp/sharp.c
View file

@ -427,7 +427,12 @@ static void init_output (sharp_job *job)
}
static void alloc_phase (sharp_job *job, int nm, int ntheta)
{ job->phase=RALLOC(dcmplx,2*job->ntrans*job->nmaps*nm*ntheta); }
{
if ((nm&1023)==0) nm+=3; // hack to avoid critical strides
job->s_m=2*job->ntrans*job->nmaps;
job->s_th=job->s_m*nm;
job->phase=RALLOC(dcmplx,2*job->ntrans*job->nmaps*nm*ntheta);
}
static void dealloc_phase (sharp_job *job)
{ DEALLOC(job->phase); }
@ -436,7 +441,7 @@ static void dealloc_phase (sharp_job *job)
static void map2phase (sharp_job *job, int mmax, int llim, int ulim)
{
if (job->type != SHARP_MAP2ALM) return;
int pstride = 2*job->ntrans*job->nmaps;
int pstride = job->s_m;
#pragma omp parallel
{
ringhelper helper;
@ -444,7 +449,7 @@ static void map2phase (sharp_job *job, int mmax, int llim, int ulim)
#pragma omp for schedule(dynamic,1)
for (int ith=llim; ith<ulim; ++ith)
{
int dim2 = pstride*(ith-llim)*(mmax+1);
int dim2 = job->s_th*(ith-llim);
for (int i=0; i<job->ntrans*job->nmaps; ++i)
ringhelper_pair2phase(&helper,mmax,&job->ginfo->pair[ith], job->map[i],
&job->phase[dim2+2*i], &job->phase[dim2+2*i+1], pstride, job->flags);
@ -589,7 +594,7 @@ static void almtmp2alm (sharp_job *job, int lmax, int mi)
static void phase2map (sharp_job *job, int mmax, int llim, int ulim)
{
if (job->type == SHARP_MAP2ALM) return;
int pstride = 2*job->ntrans*job->nmaps;
int pstride = job->s_m;
#pragma omp parallel
{
ringhelper helper;
@ -597,7 +602,7 @@ static void phase2map (sharp_job *job, int mmax, int llim, int ulim)
#pragma omp for schedule(dynamic,1)
for (int ith=llim; ith<ulim; ++ith)
{
int dim2 = pstride*(ith-llim)*(mmax+1);
int dim2 = job->s_th*(ith-llim);
for (int i=0; i<job->ntrans*job->nmaps; ++i)
ringhelper_phase2pair(&helper,mmax,&job->phase[dim2+2*i],
&job->phase[dim2+2*i+1],pstride,&job->ginfo->pair[ith],job->map[i],

8
libsharp/sharp_core_inc2.c
View file

@ -653,7 +653,7 @@ static void Z(inner_loop) (sharp_job *job, const int *ispair,
{
for (int j=0; j<njobs; ++j)
{
int phas_idx = 2*(j+njobs*(itot*job->ainfo->nm+mi));
int phas_idx = itot*job->s_th + mi*job->s_m + 2*j;
complex double r1 = p1[j].s.r[i] + p1[j].s.i[i]*_Complex_I,
r2 = p2[j].s.r[i] + p2[j].s.i[i]*_Complex_I;
job->phase[phas_idx] = r1+r2;
@ -693,7 +693,7 @@ static void Z(inner_loop) (sharp_job *job, const int *ispair,
{
for (int j=0; j<njobs; ++j)
{
int phas_idx = 4*(j+njobs*(itot*job->ainfo->nm+mi));
int phas_idx = itot*job->s_th + mi*job->s_m + 4*j;
complex double q1 = p1[j].s.qr[i] + p1[j].s.qi[i]*_Complex_I,
q2 = p2[j].s.qr[i] + p2[j].s.qi[i]*_Complex_I,
u1 = p1[j].s.ur[i] + p1[j].s.ui[i]*_Complex_I,
@ -736,7 +736,7 @@ static void Z(inner_loop) (sharp_job *job, const int *ispair,
{
for (int j=0; j<njobs; ++j)
{
int phas_idx = 2*(j+njobs*(itot*job->ainfo->nm+mi));
int phas_idx = itot*job->s_th + mi*job->s_m + 2*j;
dcmplx ph1=job->phase[phas_idx];
dcmplx ph2=ispair[itot] ? job->phase[phas_idx+1] : 0.;
p1[j].s.r[i]=creal(ph1+ph2); p1[j].s.i[i]=cimag(ph1+ph2);
@ -766,7 +766,7 @@ static void Z(inner_loop) (sharp_job *job, const int *ispair,
{
for (int j=0; j<njobs; ++j)
{
int phas_idx = 4*(j+njobs*(itot*job->ainfo->nm+mi));
int phas_idx = itot*job->s_th + mi*job->s_m + 4*j;
dcmplx p1Q=job->phase[phas_idx],
p1U=job->phase[phas_idx+2],
p2Q=ispair[itot] ? job->phase[phas_idx+1]:0.,

1
libsharp/sharp_internal.h
View file

@ -48,6 +48,7 @@ typedef struct
int flags;
void **map;
void **alm;
int s_m, s_th; // strides in m and theta direction
complex double *phase;
double *norm_l;
complex double *almtmp;