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

improve mlim computation

Browse source
This commit is contained in:
Martin Reinecke 2013-01-22 14:54:25 +01:00
parent ed0a0c9394
commit 4b3da2e693
4 changed files with 9 additions and 19 deletions
Download patch file Download diff file Expand all files Collapse all files

9
libsharp/sharp.c
View file

@ -63,9 +63,10 @@ static void get_chunk_info (int ndata, int nmult, int *nchunks, int *chunksize)
*nchunks = (ndata+(*chunksize)-1)/(*chunksize);
}
static int sharp_get_mlim (int lmax, int spin, double sth, double cth,
double ofs)
int sharp_get_mlim (int lmax, int spin, double sth, double cth)
{
double ofs=lmax*0.01;
if (ofs<100.) ofs=100.;
double b = -2*spin*fabs(cth);
double t1 = lmax*sth+ofs;
double c = (double)spin*spin-t1*t1;
@ -651,7 +652,7 @@ static void sharp_execute_job (sharp_job *job)
ispair[i] = job->ginfo->pair[i+llim].r2.nph>0;
cth[i] = job->ginfo->pair[i+llim].r1.cth;
sth[i] = job->ginfo->pair[i+llim].r1.sth;
mlim[i] = sharp_get_mlim(lmax, job->spin, sth[i], cth[i], 100.);
mlim[i] = sharp_get_mlim(lmax, job->spin, sth[i], cth[i]);
}
/* map->phase where necessary */
@ -750,7 +751,7 @@ int sharp_get_nv_max (void)
static int sharp_oracle (sharp_jobtype type, int spin, int ntrans)
{
int lmax=127;
int lmax=511;
int mmax=(lmax+1)/2;
int nrings=(lmax+1)/4;
int ppring=1;

1
libsharp/sharp_internal.h
View file

@ -61,5 +61,6 @@ typedef struct
int sharp_get_nv_max (void);
int sharp_nv_oracle (sharp_jobtype type, int spin, int ntrans);
int sharp_get_mlim (int lmax, int spin, double sth, double cth);
#endif

2
libsharp/sharp_mpi.c
View file

@ -261,7 +261,7 @@ static void sharp_execute_job_mpi (sharp_job *job, MPI_Comm comm)
{
cth[i] = cos(minfo.theta[i]);
sth[i] = sin(minfo.theta[i]);
mlim[i] = sharp_get_mlim(lmax, job->spin, sth[i], cth[i], 100.);
mlim[i] = sharp_get_mlim(lmax, job->spin, sth[i], cth[i]);
}
/* map->phase where necessary */

16
libsharp/sharp_testsuite.c
View file

@ -38,6 +38,7 @@
#include <omp.h>
#endif
#include "sharp.h"
#include "sharp_internal.h"
#include "sharp_geomhelpers.h"
#include "sharp_almhelpers.h"
#include "c_utils.h"
@ -243,19 +244,6 @@ static int good_fft_size(int n)
return bestfac;
}
static int sharp_get_mlim (int lmax, int spin, double sth, double cth,
double ofs)
{
double b = -2*spin*fabs(cth);
double t1 = lmax*sth+ofs;
double c = spin*spin-t1*t1;
double discr = b*b-4*c;
if (discr<=0) return lmax;
double res=(-b+sqrt(discr))/2.;
if (res>lmax) res=lmax;
return (int)(res+0.5);
}
static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
int *gpar2, sharp_geom_info **ginfo, sharp_alm_info **ainfo)
{
@ -319,7 +307,7 @@ static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1,
for (int i=0; i<(*ginfo)->npairs; ++i)
{
sharp_ringpair *pair=&((*ginfo)->pair[i]);
int pring=1+2*sharp_get_mlim(lmax,0,pair->r1.sth,pair->r1.cth,100.);
int pring=1+2*sharp_get_mlim(lmax,0,pair->r1.sth,pair->r1.cth);
if (pring>nlon) pring=nlon;
pring=good_fft_size(pring);
pair->r1.nph=pring;