From 78173bd83ec63c0703b288025374b2478ce6cfb8 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 21 Dec 2012 22:02:15 +0100
Subject: [PATCH] remove standalone test codes

---
 libsharp/sharp_acctest.c   | 267 ---------------------------
 libsharp/sharp_scaletest.c | 275 ----------------------------
 libsharp/sharp_test.c      | 272 ----------------------------
 libsharp/sharp_test_mpi.c  | 359 -------------------------------------
 4 files changed, 1173 deletions(-)
 delete mode 100644 libsharp/sharp_acctest.c
 delete mode 100644 libsharp/sharp_scaletest.c
 delete mode 100644 libsharp/sharp_test.c
 delete mode 100644 libsharp/sharp_test_mpi.c

diff --git a/libsharp/sharp_acctest.c b/libsharp/sharp_acctest.c
deleted file mode 100644
index 56ada7b..0000000
--- a/libsharp/sharp_acctest.c
+++ /dev/null
@@ -1,267 +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_acctest.c
-    Systematic accuracy test for libsharp.
-
-    Copyright (C) 2006-2012 Max-Planck-Society
-    \author Martin Reinecke
-*/
-
-#include <stdio.h>
-#include <string.h>
-#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 double drand (double min, double max)
-  { return min + (max-min)*rand()/(RAND_MAX+1.0); }
-
-static void random_alm (dcmplx *alm, sharp_alm_info *helper, int spin)
-  {
-  for (int mi=0;mi<helper->nm; ++mi)
-    {
-    int m=helper->mval[mi];
-    for (int l=m;l<=helper->lmax; ++l)
-      {
-      if ((l<spin)&&(m<spin))
-        alm[sharp_alm_index(helper,l,mi)] = 0.;
-      else
-        {
-        double rv = drand(-1,1);
-        double iv = (m==0) ? 0 : drand(-1,1);
-        alm[sharp_alm_index(helper,l,mi)] = rv+_Complex_I*iv;
-        }
-      }
-    }
-  }
-
-static void measure_errors (dcmplx **alm, dcmplx **alm2,
-  ptrdiff_t nalms, int ncomp)
-  {
-  for (int i=0; i<ncomp; ++i)
-    {
-    double sum=0, sum2=0, maxdiff=0;
-    for (ptrdiff_t m=0; m<nalms; ++m)
-      {
-      double x=creal(alm[i][m])-creal(alm2[i][m]),
-             y=cimag(alm[i][m])-cimag(alm2[i][m]);
-      sum+=x*x+y*y;
-      sum2+=creal(alm[i][m])*creal(alm[i][m])+cimag(alm[i][m])*cimag(alm[i][m]);
-      if (fabs(x)>maxdiff) maxdiff=fabs(x);
-      if (fabs(y)>maxdiff) maxdiff=fabs(y);
-      }
-    sum=sqrt(sum/nalms);
-    sum2=sqrt(sum2/nalms);
-    UTIL_ASSERT((maxdiff<1e-10)&&(sum/sum2<1e-10),"error");
-    }
-  }
-
-static void check_sign_scale(void)
-  {
-  int lmax=50;
-  int mmax=lmax;
-  sharp_geom_info *tinfo;
-  int nrings=lmax+1;
-  int ppring=2*lmax+2;
-  ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
-  sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-
-  /* flip theta to emulate the "old" Gaussian grid geometry */
-  for (int i=0; i<tinfo->npairs; ++i)
-    {
-    const double pi=3.141592653589793238462643383279502884197;
-    tinfo->pair[i].r1.cth=-tinfo->pair[i].r1.cth;
-    tinfo->pair[i].r2.cth=-tinfo->pair[i].r2.cth;
-    tinfo->pair[i].r1.theta=pi-tinfo->pair[i].r1.theta;
-    tinfo->pair[i].r2.theta=pi-tinfo->pair[i].r2.theta;
-    }
-
-  sharp_alm_info *alms;
-  sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
-  ptrdiff_t nalms = ((mmax+1)*(mmax+2))/2 + (mmax+1)*(lmax-mmax);
-
-  for (int ntrans=1; ntrans<10; ++ntrans)
-    {
-    double **map;
-    ALLOC2D(map,double,2*ntrans,npix);
-
-    dcmplx **alm;
-    ALLOC2D(alm,dcmplx,2*ntrans,nalms);
-    for (int i=0; i<2*ntrans; ++i)
-      for (int j=0; j<nalms; ++j)
-        alm[i][j]=1.+_Complex_I;
-
-    sharp_execute(SHARP_ALM2MAP,0,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,
-      NULL,NULL);
-    for (int it=0; it<ntrans; ++it)
-      {
-      UTIL_ASSERT(FAPPROX(map[it][0     ], 3.588246976618616912e+00,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[it][npix/2], 4.042209792157496651e+01,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[it][npix-1],-1.234675107554816442e+01,1e-12),
-        "error");
-      }
-    sharp_execute(SHARP_ALM2MAP,1,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,
-      NULL,NULL);
-    for (int it=0; it<ntrans; ++it)
-      {
-      UTIL_ASSERT(FAPPROX(map[2*it  ][0     ], 2.750897760535633285e+00,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it  ][npix/2], 3.137704477368562905e+01,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it  ][npix-1],-8.405730859837063917e+01,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it+1][0     ],-2.398026536095463346e+00,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it+1][npix/2],-4.961140548331700728e+01,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it+1][npix-1],-1.412765834230440021e+01,1e-12),
-        "error");
-      }
-
-    sharp_execute(SHARP_ALM2MAP,2,&alm[0],&map[0],tinfo,alms,ntrans,SHARP_DP,
-      NULL,NULL);
-    for (int it=0; it<ntrans; ++it)
-      {
-      UTIL_ASSERT(FAPPROX(map[2*it  ][0     ],-1.398186224727334448e+00,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it  ][npix/2],-2.456676000884031197e+01,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it  ][npix-1],-1.516249174408820863e+02,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it+1][0     ],-3.173406200299964119e+00,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it+1][npix/2],-5.831327404513146462e+01,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it+1][npix-1],-1.863257892248353897e+01,1e-12),
-        "error");
-      }
-
-    sharp_execute(SHARP_ALM2MAP_DERIV1,1,&alm[0],&map[0],tinfo,alms,ntrans,
-      SHARP_DP,NULL,NULL);
-    for (int it=0; it<ntrans; ++it)
-      {
-      UTIL_ASSERT(FAPPROX(map[2*it  ][0     ],-6.859393905369091105e-01,1e-11),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it  ][npix/2],-2.103947835973212364e+02,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it  ][npix-1],-1.092463246472086439e+03,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it+1][0     ],-1.411433220713928165e+02,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it+1][npix/2],-1.146122859381925082e+03,1e-12),
-        "error");
-      UTIL_ASSERT(FAPPROX(map[2*it+1][npix-1], 7.821618677689795049e+02,1e-12),
-        "error");
-      }
-
-    DEALLOC2D(map);
-    DEALLOC2D(alm);
-    }
-
-  sharp_destroy_alm_info(alms);
-  sharp_destroy_geom_info(tinfo);
-  }
-
-static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
-  ptrdiff_t mmax, ptrdiff_t npix, int spin, int ntrans, int nv)
-  {
-  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);
-
-  sharp_alm_info *alms;
-  sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
-
-  srand(4);
-  dcmplx **alm;
-  ALLOC2D(alm,dcmplx,ncomp,nalms);
-  for (int i=0; i<ncomp; ++i)
-    random_alm(alm[i],alms,spin);
-
-  dcmplx **alm2;
-  ALLOC2D(alm2,dcmplx,ncomp,nalms);
-
-  sharp_execute(SHARP_ALM2MAP,spin,&alm[0],&map[0],tinfo,alms,ntrans,
-    SHARP_DP|nv,NULL,NULL);
-  sharp_execute(SHARP_MAP2ALM,spin,&alm2[0],&map[0],tinfo,alms,ntrans,
-    SHARP_DP|nv,NULL,NULL);
-  measure_errors(alm,alm2,nalms,ncomp);
-
-  DEALLOC2D(map);
-  DEALLOC2D(alm);
-  DEALLOC2D(alm2);
-
-  sharp_destroy_alm_info(alms);
-  }
-
-int main(void)
-  {
-#ifdef USE_MPI
-  MPI_Init(NULL,NULL);
-#endif
-  sharp_module_startup("sharp_acctest",1,1,"",1);
-
-  int lmax=127;
-
-  printf("Checking signs and scales.\n");
-  check_sign_scale();
-  printf("Passed.\n\n");
-
-  printf("Testing map analysis accuracy.\n");
-
-  sharp_geom_info *tinfo;
-  int nrings=lmax+1;
-  int ppring=2*lmax+2;
-  ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
-  sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-  for (int nv=1; nv<=6; ++nv)
-    for (int ntrans=1; ntrans<=6; ++ntrans)
-      {
-      check_accuracy(tinfo,lmax,lmax,npix,0,ntrans,nv);
-      check_accuracy(tinfo,lmax,lmax,npix,1,ntrans,nv);
-      check_accuracy(tinfo,lmax,lmax,npix,2,ntrans,nv);
-      check_accuracy(tinfo,lmax,lmax,npix,3,ntrans,nv);
-      check_accuracy(tinfo,lmax,lmax,npix,30,ntrans,nv);
-      }
-  sharp_destroy_geom_info(tinfo);
-  printf("Passed.\n\n");
-
-#ifdef USE_MPI
-  MPI_Finalize();
-#endif
-  return 0;
-  }
diff --git a/libsharp/sharp_scaletest.c b/libsharp/sharp_scaletest.c
deleted file mode 100644
index 0acd1f3..0000000
--- a/libsharp/sharp_scaletest.c
+++ /dev/null
@@ -1,275 +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_scaletest.c
-    Copyright (C) 2012 Max-Planck-Society
-    \author Martin Reinecke
-*/
-
-#include <stdio.h>
-
-#if (defined(_OPENMP) && defined(USE_MPI))
-
-#include <stdlib.h>
-#include <string.h>
-#include <omp.h>
-#include <mpi.h>
-#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; i<ainfo->nm; 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; i<ginfo->npairs; 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; i<ainfo->nm; ++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; i<ginfo->npairs; ++i)
-    {
-    res += ginfo->pair[i].r1.nph;
-    if (ginfo->pair[i].r2.nph>0) res += ginfo->pair[i].r2.nph;
-    }
-  return res;
-  }
-
-static double drand (double min, double max)
-  { return min + (max-min)*rand()/(RAND_MAX+1.0); }
-
-static void random_alm (dcmplx *alm, sharp_alm_info *helper, int spin)
-  {
-  static int cnt=0;
-  ++cnt;
-  for (int mi=0;mi<helper->nm; ++mi)
-    {
-    int m=helper->mval[mi];
-    srand(1234567*cnt+8912*m);
-    for (int l=m;l<=helper->lmax; ++l)
-      {
-      if ((l<spin)&&(m<spin))
-        alm[sharp_alm_index(helper,l,mi)] = 0.;
-      else
-        {
-        double rv = drand(-1,1);
-        double iv = (m==0) ? 0 : drand(-1,1);
-        alm[sharp_alm_index(helper,l,mi)] = rv+_Complex_I*iv;
-        }
-      }
-    }
-  }
-
-static void measure_errors (dcmplx **alm, double *sqsum,
-  const sharp_alm_info *ainfo, int ncomp)
-  {
-  long nalms=get_nalms(ainfo), nalms_tot;
-  MPI_Allreduce(&nalms,&nalms_tot,1,MPI_LONG,MPI_SUM,MPI_COMM_WORLD);
-
-  for (int i=0; i<ncomp; ++i)
-    {
-    double sum=0, maxdiff=0, sumtot, sqsumtot, maxdifftot;
-    for (int mi=0; mi<ainfo->nm; ++mi)
-      {
-      int m=ainfo->mval[mi];
-      for (int l=m; l<=ainfo->lmax; ++l)
-        {
-        ptrdiff_t idx=sharp_alm_index(ainfo,l,mi);
-        sum+=creal(alm[i][idx])*creal(alm[i][idx])
-            +cimag(alm[i][idx])*cimag(alm[i][idx]);
-        if (fabs(creal(alm[i][idx]))>maxdiff) maxdiff=fabs(creal(alm[i][idx]));
-        if (fabs(cimag(alm[i][idx]))>maxdiff) maxdiff=fabs(cimag(alm[i][idx]));
-        }
-      }
-
-    MPI_Allreduce(&sum,&sumtot,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
-    MPI_Allreduce(&sqsum[i],&sqsumtot,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
-    MPI_Allreduce(&maxdiff,&maxdifftot,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
-    sumtot=sqrt(sumtot/nalms_tot);
-    sqsumtot=sqrt(sqsumtot/nalms_tot);
-    if (mytask==0)
-      printf("component %i: rms %e, maxerr %e\n",i,sumtot/sqsumtot, maxdifftot);
-    }
-  }
-
-int main(int argc, char **argv)
-  {
-  MPI_Init(NULL,NULL);
-  MPI_Comm_size(MPI_COMM_WORLD,&ntasks);
-  MPI_Comm_rank(MPI_COMM_WORLD,&mytask);
-
-  sharp_module_startup("sharp_scaletest",argc,4,"<lmax> <nphi> <spin>",
-    mytask==0);
-
-  int lmax=atoi(argv[1]);
-  int spin=atoi(argv[3]);
-
-  sharp_geom_info *tinfo;
-  ptrdiff_t npix=0;
-  int nrings=lmax+1;
-  int ppring=atoi(argv[2]);
-  sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-
-  reduce_geom_info(tinfo);
-  npix=get_npix(tinfo);
-
-  int mmax=lmax;
-  int ncomp = (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);
-
-  if (mytask==0) printf("Testing map analysis accuracy.\n");
-  if (mytask==0) printf("lmax=%d, spin=%d\n", lmax, spin);
-  if (mytask==0)
-    printf("\nTesting Gaussian grid (%d rings, %d pixels/ring, %ld pixels)\n",
-      nrings,ppring,(long)npix);
-
-  for (int n=0; n<ncomp; ++n)
-    random_alm(alm[n],alms,spin);
-
-  double time=0;
-  unsigned long long opcnt=0;
-  sharp_execute_mpi(MPI_COMM_WORLD,SHARP_ALM2MAP,spin,&alm[0],&map[0],
-    tinfo,alms,1,SHARP_DP,&time,&opcnt);
-
-  time=maxTime(time);
-  opcnt=totalops(opcnt);
-  if (mytask==0) printf("wall time for alm2map: %fs\n",time);
-  if (mytask==0) printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);
-
-  double *sqsum=RALLOC(double,ncomp);
-  for (int i=0; i<ncomp; ++i)
-    {
-    sqsum[i]=0;
-    for (ptrdiff_t j=0; j<nalms; ++j)
-      {
-      sqsum[i]+=creal(alm[i][j])*creal(alm[i][j])
-               +cimag(alm[i][j])*cimag(alm[i][j]);
-      alm[i][j]=-alm[i][j];
-      }
-    }
-
-  sharp_execute_mpi(MPI_COMM_WORLD,SHARP_MAP2ALM,spin,&alm[0],&map[0],
-    tinfo,alms,1,SHARP_DP|SHARP_ADD,&time,&opcnt);
-
-  time=maxTime(time);
-  opcnt=totalops(opcnt);
-  if (mytask==0) printf("wall time for map2alm: %fs\n",time);
-  if (mytask==0) printf("Performance: %fGFLOPs/s\n",1e-9*opcnt/time);
-
-  measure_errors(alm,sqsum,alms,ncomp);
-
-  DEALLOC(sqsum);
-  DEALLOC2D(map);
-  DEALLOC2D(alm);
-
-  sharp_destroy_alm_info(alms);
-  sharp_destroy_geom_info(tinfo);
-
-  double mHWM=totalMem(VmHWM());
-  if (mytask==0) printf("\nMemory high water mark: %.2f MB\n",mHWM/(1<<20));
-
-  MPI_Finalize();
-  return 0;
-  }
-
-#else
-
-#include "c_utils.h"
-
-int main(void)
-  { UTIL_FAIL("Need OpenMP and MPI"); return 1; }
-
-#endif
diff --git a/libsharp/sharp_test.c b/libsharp/sharp_test.c
deleted file mode 100644
index dc9bbff..0000000
--- a/libsharp/sharp_test.c
+++ /dev/null
@@ -1,272 +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_test.c
-    Accuracy test for libsharp's map analysis.
-
-    This program first generates a_lm coefficients up to
-    a user-specified lmax (with mmax=lmax); where applicable, the
-    real and imaginary parts of the coefficients are uniform
-    random numbers of the interval [-1;1[.
-    Afterwards, the random a_lm are converted to a map.
-    This map is analyzed (optionally using an iterative scheme
-    with a user-supplied number of steps).
-    After every iteration, the code then outputs the RMS of the residual a_lm
-    (i.e. the difference between the current and original a_lm), divided by
-    the RMS of the original a_lm, as well as the maximum absolute change of any
-    real or imaginary part between the current and original a_lm.
-
-    This operation can be performed for several different pixelisations:
-      - a Gaussian with the minimal number of rings for exact analysis
-        and a user-defined ring resolution
-      - an ECP grid with the minimal number of rings for exact analysis
-        and a user-defined ring resolution
-      - a Healpix grid with a user-defined Nside parameter.
-
-    The user can specify the spin of the desired transform.
-
-    Copyright (C) 2006-2012 Max-Planck-Society
-    \author Martin Reinecke
-*/
-
-#include <stdio.h>
-#include <string.h>
-#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"
-#include "memusage.h"
-
-typedef complex double dcmplx;
-
-static double drand (double min, double max)
-  { return min + (max-min)*rand()/(RAND_MAX+1.0); }
-
-static void random_alm (dcmplx *alm, sharp_alm_info *helper, int spin)
-  {
-  for (int mi=0;mi<helper->nm; ++mi)
-    {
-    int m=helper->mval[mi];
-    for (int l=m;l<=helper->lmax; ++l)
-      {
-      if ((l<spin)&&(m<spin))
-        alm[sharp_alm_index(helper,l,mi)] = 0.;
-      else
-        {
-        double rv = drand(-1,1);
-        double iv = (m==0) ? 0 : drand(-1,1);
-        alm[sharp_alm_index(helper,l,mi)] = rv+_Complex_I*iv;
-        }
-      }
-    }
-  }
-
-static void measure_errors (dcmplx **alm, dcmplx **alm2,
-  ptrdiff_t nalms, int ncomp)
-  {
-  for (int i=0; i<ncomp; ++i)
-    {
-    double sum=0, sum2=0, maxdiff=0;
-    for (ptrdiff_t m=0; m<nalms; ++m)
-      {
-      double x=creal(alm[i][m])-creal(alm2[i][m]),
-             y=cimag(alm[i][m])-cimag(alm2[i][m]);
-      sum+=x*x+y*y;
-      sum2+=creal(alm[i][m])*creal(alm[i][m])+cimag(alm[i][m])*cimag(alm[i][m]);
-      if (fabs(x)>maxdiff) maxdiff=fabs(x);
-      if (fabs(y)>maxdiff) maxdiff=fabs(y);
-      }
-    sum=sqrt(sum/nalms);
-    sum2=sqrt(sum2/nalms);
-    printf("component %i: rms %e, maxerr %e\n",i, sum/sum2, maxdiff);
-    }
-  }
-
-static void map2alm_iter (sharp_geom_info *tinfo, double **map,
-  dcmplx **alm_orig, dcmplx **alm, int lmax, int mmax,
-  ptrdiff_t npix, ptrdiff_t nalms, int spin, int ntrans, int niter)
-  {
-  int ncomp = ntrans*((spin==0) ? 1 : 2);
-
-  sharp_alm_info *alms;
-  sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
-
-  double time;
-  unsigned long long opcnt;
-  sharp_execute(SHARP_MAP2ALM,spin,&alm[0],&map[0],tinfo,alms,ntrans,
-    SHARP_DP,&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)
-    {
-    double **map2;
-    ALLOC2D(map2,double,ncomp,npix);
-    printf ("\niteration %i:\n", iter+1);
-    sharp_execute(SHARP_ALM2MAP,spin,&alm[0],&map2[0],tinfo,alms,ntrans,
-      SHARP_DP,&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];
-
-    sharp_execute(SHARP_MAP2ALM,spin,&alm[0],&map2[0],tinfo,alms,ntrans,
-      SHARP_DP|SHARP_ADD,&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);
-    }
-
-  sharp_destroy_alm_info(alms);
-  }
-
-static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
-  ptrdiff_t mmax, ptrdiff_t npix, int spin, int ntrans, int niter)
-  {
-  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);
-
-  sharp_alm_info *alms;
-  sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
-
-  srand(4);
-  dcmplx **alm;
-  ALLOC2D(alm,dcmplx,ncomp,nalms);
-  for (int i=0; i<ncomp; ++i)
-    random_alm(alm[i],alms,spin);
-
-  dcmplx **alm2;
-  ALLOC2D(alm2,dcmplx,ncomp,nalms);
-
-  double time;
-  unsigned long long opcnt;
-  printf ("\niteration 0:\n");
-  sharp_execute(SHARP_ALM2MAP,spin,&alm[0],&map[0],tinfo,alms,ntrans,
-    SHARP_DP,&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);
-
-  DEALLOC2D(map);
-  DEALLOC2D(alm);
-  DEALLOC2D(alm2);
-
-  sharp_destroy_alm_info(alms);
-  }
-
-int main(int argc, char **argv)
-  {
-#ifdef USE_MPI
-  MPI_Init(NULL,NULL);
-#endif
-  sharp_module_startup("sharp_test",argc,7,
-    "<healpix|ecp|gauss|hw|fejer2> <lmax> <nside|nphi> <niter> <spin> <ntrans>",
-    1);
-
-  int lmax=atoi(argv[2]);
-  int niter=atoi(argv[4]);
-  int spin=atoi(argv[5]);
-  int ntrans=atoi(argv[6]);
-
-  printf("Testing map analysis accuracy.\n");
-  printf("lmax=%d, %d iterations, spin=%d\n", lmax, niter, spin);
-
-  sharp_geom_info *tinfo;
-  if (strcmp(argv[1],"gauss")==0)
-    {
-    int nrings=lmax+1;
-    int ppring=atoi(argv[3]);
-    ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
-    printf("\nTesting Gaussian grid (%d rings, %d pixels/ring, %ld pixels)\n",
-          nrings,ppring,(long)npix);
-    sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-    check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
-    sharp_destroy_geom_info(tinfo);
-    }
-  else if (strcmp(argv[1],"ecp")==0)
-    {
-    int nrings=2*lmax+1;
-    int ppring=atoi(argv[3]);
-    ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
-    printf("\nTesting ECP grid (%d rings, %d pixels/ring, %ld pixels)\n",
-          nrings,ppring,(long)npix);
-    sharp_make_ecp_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-    check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
-    sharp_destroy_geom_info(tinfo);
-    }
-  else if (strcmp(argv[1],"healpix")==0)
-    {
-    int nside=atoi(argv[3]);
-    if (nside<1) nside=1;
-    ptrdiff_t npix=12*(ptrdiff_t)nside*nside;
-    printf("\nTesting Healpix grid (nside=%d, %ld pixels)\n",
-          nside,(long)npix);
-    sharp_make_healpix_geom_info (nside, 1, &tinfo);
-    check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
-    sharp_destroy_geom_info(tinfo);
-    }
-  else if (strcmp(argv[1],"hw")==0)
-    {
-    int nrings=2*lmax+1;
-    int ppring=atoi(argv[3]);
-    ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
-    printf("\nTesting HW grid (%d rings, %d pixels/ring, %ld pixels)\n",
-          nrings,ppring,(long)npix);
-    sharp_make_hw_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-    check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
-    sharp_destroy_geom_info(tinfo);
-    }
-  else if (strcmp(argv[1],"fejer2")==0)
-    {
-    int nrings=2*lmax+1;
-    int ppring=atoi(argv[3]);
-    ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
-    printf("\nTesting Fejer2 grid (%d rings, %d pixels/ring, %ld pixels)\n",
-          nrings,ppring,(long)npix);
-    sharp_make_fejer2_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-    check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
-    sharp_destroy_geom_info(tinfo);
-    }
-  else
-    UTIL_FAIL("unknown grid geometry");
-
-  printf("\nMemory high water mark: %.2f MB\n",VmHWM()/(1<<20));
-
-#ifdef USE_MPI
-  MPI_Finalize();
-#endif
-  return 0;
-  }
diff --git a/libsharp/sharp_test_mpi.c b/libsharp/sharp_test_mpi.c
deleted file mode 100644
index 2c87fde..0000000
--- a/libsharp/sharp_test_mpi.c
+++ /dev/null
@@ -1,359 +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_test_mpi.c
-    Accuracy test for libsharp's map analysis with MPI support.
-
-    This program first generates a_lm coefficients up to
-    a user-specified lmax (with mmax=lmax); where applicable, the
-    real and imaginary parts of the coefficients are uniform
-    random numbers of the interval [-1;1[.
-    Afterwards, the random a_lm are converted to a map.
-    This map is analyzed (optionally using an iterative scheme
-    with a user-supplied number of steps).
-    After every iteration, the code then outputs the RMS of the residual a_lm
-    (i.e. the difference between the current and original a_lm), divided by
-    the RMS of the original a_lm, as well as the maximum absolute change of any
-    real or imaginary part between the current and original a_lm.
-
-    This operation can be performed for several different pixelisations:
-      - a Gaussian with the minimal number of rings for exact analysis
-        and a user-defined ring resolution
-      - an ECP grid with the minimal number of rings for exact analysis
-        and a user-defined ring resolution
-      - a Healpix grid with a user-defined Nside parameter.
-
-    The user can specify the spin of the desired transform.
-
-    Copyright (C) 2006-2012 Max-Planck-Society
-    \author Martin Reinecke
-*/
-
-#ifdef USE_MPI
-
-#include <stdio.h>
-#include <string.h>
-#include "sharp_mpi.h"
-#include "sharp_geomhelpers.h"
-#include "sharp_almhelpers.h"
-#include "c_utils.h"
-#include "walltime_c.h"
-#include "sharp_announce.h"
-#include "sharp_core.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 drand (double min, double max)
-  { return min + (max-min)*rand()/(RAND_MAX+1.0); }
-
-static ptrdiff_t get_nalms(const sharp_alm_info *ainfo)
-  {
-  ptrdiff_t res=0;
-  for (int i=0; i<ainfo->nm; ++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; i<ginfo->npairs; ++i)
-    {
-    res += ginfo->pair[i].r1.nph;
-    if (ginfo->pair[i].r2.nph>0) res += ginfo->pair[i].r2.nph;
-    }
-  return res;
-  }
-
-static void reduce_alm_info(sharp_alm_info *ainfo)
-  {
-  int nmnew=0;
-  ptrdiff_t ofs = 0;
-  for (int i=mytask; i<ainfo->nm; 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; i<ginfo->npairs; 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 void random_alm (dcmplx *alm, sharp_alm_info *helper, int spin)
-  {
-  static int cnt=0;
-  ++cnt;
-  for (int mi=0;mi<helper->nm; ++mi)
-    {
-    int m=helper->mval[mi];
-    srand(1234567*cnt+8912*m);
-    for (int l=m;l<=helper->lmax; ++l)
-      {
-      if ((l<spin)&&(m<spin))
-        alm[sharp_alm_index(helper,l,mi)] = 0.;
-      else
-        {
-        double rv = drand(-1,1);
-        double iv = (m==0) ? 0 : drand(-1,1);
-        alm[sharp_alm_index(helper,l,mi)] = rv+_Complex_I*iv;
-        }
-      }
-    }
-  }
-
-static void measure_errors (dcmplx **alm, dcmplx **alm2,
-  const sharp_alm_info *ainfo, int ncomp)
-  {
-  long nalms=get_nalms(ainfo), nalms_tot;
-  MPI_Allreduce(&nalms,&nalms_tot,1,MPI_LONG,MPI_SUM,MPI_COMM_WORLD);
-
-  for (int i=0; i<ncomp; ++i)
-    {
-    double sum=0, sum2=0, maxdiff=0, sumtot, sum2tot, maxdifftot;
-    for (int mi=0; mi<ainfo->nm; ++mi)
-      {
-      int m=ainfo->mval[mi];
-      for (int l=m; l<=ainfo->lmax; ++l)
-        {
-        ptrdiff_t idx=sharp_alm_index(ainfo,l,mi);
-        double x=creal(alm[i][idx])-creal(alm2[i][idx]),
-               y=cimag(alm[i][idx])-cimag(alm2[i][idx]);
-        sum+=x*x+y*y;
-        sum2+=creal(alm[i][idx])*creal(alm[i][idx])
-             +cimag(alm[i][idx])*cimag(alm[i][idx]);
-        if (fabs(x)>maxdiff) maxdiff=fabs(x);
-        if (fabs(y)>maxdiff) maxdiff=fabs(y);
-        }
-      }
-
-    MPI_Allreduce(&sum,&sumtot,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
-    MPI_Allreduce(&sum2,&sum2tot,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
-    MPI_Allreduce(&maxdiff,&maxdifftot,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
-    sumtot=sqrt(sumtot/nalms_tot);
-    sum2tot=sqrt(sum2tot/nalms_tot);
-    if (mytask==0)
-      printf("component %i: rms %e, maxerr %e\n",i, sumtot/sum2tot, maxdifftot);
-    }
-  }
-
-static void map2alm_iter (sharp_geom_info *tinfo, double **map,
-  dcmplx **alm_orig, dcmplx **alm, int lmax, int mmax,
-  ptrdiff_t npix, int spin, int ntrans, int niter)
-  {
-  int ncomp = ntrans*((spin==0) ? 1 : 2);
-
-  sharp_alm_info *alms;
-  sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
-  reduce_alm_info(alms);
-
-  double jtime;
-  unsigned long long jopcnt;
-
-  sharp_execute_mpi(MPI_COMM_WORLD,SHARP_MAP2ALM,spin,&alm[0],&map[0],
-    tinfo,alms,ntrans,SHARP_DP,&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);
-
-  for (int iter=0; iter<niter; ++iter)
-    {
-    double **map2;
-    ALLOC2D(map2,double,ncomp,npix);
-    if (mytask==0) printf ("\niteration %i:\n", iter+1);
-    sharp_execute_mpi(MPI_COMM_WORLD,SHARP_ALM2MAP,spin,&alm[0],&map2[0],
-      tinfo,alms,ntrans,SHARP_DP,&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];
-
-    sharp_execute_mpi(MPI_COMM_WORLD,SHARP_MAP2ALM,spin,&alm[0],&map2[0],
-      tinfo,alms,ntrans,SHARP_DP|SHARP_ADD,&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);
-    measure_errors(alm_orig,alm,alms,ncomp);
-    }
-
-  sharp_destroy_alm_info(alms);
-  }
-
-static void check_accuracy (sharp_geom_info *tinfo, ptrdiff_t lmax,
-  ptrdiff_t mmax, ptrdiff_t npix, int spin, int ntrans, int niter)
-  {
-  int ncomp = ntrans*((spin==0) ? 1 : 2);
-
-  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);
-  reduce_alm_info(alms);
-  nalms=get_nalms(alms);
-
-  dcmplx **alm;
-  ALLOC2D(alm,dcmplx,ncomp,nalms);
-  srand(4);
-  for (int i=0; i<ncomp; ++i)
-    random_alm(alm[i],alms,spin);
-
-  dcmplx **alm2;
-  ALLOC2D(alm2,dcmplx,ncomp,nalms);
-
-  if (mytask==0) printf ("\niteration 0:\n");
-  sharp_execute_mpi(MPI_COMM_WORLD,SHARP_ALM2MAP,spin,&alm[0],&map[0],
-    tinfo,alms,ntrans,SHARP_DP,&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);
-
-  map2alm_iter(tinfo, map, alm, alm2, lmax, mmax, npix, spin, ntrans, niter);
-
-  DEALLOC2D(map);
-  DEALLOC2D(alm);
-  DEALLOC2D(alm2);
-
-  sharp_destroy_alm_info(alms);
-  }
-
-int main(int argc, char **argv)
-  {
-  MPI_Init(NULL,NULL);
-  MPI_Comm_size(MPI_COMM_WORLD,&ntasks);
-  MPI_Comm_rank(MPI_COMM_WORLD,&mytask);
-
-  sharp_module_startup("sharp_test_mpi",argc,7,
-    "<healpix|ecp|gauss> <lmax> <nside|nphi> <niter> <spin> <ntrans>",
-    mytask==0);
-  int lmax=atoi(argv[2]);
-  int niter=atoi(argv[4]);
-  int spin=atoi(argv[5]);
-  int ntrans=atoi(argv[6]);
-
-  if (mytask==0)
-    {
-    printf("Testing map analysis accuracy.\n");
-    printf("lmax=%d, %d iterations, spin=%d\n", lmax, niter, spin);
-    }
-
-  sharp_geom_info *tinfo;
-  if (strcmp(argv[1],"gauss")==0)
-    {
-    int nrings=lmax+1;
-    int ppring=atoi(argv[3]);
-    ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
-    if (mytask==0)
-      printf("\nTesting Gaussian grid (%d rings, %d pixels/ring, %ld pixels)\n",
-             nrings,ppring,(long)npix);
-    sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-    reduce_geom_info(tinfo);
-    npix=get_npix(tinfo);
-    check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
-    sharp_destroy_geom_info(tinfo);
-    }
-  else if (strcmp(argv[1],"ecp")==0)
-    {
-    int nrings=2*lmax+2;
-    int ppring=atoi(argv[3]);
-    ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
-    if (mytask==0)
-      printf("\nTesting ECP grid (%d rings, %d pixels/ring, %ld pixels)\n",
-             nrings,ppring,(long)npix);
-    sharp_make_ecp_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
-    reduce_geom_info(tinfo);
-    npix=get_npix(tinfo);
-    check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
-    sharp_destroy_geom_info(tinfo);
-    }
-  else if (strcmp(argv[1],"healpix")==0)
-    {
-    int nside=atoi(argv[3]);
-    if (nside<1) nside=1;
-    ptrdiff_t npix=12*(ptrdiff_t)nside*nside;
-    if (mytask==0)
-      printf("\nTesting Healpix grid (nside=%d, %ld pixels)\n",
-             nside,(long)npix);
-    sharp_make_healpix_geom_info (nside, 1, &tinfo);
-    reduce_geom_info(tinfo);
-    npix=get_npix(tinfo);
-    check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
-    sharp_destroy_geom_info(tinfo);
-    }
-  else
-    UTIL_FAIL("unknown grid geometry");
-
-  MPI_Finalize();
-  return 0;
-  }
-
-#else
-
-#include "c_utils.h"
-
-int main(void)
-  { UTIL_FAIL("MPI support not enabled."); return 1; }
-
-#endif