diff --git a/libsharp/libsharp.dox b/libsharp/libsharp.dox
index c0c1a59..4c441f1 100644
--- a/libsharp/libsharp.dox
+++ b/libsharp/libsharp.dox
@@ -82,12 +82,4 @@
single-precision transforms will most likely not be faster than their
double-precision counterparts, but they will require significantly less
memory.
-
- Two example and benchmark programs are distributed with libsharp:
-
- - sharp_test.c checks the accuracy of the (iterative) map analysis
- algorithm
- - sharp_bench.c determines the quickest transform strategy for a given
- SHT
-
*/
diff --git a/libsharp/sharp_bench.c b/libsharp/sharp_bench.c
deleted file mode 100644
index e7f8f23..0000000
--- a/libsharp/sharp_bench.c
+++ /dev/null
@@ -1,149 +0,0 @@
-/*
- * 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_bench.c
- Copyright (C) 2012 Max-Planck-Society
- \author Martin Reinecke
-*/
-
-#include
-#include
-#ifdef USE_MPI
-#include "mpi.h"
-#endif
-#include "sharp.h"
-#include "sharp_geomhelpers.h"
-#include "sharp_almhelpers.h"
-#include "c_utils.h"
-#include "sharp_announce.h"
-#include "sharp_core.h"
-
-typedef complex double dcmplx;
-
-static void bench_sht (int spin, int nv, sharp_jobtype type,
- int ntrans, double *time, unsigned long long *opcnt)
- {
- int lmax=2047;
- int mmax=128;
- int nrings=512;
- int ppring=1024;
- ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
- sharp_geom_info *tinfo;
- sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-
- ptrdiff_t nalms = ((mmax+1)*(mmax+2))/2 + (mmax+1)*(lmax-mmax);
- int ncomp = ntrans*((spin==0) ? 1 : 2);
-
- double **map;
- ALLOC2D(map,double,ncomp,npix);
- SET_ARRAY(map[0],0,npix*ncomp,0.);
-
- sharp_alm_info *alms;
- sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
-
- dcmplx **alm;
- ALLOC2D(alm,dcmplx,ncomp,nalms);
- SET_ARRAY(alm[0],0,nalms*ncomp,0.);
-
- int nruns=0;
- *time=1e30;
- *opcnt=1000000000000000;
- do
- {
- double jtime;
- unsigned long long jopcnt;
- sharp_execute(type,spin,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP|nv,
- &jtime,&jopcnt);
-
- if (jopcnt<*opcnt) *opcnt=jopcnt;
- if (jtime<*time) *time=jtime;
- }
- while (++nruns < 4);
-
- DEALLOC2D(map);
- DEALLOC2D(alm);
-
- sharp_destroy_alm_info(alms);
- sharp_destroy_geom_info(tinfo);
- }
-
-int main(void)
- {
-#ifdef USE_MPI
- MPI_Init(NULL,NULL);
-#endif
- sharp_module_startup("sharp_bench",1,1,"",1);
-
- printf("Benchmarking SHTs.\n\n");
- FILE *fp=fopen("sharp_oracle.inc","w");
- UTIL_ASSERT(fp, "failed to open oracle file for writing");
- fprintf(fp,"static const int maxtr = 6;\n");
- fprintf(fp,"static const int nv_opt[6][2][3] = {\n");
-
- const char *shtname[]={"map2alm","alm2map","a2mder1"};
-
- for (int ntr=1; ntr<=6; ++ntr)
- {
- fprintf(fp,"{");
- for (int spin=0; spin<=2; spin+=2)
- {
- fprintf(fp,"{");
- for (sharp_jobtype type=SHARP_MAP2ALM; type<=SHARP_ALM2MAP_DERIV1; ++type)
- {
- if ((type==SHARP_ALM2MAP_DERIV1) && (spin==0))
- fprintf(fp,"-1");
- else
- {
- int nvbest=-1, nvoracle=sharp_nv_oracle(type,spin,ntr);
- unsigned long long opmin=1000000000000000, op;
- double tmin=1e30;
- double *time=RALLOC(double,sharp_get_nv_max()+1);
- for (int nv=1; nv<=sharp_get_nv_max(); ++nv)
- {
- bench_sht (spin,nv,type,ntr,&time[nv],&op);
- if (op
-
-#if (defined(_OPENMP) && defined(USE_MPI))
-
-#include
-#include
-#include
-#include
-#include "sharp_mpi.h"
-#include "sharp.h"
-#include "sharp_vecutil.h"
-#include "sharp_geomhelpers.h"
-#include "sharp_almhelpers.h"
-#include "c_utils.h"
-#include "sharp_announce.h"
-#include "sharp_core.h"
-#include "memusage.h"
-
-typedef complex double dcmplx;
-
-int ntasks, mytask;
-
-static unsigned long long totalops (unsigned long long val)
- {
- unsigned long long tmp;
- MPI_Allreduce (&val, &tmp,1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
- return tmp;
- }
-
-static double maxTime (double val)
- {
- double tmp;
- MPI_Allreduce (&val, &tmp,1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
- return tmp;
- }
-
-static double totalMem (double val)
- {
- double tmp;
- MPI_Allreduce (&val, &tmp,1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
- return tmp;
- }
-
-static void reduce_alm_info(sharp_alm_info *ainfo)
- {
- int nmnew=0;
- ptrdiff_t ofs = 0;
- for (int i=mytask; inm; i+=ntasks,++nmnew)
- {
- ainfo->mval[nmnew]=ainfo->mval[i];
- ainfo->mvstart[nmnew]=ofs-ainfo->mval[nmnew];
- ofs+=ainfo->lmax-ainfo->mval[nmnew]+1;
- }
- ainfo->nm=nmnew;
- }
-
-static void reduce_geom_info(sharp_geom_info *ginfo)
- {
- int npairsnew=0;
- ptrdiff_t ofs = 0;
- for (int i=mytask; inpairs; i+=ntasks,++npairsnew)
- {
- ginfo->pair[npairsnew]=ginfo->pair[i];
- ginfo->pair[npairsnew].r1.ofs=ofs;
- ofs+=ginfo->pair[npairsnew].r1.nph;
- ginfo->pair[npairsnew].r2.ofs=ofs;
- if (ginfo->pair[npairsnew].r2.nph>0) ofs+=ginfo->pair[npairsnew].r2.nph;
- }
- ginfo->npairs=npairsnew;
- }
-
-static ptrdiff_t get_nalms(const sharp_alm_info *ainfo)
- {
- ptrdiff_t res=0;
- for (int i=0; inm; ++i)
- res += ainfo->lmax-ainfo->mval[i]+1;
- return res;
- }
-
-static ptrdiff_t get_npix(const sharp_geom_info *ginfo)
- {
- ptrdiff_t res=0;
- for (int i=0; inpairs; ++i)
- {
- res += ginfo->pair[i].r1.nph;
- if (ginfo->pair[i].r2.nph>0) res += ginfo->pair[i].r2.nph;
- }
- return res;
- }
-
-int main(int argc, char **argv)
- {
- MPI_Init(NULL,NULL);
- MPI_Comm_size(MPI_COMM_WORLD,&ntasks);
- MPI_Comm_rank(MPI_COMM_WORLD,&mytask);
- int master=(mytask==0);
-
- sharp_module_startup("sharp_bench2",argc,7,
- " ",0);
-
- int lmax=atoi(argv[2]);
- sharp_jobtype jtype = (strcmp(argv[4],"a2m")==0) ?
- SHARP_ALM2MAP : SHARP_MAP2ALM;
- int spin=atoi(argv[5]);
- int ntrans=atoi(argv[6]);
-
- sharp_geom_info *tinfo;
- ptrdiff_t npix=0;
- int geom2=0;
- if (strcmp(argv[1],"gauss")==0)
- {
- int nrings=geom2=lmax+1;
- int ppring=atoi(argv[3]);
- sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
- }
- else if (strcmp(argv[1],"ecp")==0)
- {
- int nrings=geom2=2*lmax+2;
- int ppring=atoi(argv[3]);
- sharp_make_ecp_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
- }
- else if (strcmp(argv[1],"healpix")==0)
- {
- int nside=atoi(argv[3]);
- if (nside<1) nside=1;
- geom2=4*nside-1;
- sharp_make_healpix_geom_info (nside, 1, &tinfo);
- }
- else
- UTIL_FAIL("unknown grid geometry");
-
- reduce_geom_info(tinfo);
- npix=get_npix(tinfo);
-
- int mmax=lmax;
- int ncomp = ntrans*((spin==0) ? 1 : 2);
-
- double **map;
- ALLOC2D(map,double,ncomp,npix);
-
- sharp_alm_info *alms;
- sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
-
- reduce_alm_info(alms);
- ptrdiff_t nalms=get_nalms(alms);
-
- dcmplx **alm;
- ALLOC2D(alm,dcmplx,ncomp,nalms);
-
- for (int n=0; n