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

simplify the interface as much as possible; nicer interfaces

Browse source 
can be added again later
This commit is contained in:
Martin Reinecke 2012-09-17 15:35:38 +02:00
parent 6beb0e027d
commit c459e08b48
12 changed files with 189 additions and 196 deletions
Download patch file Download diff file Expand all files Collapse all files

74
libsharp/sharp.c
View file

@ -32,7 +32,7 @@
#include <math.h>
#include "ls_fft.h"
#include "sharp_ylmgen_c.h"
#include "sharp.h"
#include "sharp_internal.h"
#include "c_utils.h"
#include "sharp_core.h"
#include "sharp_vecutil.h"
@ -458,7 +458,7 @@ static void phase2map (sharp_job *job, int mmax, int llim, int ulim)
} /* end of parallel region */
}
void sharp_execute_job (sharp_job *job)
static void sharp_execute_job (sharp_job *job)
{
double timer=wallTime();
job->opcnt=0;
@ -541,9 +541,10 @@ void sharp_execute_job (sharp_job *job)
job->time=wallTime()-timer;
}
static void sharp_build_job_common (sharp_job *job, sharp_jobtype type, int spin,
int add_output, const sharp_geom_info *geom_info,
const sharp_alm_info *alm_info, int ntrans)
static void sharp_build_job_common (sharp_job *job, sharp_jobtype type,
int spin, int add_output, void **alm, void **map,
const sharp_geom_info *geom_info, const sharp_alm_info *alm_info, int ntrans,
int dp, int nv)
{
UTIL_ASSERT((ntrans>0),"bad number of simultaneous transforms");
if (type==SHARP_ALM2MAP_DERIV1) spin=1;
@ -556,47 +557,26 @@ static void sharp_build_job_common (sharp_job *job, sharp_jobtype type, int spin
job->nalm = (type==SHARP_ALM2MAP_DERIV1) ? 1 : ((spin>0) ? 2 : 1);
job->ginfo = geom_info;
job->ainfo = alm_info;
job->nv = sharp_nv_oracle (type, spin, ntrans);
job->nv = (nv==0) ? sharp_nv_oracle (type, spin, ntrans) : nv;
job->time = 0.;
job->opcnt = 0;
job->ntrans = ntrans;
job->alm=alm;
job->map=map;
job->fde=dp ? DOUBLE : FLOAT;
}
void sharpd_build_job (sharp_job *job, sharp_jobtype type, int spin,
int add_output, dcmplx **alm, double **map, const sharp_geom_info *geom_info,
const sharp_alm_info *alm_info, int ntrans)
{
sharp_build_job_common (job, type, spin, add_output, geom_info, alm_info,
ntrans);
job->alm=(void **)alm;
job->map=(void **)map;
job->fde=DOUBLE;
}
void sharps_build_job (sharp_job *job, sharp_jobtype type, int spin,
int add_output, fcmplx **alm, float **map, const sharp_geom_info *geom_info,
const sharp_alm_info *alm_info, int ntrans)
{
sharp_build_job_common (job, type, spin, add_output, geom_info, alm_info,
ntrans);
job->alm=(void **)alm;
job->map=(void **)map;
job->fde=FLOAT;
}
void sharp_execute_ll (sharp_jobtype type, int spin, int add_output, void **alm,
void sharp_execute (sharp_jobtype type, int spin, int add_output, void **alm,
void **map, const sharp_geom_info *geom_info, const sharp_alm_info *alm_info,
int ntrans, int dp)
int ntrans, int dp, int nv, double *time, unsigned long long *opcnt)
{
sharp_job job;
if (dp)
sharpd_build_job (&job, type, spin, add_output, (dcmplx **)alm,
(double **)map, geom_info, alm_info, ntrans);
else
sharps_build_job (&job, type, spin, add_output, (fcmplx **)alm,
(float **)map, geom_info, alm_info, ntrans);
sharp_build_job_common (&job, type, spin, add_output, alm, map, geom_info,
alm_info, ntrans, dp, nv);
sharp_execute_job (&job);
if (time!=NULL) *time = job.time;
if (opcnt!=NULL) *opcnt = job.opcnt;
}
int sharp_get_nv_max (void)
@ -633,17 +613,18 @@ static int sharp_oracle (sharp_jobtype type, int spin, int ntrans)
for (int nv=1; nv<=sharp_get_nv_max(); ++nv)
{
double time_acc=0.;
sharp_job job;
sharpd_build_job(&job,type,spin,0,&alm[0],&map[0],tinfo,alms,ntrans);
job.nv=nv;
double jtime;
int ntries=0;
do
{
sharp_execute_job(&job);
sharp_execute(type,spin,0,(void **)(&alm[0]),(void **)(&map[0]),tinfo,
alms,ntrans,1,nv,&jtime,NULL);
if (job.time<time) { time=job.time; nvbest=nv; }
time_acc+=job.time;
if (jtime<time) { time=jtime; nvbest=nv; }
time_acc+=jtime;
++ntries;
}
while (time_acc<0.02);
while ((time_acc<0.02)&&(ntries<2));
}
DEALLOC2D(map);
@ -665,20 +646,13 @@ int sharp_nv_oracle (sharp_jobtype type, int spin, int ntrans)
{{0,0,0},{0,0,0}},
{{0,0,0},{0,0,0}} };
static int in_oracle=0;
if (in_oracle) return -20;
if (type==SHARP_ALM2MAP_DERIV1) spin=1;
UTIL_ASSERT((ntrans>0),"bad number of simultaneous transforms");
UTIL_ASSERT((spin>=0)&&(spin<=30), "bad spin");
ntrans=IMIN(ntrans,maxtr);
if (nv_opt[ntrans-1][spin!=0][type]==0)
{
in_oracle=1;
nv_opt[ntrans-1][spin!=0][type]=sharp_oracle(type,spin,ntrans);
in_oracle=0;
}
return nv_opt[ntrans-1][spin!=0][type];
}

64
libsharp/sharp.h
View file

@ -40,68 +40,4 @@
#include "sharp_lowlevel.h"
typedef enum { FLOAT, DOUBLE } sharp_fde;
/*! \internal
Type holding all required information about an SHT job. */
typedef struct
{
sharp_jobtype type;
int spin;
int add_output;
int nmaps, nalm;
sharp_fde fde;
void **map;
void **alm;
complex double *phase;
double *norm_l;
complex double *almtmp;
const sharp_geom_info *ginfo;
const sharp_alm_info *ainfo;
int nv;
double time;
int ntrans;
unsigned long long opcnt;
} sharp_job;
/*! \defgroup jobgroup Functionality for defining and executing SHTs */
/*! \{ */
/*! Initializes \a job with the appropriate parameters to perform the required
SHT.
\param type the type of SHT
\param spin the spin of the quantities to be transformed
\param add_output if 0, the output arrays will be overwritten,
else the result will be added to the output arrays.
\param ntrans the number of simultaneous SHTs
\param alm contains pointers to the a_lm coefficients. If \a spin==0,
alm[0] points to the a_lm of the first SHT, alm[1] to those of the second
etc. If \a spin>0, alm[0] and alm[1] point to the a_lm of the first SHT,
alm[2] and alm[3] to those of the second, etc.
\param map contains pointers to the maps. If \a spin==0,
map[0] points to the map of the first SHT, map[1] to that of the second
etc. If \a spin>0, map[0] and map[1] point to the maps of the first SHT,
map[2] and map[3] to those of the second, etc.
\note \a map and \a a_lm must not be de-allocated until after the last call of
sharp_execute_job()! This is because the library does not copy the input
data, but only stores the pointers to the supplied maps and a_lm. */
void sharpd_build_job (sharp_job *job, sharp_jobtype type, int spin,
int add_output, complex double **alm, double **map,
const sharp_geom_info *geom_info, const sharp_alm_info *alm_info, int ntrans);
void sharps_build_job (sharp_job *job, sharp_jobtype type, int spin,
int add_output, complex float **alm, float **map,
const sharp_geom_info *geom_info, const sharp_alm_info *alm_info, int ntrans);
/*! Execute the SHT job previously constructed by sharpd_build_job() or
sharps_build_job(). */
void sharp_execute_job (sharp_job *job);
/*! \} */
/*! Internal */
int sharp_get_nv_max (void);
/*! Internal */
int sharp_nv_oracle (sharp_jobtype type, int spin, int ntrans);
#endif

30
libsharp/sharp_acctest.c
View file

@ -111,9 +111,8 @@ static void check_sign_scale(void)
for (int j=0; j<nalms; ++j)
alm[i][j]=1.+_Complex_I;
sharp_job job;
sharpd_build_job(&job,SHARP_ALM2MAP,0,0,&alm[0],&map[0],tinfo,alms,ntrans);
sharp_execute_job(&job);
sharp_execute(SHARP_ALM2MAP,0,0,(void **)&alm[0],(void **)&map[0],tinfo,
alms,ntrans,1,0,NULL,NULL);
for (int it=0; it<ntrans; ++it)
{
UTIL_ASSERT(FAPPROX(map[it][0 ], 3.588246976618616912e+00,1e-12),
@ -123,8 +122,8 @@ static void check_sign_scale(void)
UTIL_ASSERT(FAPPROX(map[it][npix-1],-1.234675107554816442e+01,1e-12),
"error");
}
sharpd_build_job(&job,SHARP_ALM2MAP,1,0,&alm[0],&map[0],tinfo,alms,ntrans);
sharp_execute_job(&job);
sharp_execute(SHARP_ALM2MAP,1,0,(void **)&alm[0],(void **)&map[0],tinfo,
alms,ntrans,1,0,NULL,NULL);
for (int it=0; it<ntrans; ++it)
{
UTIL_ASSERT(FAPPROX(map[2*it ][0 ], 2.750897760535633285e+00,1e-12),
@ -141,8 +140,8 @@ static void check_sign_scale(void)
"error");
}
sharpd_build_job(&job,SHARP_ALM2MAP,2,0,&alm[0],&map[0],tinfo,alms,ntrans);
sharp_execute_job(&job);
sharp_execute(SHARP_ALM2MAP,2,0,(void **)&alm[0],(void **)&map[0],tinfo,
alms,ntrans,1,0,NULL,NULL);
for (int it=0; it<ntrans; ++it)
{
UTIL_ASSERT(FAPPROX(map[2*it ][0 ],-1.398186224727334448e+00,1e-12),
@ -159,9 +158,8 @@ static void check_sign_scale(void)
"error");
}
sharpd_build_job(&job,SHARP_ALM2MAP_DERIV1,1,0,&alm[0],&map[0],tinfo,alms,
ntrans);
sharp_execute_job(&job);
sharp_execute(SHARP_ALM2MAP_DERIV1,1,0,(void **)&alm[0],(void **)&map[0],
tinfo,alms,ntrans,1,0,NULL,NULL);
for (int it=0; it<ntrans; ++it)
{
UTIL_ASSERT(FAPPROX(map[2*it ][0 ],-6.859393905369091105e-01,1e-11),
@ -207,14 +205,10 @@ static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
dcmplx **alm2;
ALLOC2D(alm2,dcmplx,ncomp,nalms);
sharp_job job;
sharpd_build_job(&job,SHARP_ALM2MAP,spin,0,&alm[0],&map[0],tinfo,alms,ntrans);
job.nv=nv;
sharp_execute_job(&job);
sharpd_build_job(&job,SHARP_MAP2ALM,spin,0,&alm2[0],&map[0],tinfo,alms,ntrans);
job.nv=nv;
sharp_execute_job(&job);
sharp_execute(SHARP_ALM2MAP,spin,0,(void **)(&alm[0]),(void **)(&map[0]),
tinfo,alms,ntrans,1,nv,NULL,NULL);
sharp_execute(SHARP_MAP2ALM,spin,0,(void **)(&alm2[0]),(void **)(&map[0]),
tinfo,alms,ntrans,1,nv,NULL,NULL);
measure_errors(alm,alm2,nalms,ncomp);
DEALLOC2D(map);

14
libsharp/sharp_bench.c
View file

@ -67,19 +67,17 @@ static void bench_sht (int spin, int nv, sharp_jobtype type,
SET_ARRAY(alm[0],0,nalms*ncomp,0.);
int nruns=0;
sharp_job job;
sharpd_build_job(&job,type,spin,0,&alm[0],&map[0],tinfo,alms,ntrans);
job.nv=nv;
*time=1e30;
*opcnt=1000000000000000;
do
{
sharpd_build_job(&job,type,spin,0,&alm[0],&map[0],tinfo,alms,ntrans);
job.nv=nv;
sharp_execute_job(&job);
double jtime;
unsigned long long jopcnt;
sharp_execute(type,spin,0,(void **)(&alm[0]),(void **)(&map[0]),tinfo,alms,
ntrans,1,nv,&jtime,&jopcnt);
if (job.opcnt<*opcnt) *opcnt=job.opcnt;
if (job.time<*time) *time=job.time;
if (jopcnt<*opcnt) *opcnt=jopcnt;
if (jtime<*time) *time=jtime;
}
while (++nruns < 4);

2
libsharp/sharp_core.h
View file

@ -32,7 +32,7 @@
#ifndef PLANCK_SHARP_CORE_H
#define PLANCK_SHARP_CORE_H
#include "sharp.h"
#include "sharp_internal.h"
#include "sharp_ylmgen_c.h"
#ifdef __cplusplus

20
libsharp/sharp_cxx.h
View file

@ -98,8 +98,8 @@ template<typename T> class sharp_cxxjob: public sharp_base
void alm2map (const xcomplex<T> *alm, T *map, bool add)
{
void *aptr=conv(alm), *mptr=conv(map);
sharp_execute_ll (SHARP_ALM2MAP, 0, add, &aptr, &mptr, ginfo, ainfo, 1,
cxxjobhelper__<T>::val);
sharp_execute (SHARP_ALM2MAP, 0, add, &aptr, &mptr, ginfo, ainfo, 1,
cxxjobhelper__<T>::val,0,0,0);
}
void alm2map_spin (const xcomplex<T> *alm1, const xcomplex<T> *alm2,
T *map1, T *map2, int spin, bool add)
@ -107,21 +107,21 @@ template<typename T> class sharp_cxxjob: public sharp_base
void *aptr[2], *mptr[2];
aptr[0]=conv(alm1); aptr[1]=conv(alm2);
mptr[0]=conv(map1); mptr[1]=conv(map2);
sharp_execute_ll (SHARP_ALM2MAP, spin, add, aptr, mptr, ginfo, ainfo, 1,
cxxjobhelper__<T>::val);
sharp_execute (SHARP_ALM2MAP, spin, add, aptr, mptr, ginfo, ainfo, 1,
cxxjobhelper__<T>::val,0,0,0);
}
void alm2map_der1 (const xcomplex<T> *alm, T *map1, T *map2, bool add)
{
void *aptr=conv(alm), *mptr[2];
mptr[0]=conv(map1); mptr[1]=conv(map2);
sharp_execute_ll (SHARP_ALM2MAP_DERIV1, 1, add,&aptr, mptr, ginfo, ainfo,
1, cxxjobhelper__<T>::val);
sharp_execute (SHARP_ALM2MAP_DERIV1, 1, add,&aptr, mptr, ginfo, ainfo,
1, cxxjobhelper__<T>::val,0,0,0);
}
void map2alm (const T *map, xcomplex<T> *alm, bool add)
{
void *aptr=conv(alm), *mptr=conv(map);
sharp_execute_ll (SHARP_MAP2ALM, 0, add, &aptr, &mptr, ginfo, ainfo, 1,
cxxjobhelper__<T>::val);
sharp_execute (SHARP_MAP2ALM, 0, add, &aptr, &mptr, ginfo, ainfo, 1,
cxxjobhelper__<T>::val,0,0,0);
}
void map2alm_spin (const T *map1, const T *map2, xcomplex<T> *alm1,
xcomplex<T> *alm2, int spin, bool add)
@ -129,8 +129,8 @@ template<typename T> class sharp_cxxjob: public sharp_base
void *aptr[2], *mptr[2];
aptr[0]=conv(alm1); aptr[1]=conv(alm2);
mptr[0]=conv(map1); mptr[1]=conv(map2);
sharp_execute_ll (SHARP_MAP2ALM, spin, add, aptr, mptr, ginfo, ainfo, 1,
cxxjobhelper__<T>::val);
sharp_execute (SHARP_MAP2ALM, spin, add, aptr, mptr, ginfo, ainfo, 1,
cxxjobhelper__<T>::val,0,0,0);
}
};

66
libsharp/sharp_internal.h Normal file
View file

@ -0,0 +1,66 @@
/*
* This file is part of libsharp.
*
* libsharp is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libsharp is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libsharp; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libsharp is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*! \file sharp_internal.h
* Internally used functionality for the spherical transform library.
*
* Copyright (C) 2006-2012 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef PLANCK_SHARP_INTERNAL_H
#define PLANCK_SHARP_INTERNAL_H
#ifdef __cplusplus
#error This header file cannot be included from C++, only from C
#endif
#include "sharp.h"
typedef enum { FLOAT, DOUBLE } sharp_fde;
typedef struct
{
sharp_jobtype type;
int spin;
int add_output;
int nmaps, nalm;
sharp_fde fde;
void **map;
void **alm;
complex double *phase;
double *norm_l;
complex double *almtmp;
const sharp_geom_info *ginfo;
const sharp_alm_info *ainfo;
int nv;
double time;
int ntrans;
unsigned long long opcnt;
} sharp_job;
int sharp_get_nv_max (void);
int sharp_nv_oracle (sharp_jobtype type, int spin, int ntrans);
#endif

12
libsharp/sharp_lowlevel.h
View file

@ -163,10 +163,16 @@ typedef enum { SHARP_MAP2ALM, /*!< analysis */
\param ntrans the number of simultaneous SHTs
\param dp if 0, the \a alm is expected to have the type "complex float **"
and \a map is expected to have the type "float **"; otherwise the expected
types are "complex double **" and "double **", respectively. */
void sharp_execute_ll (sharp_jobtype type, int spin, int add_output, void **alm,
types are "complex double **" and "double **", respectively.
\param nv Internally used SHT parameter. Set to 0 unless you know what you are
doing.
\param time If not NULL, the wall clock time required for this SHT
(in seconds)will be written here.
\param opcnt If not NULL, a conservative estimate of the total floating point
operation count for this SHT will be written here. */
void sharp_execute (sharp_jobtype type, int spin, int add_output, void **alm,
void **map, const sharp_geom_info *geom_info, const sharp_alm_info *alm_info,
int ntrans, int dp);
int ntrans, int dp, int nv, double *time, unsigned long long *opcnt);
/*! \} */

16
libsharp/sharp_mpi.c
View file

@ -199,7 +199,7 @@ static void map2alm_comm (sharp_job *job, const sharp_mpi_info *minfo)
sharp_communicate_map2alm (minfo,&job->phase);
}
void sharp_execute_job_mpi (sharp_job *job, MPI_Comm comm)
static void sharp_execute_job_mpi (sharp_job *job, MPI_Comm comm)
{
double timer=wallTime();
int ntasks;
@ -283,4 +283,18 @@ void sharp_execute_job_mpi (sharp_job *job, MPI_Comm comm)
job->time=wallTime()-timer;
}
void sharp_execute_mpi (MPI_Comm comm, sharp_jobtype type, int spin,
int add_output, void **alm, void **map, const sharp_geom_info *geom_info,
const sharp_alm_info *alm_info, int ntrans, int dp, int nv, double *time,
unsigned long long *opcnt)
{
sharp_job job;
sharp_build_job_common (&job, type, spin, add_output, alm, map, geom_info,
alm_info, ntrans, dp, nv);
sharp_execute_job_mpi (&job, comm);
if (time!=NULL) *time = job.time;
if (opcnt!=NULL) *opcnt = job.opcnt;
}
#endif

5
libsharp/sharp_mpi.h
View file

@ -39,7 +39,10 @@
extern "C" {
#endif
void sharp_execute_job_mpi (sharp_job *job, MPI_Comm comm);
void sharp_execute_mpi (MPI_Comm comm, sharp_jobtype type, int spin,
int add_output, void **alm, void **map, const sharp_geom_info *geom_info,
const sharp_alm_info *alm_info, int ntrans, int dp, int nv, double *time,
unsigned long long *opcnt);
#ifdef __cplusplus
}

40
libsharp/sharp_test.c
View file

@ -116,11 +116,12 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
sharp_alm_info *alms;
sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
sharp_job job;
sharpd_build_job(&job,SHARP_MAP2ALM,spin,0,&alm[0],&map[0],tinfo,alms,ntrans);
sharp_execute_job(&job);
printf("wall time for map2alm: %fs\n",job.time);
printf("Performance: %fGFLOPs/s\n",1e-9*job.opcnt/job.time);
double time;
unsigned long long opcnt;
sharp_execute(SHARP_MAP2ALM,spin,0,(void **)&alm[0],(void **)&map[0],tinfo,alms,ntrans,1,0,
&time,&opcnt);
printf("wall time for map2alm: %fs\n",time);
printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);
measure_errors(alm_orig,alm,nalms,ncomp);
for (int iter=0; iter<niter; ++iter)
@ -128,20 +129,18 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
double **map2;
ALLOC2D(map2,double,ncomp,npix);
printf ("\niteration %i:\n", iter+1);
sharpd_build_job(&job,SHARP_ALM2MAP,spin,0,&alm[0],&map2[0],tinfo,alms,
ntrans);
sharp_execute_job(&job);
printf("wall time for alm2map: %fs\n",job.time);
printf("Performance: %fGFLOPs/s\n",1e-9*job.opcnt/job.time);
sharp_execute(SHARP_ALM2MAP,spin,0,(void **)&alm[0],(void **)&map2[0],tinfo,alms,ntrans,1,0,
&time,&opcnt);
printf("wall time for alm2map: %fs\n",time);
printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);
for (int i=0; i<ncomp; ++i)
for (ptrdiff_t m=0; m<npix; ++m)
map2[i][m] = map[i][m]-map2[i][m];
sharpd_build_job(&job,SHARP_MAP2ALM,spin,1,&alm[0],&map2[0],tinfo,alms,
ntrans);
sharp_execute_job(&job);
printf("wall time for map2alm: %fs\n",job.time);
printf("Performance: %fGFLOPs/s\n",1e-9*job.opcnt/job.time);
sharp_execute(SHARP_MAP2ALM,spin,1,(void **)&alm[0],(void **)&map2[0],tinfo,alms,ntrans,1,0,
&time,&opcnt);
printf("wall time for map2alm: %fs\n",time);
printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);
DEALLOC2D(map2);
measure_errors(alm_orig,alm,nalms,ncomp);
}
@ -170,12 +169,13 @@ static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
dcmplx **alm2;
ALLOC2D(alm2,dcmplx,ncomp,nalms);
sharp_job job;
double time;
unsigned long long opcnt;
printf ("\niteration 0:\n");
sharpd_build_job(&job,SHARP_ALM2MAP,spin,0,&alm[0],&map[0],tinfo,alms,ntrans);
sharp_execute_job(&job);
printf("wall time for alm2map: %fs\n",job.time);
printf("Performance: %fGFLOPs/s\n",1e-9*job.opcnt/job.time);
sharp_execute(SHARP_ALM2MAP,spin,0,(void **)&alm[0],(void **)&map[0],tinfo,alms,ntrans,1,0,
&time,&opcnt);
printf("wall time for alm2map: %fs\n",time);
printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);
map2alm_iter(tinfo,map,alm,alm2,lmax,mmax,npix,nalms,spin,ntrans,niter);

42
libsharp/sharp_test_mpi.c
View file

@ -197,11 +197,13 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
reduce_alm_info(alms);
sharp_job job;
sharpd_build_job(&job,SHARP_MAP2ALM,spin,0,&alm[0],&map[0],tinfo,alms,ntrans);
sharp_execute_job_mpi(&job,MPI_COMM_WORLD);
unsigned long long opcnt=totalops(job.opcnt);
double timer=maxTime(job.time);
double jtime;
unsigned long long jopcnt;
sharp_execute_mpi(MPI_COMM_WORLD,SHARP_MAP2ALM,spin,0,(void **)&alm[0],
(void **)&map[0],tinfo,alms,ntrans,1,0,&jtime,&jopcnt);
unsigned long long opcnt=totalops(jopcnt);
double timer=maxTime(jtime);
if (mytask==0) printf("wall time for map2alm: %fs\n",timer);
if (mytask==0) printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/timer);
measure_errors(alm_orig,alm,alms,ncomp);
@ -211,22 +213,20 @@ static void map2alm_iter (sharp_geom_info *tinfo, double **map,
double **map2;
ALLOC2D(map2,double,ncomp,npix);
if (mytask==0) printf ("\niteration %i:\n", iter+1);
sharpd_build_job(&job,SHARP_ALM2MAP,spin,0,&alm[0],&map2[0],tinfo,alms,
ntrans);
sharp_execute_job_mpi(&job,MPI_COMM_WORLD);
opcnt=totalops(job.opcnt);
timer=maxTime(job.time);
sharp_execute_mpi(MPI_COMM_WORLD,SHARP_ALM2MAP,spin,0,(void **)&alm[0],
(void **)&map2[0],tinfo,alms,ntrans,1,0,&jtime,&jopcnt);
opcnt=totalops(jopcnt);
timer=maxTime(jtime);
if (mytask==0) printf("wall time for alm2map: %fs\n",timer);
if (mytask==0) printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/timer);
for (int i=0; i<ncomp; ++i)
for (ptrdiff_t m=0; m<npix; ++m)
map2[i][m] = map[i][m]-map2[i][m];
sharpd_build_job(&job,SHARP_MAP2ALM,spin,1,&alm[0],&map2[0],tinfo,alms,
ntrans);
sharp_execute_job_mpi(&job,MPI_COMM_WORLD);
opcnt=totalops(job.opcnt);
timer=maxTime(job.time);
sharp_execute_mpi(MPI_COMM_WORLD,SHARP_MAP2ALM,spin,1,(void **)&alm[0],
(void **)&map2[0],tinfo,alms,ntrans,1,0,&jtime,&jopcnt);
opcnt=totalops(jopcnt);
timer=maxTime(jtime);
if (mytask==0) printf("wall time for map2alm: %fs\n",wallTime()-timer);
if (mytask==0) printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/timer);
DEALLOC2D(map2);
@ -244,6 +244,9 @@ static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
double **map;
ALLOC2D(map,double,ncomp,npix);
double jtime;
unsigned long long jopcnt;
sharp_alm_info *alms;
ptrdiff_t nalms;
sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
@ -260,11 +263,10 @@ static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
ALLOC2D(alm2,dcmplx,ncomp,nalms);
if (mytask==0) printf ("\niteration 0:\n");
sharp_job job;
sharpd_build_job(&job,SHARP_ALM2MAP,spin,0,&alm[0],&map[0],tinfo,alms,ntrans);
sharp_execute_job_mpi(&job,MPI_COMM_WORLD);
unsigned long long opcnt=totalops(job.opcnt);
double timer=maxTime(job.time);
sharp_execute_mpi(MPI_COMM_WORLD,SHARP_ALM2MAP,spin,0,(void **)&alm[0],
(void **)&map[0],tinfo,alms,ntrans,1,0,&jtime,&jopcnt);
unsigned long long opcnt=totalops(jopcnt);
double timer=maxTime(jtime);
if (mytask==0) printf("wall time for alm2map: %fs\n",timer);
if (mytask==0) printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/timer);