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

tweak alm <-> alm_tmp

Browse source
This commit is contained in:
Martin Reinecke 2012-11-05 13:13:20 +01:00
parent bea64bd65f
commit 549d1c35e1
1 changed files with 56 additions and 19 deletions
Download patch file Download diff file Expand all files Collapse all files

75
libsharp/sharp.c
View file

@ -415,18 +415,36 @@ static void dealloc_almtmp (sharp_job *job)
static void alm2almtmp (sharp_job *job, int lmax, int mi) static void alm2almtmp (sharp_job *job, int lmax, int mi)
{ {
if (job->type!=SHARP_MAP2ALM) if (job->type!=SHARP_MAP2ALM)
for (int l=job->ainfo->mval[mi]; l<=lmax; ++l) {
ptrdiff_t ofs=job->ainfo->mvstart[mi];
int stride=job->ainfo->stride;
if (job->spin==0)
{ {
ptrdiff_t aidx = sharp_alm_index(job->ainfo,l,mi); if (job->fde==DOUBLE)
double fct = job->norm_l[l]; for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
for (int i=0; i<job->ntrans*job->nalm; ++i) for (int i=0; i<job->ntrans*job->nalm; ++i)
if (job->fde==DOUBLE) job->almtmp[job->ntrans*job->nalm*l+i]
job->almtmp[job->ntrans*job->nalm*l+i] = ((dcmplx *)job->alm[i])[ofs+l*stride];
= ((dcmplx *)job->alm[i])[aidx]*fct; else
else for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
job->almtmp[job->ntrans*job->nalm*l+i] for (int i=0; i<job->ntrans*job->nalm; ++i)
= ((fcmplx *)job->alm[i])[aidx]*fct; job->almtmp[job->ntrans*job->nalm*l+i]
= ((fcmplx *)job->alm[i])[ofs+l*stride];
} }
else
{
if (job->fde==DOUBLE)
for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
for (int i=0; i<job->ntrans*job->nalm; ++i)
job->almtmp[job->ntrans*job->nalm*l+i]
= ((dcmplx *)job->alm[i])[ofs+l*stride]*job->norm_l[l];
else
for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
for (int i=0; i<job->ntrans*job->nalm; ++i)
job->almtmp[job->ntrans*job->nalm*l+i]
= ((fcmplx *)job->alm[i])[ofs+l*stride]*job->norm_l[l];
}
}
else else
SET_ARRAY(job->almtmp,job->ntrans*job->nalm*job->ainfo->mval[mi], SET_ARRAY(job->almtmp,job->ntrans*job->nalm*job->ainfo->mval[mi],
job->ntrans*job->nalm*(lmax+1),0.); job->ntrans*job->nalm*(lmax+1),0.);
@ -435,16 +453,33 @@ static void alm2almtmp (sharp_job *job, int lmax, int mi)
static void almtmp2alm (sharp_job *job, int lmax, int mi) static void almtmp2alm (sharp_job *job, int lmax, int mi)
{ {
if (job->type != SHARP_MAP2ALM) return; if (job->type != SHARP_MAP2ALM) return;
for (int l=job->ainfo->mval[mi]; l<=lmax; ++l) ptrdiff_t ofs=job->ainfo->mvstart[mi];
int stride=job->ainfo->stride;
if (job->spin==0)
{ {
ptrdiff_t aidx = sharp_alm_index(job->ainfo,l,mi); if (job->fde==DOUBLE)
for (int i=0;i<job->ntrans*job->nalm;++i) for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
if (job->fde==DOUBLE) for (int i=0;i<job->ntrans*job->nalm;++i)
((dcmplx *)job->alm[i])[aidx] += ((dcmplx *)job->alm[i])[ofs+l*stride] +=
job->almtmp[job->ntrans*job->nalm*l+i]*job->norm_l[l]; job->almtmp[job->ntrans*job->nalm*l+i];
else else
((fcmplx *)job->alm[i])[aidx] += for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
(fcmplx)(job->almtmp[job->ntrans*job->nalm*l+i]*job->norm_l[l]); for (int i=0;i<job->ntrans*job->nalm;++i)
((fcmplx *)job->alm[i])[ofs+l*stride] +=
(fcmplx)(job->almtmp[job->ntrans*job->nalm*l+i]);
}
else
{
if (job->fde==DOUBLE)
for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
for (int i=0;i<job->ntrans*job->nalm;++i)
((dcmplx *)job->alm[i])[ofs+l*stride] +=
job->almtmp[job->ntrans*job->nalm*l+i]*job->norm_l[l];
else
for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
for (int i=0;i<job->ntrans*job->nalm;++i)
((fcmplx *)job->alm[i])[ofs+l*stride] +=
(fcmplx)(job->almtmp[job->ntrans*job->nalm*l+i]*job->norm_l[l]);
} }
} }
@ -600,6 +635,8 @@ static int sharp_oracle (sharp_jobtype type, int spin, int ntrans)
int nrings=(lmax+1)/4; int nrings=(lmax+1)/4;
int ppring=1; int ppring=1;
spin = (spin!=0) ? 2 : 0;
ptrdiff_t npix=(ptrdiff_t)nrings*ppring; ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
sharp_geom_info *tinfo; sharp_geom_info *tinfo;
sharp_make_gauss_geom_info (nrings, ppring, 1, ppring, &tinfo); sharp_make_gauss_geom_info (nrings, ppring, 1, ppring, &tinfo);