more c++ification

Makefile.am

@@ -3,8 +3,6 @@ ACLOCAL_AMFLAGS = -I m4
 lib_LTLIBRARIES = libmrutil.la
 
 libmrutil_la_SOURCES = \
-  libsharp2/sharp_utils.h \
-  libsharp2/sharp_utils.cc \
   libsharp2/sharp.cc \
   libsharp2/sharp_almhelpers.cc \
   libsharp2/sharp_core.cc \
@@ -12,6 +10,7 @@ libmrutil_la_SOURCES = \
   libsharp2/sharp_ylmgen.cc \
   libsharp2/sharp_internal.h \
   libsharp2/sharp_ylmgen.h \
+  mr_util/math_utils.h \
   mr_util/cmplx.h \
   mr_util/error_handling.cc \
   mr_util/error_handling.h \

libsharp2/sharp.cc

@@ -26,13 +26,13 @@
  */
 
 #include <cmath>
-#include <cstring>
 #include <atomic>
 #include <memory>
+#include "mr_util/math_utils.h"
 #include "mr_util/fft.h"
+#include "mr_util/aligned_array.h"
 #include "libsharp2/sharp_ylmgen.h"
 #include "libsharp2/sharp_internal.h"
-#include "libsharp2/sharp_utils.h"
 #include "libsharp2/sharp_almhelpers.h"
 #include "libsharp2/sharp_geomhelpers.h"
 #include "mr_util/threading.h"
@@ -41,6 +41,7 @@
 #include "mr_util/timers.h"
 
 using namespace std;
+using namespace mr;
 
 using dcmplx = complex<double>;
 using fcmplx = complex<float>;
@@ -92,7 +93,7 @@ struct ringhelper
     {
     norot = (fabs(phi0)<1e-14);
     if (!(norot))
-      if ((mmax!=s_shift-1) || (!FAPPROX(phi0,phi0_,1e-12)))
+      if ((mmax!=s_shift-1) || (!approx(phi0,phi0_,1e-12)))
       {
       shiftarr.resize(mmax+1);
       s_shift = mmax+1;
@@ -111,22 +112,6 @@ struct ringhelper
     }
   };
 
-static int ringinfo_compare (const void *xa, const void *xb)
-  {
-  const sharp_ringinfo *a=(const sharp_ringinfo *)xa, *b=(const sharp_ringinfo *)xb;
-  return (a->sth < b->sth) ? -1 : (a->sth > b->sth) ? 1 : 0;
-  }
-static int ringpair_compare (const void *xa, const void *xb)
-  {
-  const sharp_ringpair *a=(const sharp_ringpair *)xa, *b=(const sharp_ringpair *)xb;
-//  return (a->r1.sth < b->r1.sth) ? -1 : (a->r1.sth > b->r1.sth) ? 1 : 0;
-  if (a->r1.nph==b->r1.nph)
-    return (a->r1.phi0 < b->r1.phi0) ? -1 :
-      ((a->r1.phi0 > b->r1.phi0) ? 1 :
-        (a->r1.cth>b->r1.cth ? -1 : 1));
-  return (a->r1.nph<b->r1.nph) ? -1 : 1;
-  }
-
 void sharp_make_general_alm_info (int lmax, int nm, int stride, const int *mval,
   const ptrdiff_t *mstart, int flags, sharp_alm_info **alm_info)
   {
@@ -204,11 +189,12 @@ void sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
     infos[m].nph = nph[m];
     if (info->nphmax<nph[m]) info->nphmax=nph[m];
     }
-  qsort(infos.data(),nrings,sizeof(sharp_ringinfo),ringinfo_compare);
+  sort(infos.begin(), infos.end(),[](const sharp_ringinfo &a, const sharp_ringinfo &b)
+    { return (a.sth<b.sth); });
   while (pos<nrings)
     {
     info->pair[npairs].r1=infos[pos];
-    if ((pos<nrings-1) && FAPPROX(infos[pos].cth,-infos[pos+1].cth,1e-12))
+    if ((pos<nrings-1) && approx(infos[pos].cth,-infos[pos+1].cth,1e-12))
       {
       if (infos[pos].cth>0)  // make sure northern ring is in r1
         info->pair[npairs].r2=infos[pos+1];
@@ -224,9 +210,16 @@ void sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
     ++pos;
     ++npairs;
     }
-
-  qsort(info->pair.data(),npairs,sizeof(sharp_ringpair),ringpair_compare);
   info->pair.resize(npairs);
+
+  sort(info->pair.begin(), info->pair.end(), [] (const sharp_ringpair &a, const sharp_ringpair &b)
+    {
+    if (a.r1.nph==b.r1.nph)
+    return (a.r1.phi0 < b.r1.phi0) ? true :
+      ((a.r1.phi0 > b.r1.phi0) ? false :
+        (a.r1.cth>b.r1.cth));
+    return a.r1.nph<b.r1.nph;
+    });
   }
 
 ptrdiff_t sharp_map_size(const sharp_geom_info *info)
@@ -476,7 +469,7 @@ MRUTIL_NOINLINE static void init_output (sharp_job *job)
       clear_map (job->ginfo,job->map[i],job->flags);
   }
 
-MRUTIL_NOINLINE static void alloc_phase (sharp_job *job, int nm, int ntheta)
+MRUTIL_NOINLINE static void alloc_phase (sharp_job *job, int nm, int ntheta, aligned_array<dcmplx> &data)
   {
   if (job->type==SHARP_MAP2ALM)
     {
@@ -490,17 +483,15 @@ MRUTIL_NOINLINE static void alloc_phase (sharp_job *job, int nm, int ntheta)
     if (((job->s_th*16*ntheta)&1023)==0) ntheta+=3; // hack to avoid critical strides
     job->s_m=job->s_th*ntheta;
     }
-  job->phase=RALLOC(dcmplx,2*job->nmaps*nm*ntheta);
+  data.resize(2*job->nmaps*nm*ntheta);
+  job->phase=data.data();
   }
 
-static void dealloc_phase (sharp_job *job)
-  { DEALLOC(job->phase); }
-
-static void alloc_almtmp (sharp_job *job, int lmax)
-  { job->almtmp=RALLOC(dcmplx,job->nalm*(lmax+2)); }
-
-static void dealloc_almtmp (sharp_job *job)
-  { DEALLOC(job->almtmp); }
+static void alloc_almtmp (sharp_job *job, int lmax, aligned_array<dcmplx> &data)
+  {
+  data.resize(job->nalm*(lmax+2));
+  job->almtmp=data.data();
+  }
 
 MRUTIL_NOINLINE static void alm2almtmp (sharp_job *job, int lmax, int mi)
   {
@@ -832,8 +823,9 @@ MRUTIL_NOINLINE static void sharp_execute_job (sharp_job *job)
   int nchunks, chunksize;
   get_chunk_info(job->ginfo->pair.size(),sharp_veclen()*sharp_max_nvec(job->spin),
                  &nchunks,&chunksize);
+  aligned_array<dcmplx> phasebuffer;
 //FIXME: needs to be changed to "nm"
-  alloc_phase (job,mmax+1,chunksize);
+  alloc_phase (job,mmax+1,chunksize, phasebuffer);
   std::atomic<size_t> opcnt = 0;
 
 /* chunk loop */
@@ -859,7 +851,8 @@ MRUTIL_NOINLINE static void sharp_execute_job (sharp_job *job)
       sharp_job ljob = *job;
       ljob.opcnt=0;
       sharp_Ylmgen generator(lmax,mmax,ljob.spin);
-      alloc_almtmp(&ljob,lmax);
+      aligned_array<dcmplx> almbuffer;
+      alloc_almtmp(&ljob,lmax,almbuffer);
 
       while (auto rng=sched.getNext()) for(auto mi=rng.lo; mi<rng.hi; ++mi)
         {
@@ -872,8 +865,6 @@ MRUTIL_NOINLINE static void sharp_execute_job (sharp_job *job)
         almtmp2alm (&ljob, lmax, mi);
         }
 
-      dealloc_almtmp(&ljob);
-
       opcnt+=ljob.opcnt;
       }); /* end of parallel region */
 
@@ -881,7 +872,6 @@ MRUTIL_NOINLINE static void sharp_execute_job (sharp_job *job)
     phase2map (job, mmax, llim, ulim);
     } /* end of chunk loop */
 
-  dealloc_phase (job);
   job->opcnt = opcnt;
   job->time=timer();
   }

libsharp2/sharp_almhelpers.cc

@@ -26,7 +26,6 @@
  */
 
 #include "libsharp2/sharp_almhelpers.h"
-#include "libsharp2/sharp_utils.h"
 
 void sharp_make_triangular_alm_info (int lmax, int mmax, int stride,
   sharp_alm_info **alm_info)

libsharp2/sharp_core_inc.cc

@@ -35,11 +35,10 @@
 #define XARCH(a) XCONCATX2(a,ARCH)
 
 #include <complex>
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
 #include "libsharp2/sharp.h"
 #include "libsharp2/sharp_internal.h"
-#include "libsharp2/sharp_utils.h"
 #include "mr_util/error_handling.h"
 #include "mr_util/useful_macros.h"
 #include "mr_util/simd.h"
@@ -84,8 +83,8 @@ static inline void vhsum_cmplx_special (Tv a, Tv b, Tv c, Tv d,
 
 using dcmplx = complex<double>;
 
-#define nv0 (128/VLEN)
-#define nvx (64/VLEN)
+constexpr size_t nv0 = 128/VLEN;
+constexpr size_t nvx = 64/VLEN;
 
 using Tbv0 = Tv[nv0];
 using Tbs0 = double[nv0*VLEN];

libsharp2/sharp_utils.cc

-/*
- *  This file is part of libsharp2.
- *
- *  libsharp2 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.
- *
- *  libsharp2 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 libsharp2; if not, write to the Free Software
- *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
- */
-
-/* libsharp2 is being developed at the Max-Planck-Institut fuer Astrophysik */
-
-/*
- *  Convenience functions
- *
- *  Copyright (C) 2008-2019 Max-Planck-Society
- *  Author: Martin Reinecke
- */
-
-#include <stdio.h>
-#include "libsharp2/sharp_utils.h"
-#include "mr_util/error_handling.h"
-
-#pragma GCC visibility push(hidden)
-
-/* This function tries to avoid allocations with a total size close to a high
-   power of two (called the "critical stride" here), by adding a few more bytes
-   if necessary. This lowers the probability that two arrays differ by a multiple
-   of the critical stride in their starting address, which in turn lowers the
-   risk of cache line contention. */
-static size_t manipsize(size_t sz)
-  {
-  const size_t critical_stride=4096, cacheline=64, overhead=32;
-  if (sz < (critical_stride/2)) return sz;
-  if (((sz+overhead)%critical_stride)>(2*cacheline)) return sz;
-  return sz+2*cacheline;
-  }
-#pragma GCC visibility pop
-
-#ifdef __SSE__
-#include <xmmintrin.h>
-#pragma GCC visibility push(hidden)
-void *sharp_malloc_ (size_t sz)
-  {
-  void *res;
-  if (sz==0) return NULL;
-  res = _mm_malloc(manipsize(sz),32);
-  MR_assert(res,"_mm_malloc() failed");
-  return res;
-  }
-void sharp_free_ (void *ptr)
-  { if ((ptr)!=NULL) _mm_free(ptr); }
-#pragma GCC visibility pop
-#else
-#pragma GCC visibility push(hidden)
-void *sharp_malloc_ (size_t sz)
-  {
-  void *res;
-  if (sz==0) return NULL;
-  res = malloc(manipsize(sz));
-  MR_assert(res,"malloc() failed");
-  return res;
-  }
-void sharp_free_ (void *ptr)
-  { if ((ptr)!=NULL) free(ptr); }
-#pragma GCC visibility pop
-#endif

libsharp2/sharp_utils.h

-/*
- *  This file is part of libc_utils.
- *
- *  libc_utils 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.
- *
- *  libc_utils 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 libc_utils; if not, write to the Free Software
- *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
- */
-
-/* libc_utils is being developed at the Max-Planck-Institut fuer Astrophysik */
-
-/*! \file c_utils.h
- *  Convenience functions
- *
- *  Copyright (C) 2008-2019 Max-Planck-Society
- *  \author Martin Reinecke
- *  \note This file should only be included from .c files, NOT from .h files.
- */
-
-#ifndef SHARP_UTILS_H
-#define SHARP_UTILS_H
-
-#include <math.h>
-#include <stdlib.h>
-#include <stddef.h>
-
-void *sharp_malloc_ (size_t sz);
-void sharp_free_ (void *ptr);
-
-/*! \def ALLOC(ptr,type,num)
-    Allocate space for \a num objects of type \a type. Make sure that the
-    allocation succeeded, else stop the program with an error. Return the
-    resulting pointer in \a ptr. */
-#define ALLOC(ptr,type,num) \
-  do { (ptr)=(type *)sharp_malloc_((num)*sizeof(type)); } while (0)
-/*! \def RALLOC(type,num)
-    Allocate space for \a num objects of type \a type. Make sure that the
-    allocation succeeded, else stop the program with an error. Cast the
-    resulting pointer to \a (type*). */
-#define RALLOC(type,num) \
-  ((type *)sharp_malloc_((num)*sizeof(type)))
-/*! \def DEALLOC(ptr)
-    Deallocate \a ptr. It must have been allocated using \a ALLOC or
-    \a RALLOC. */
-#define DEALLOC(ptr) \
-  do { sharp_free_(ptr); (ptr)=NULL; } while(0)
-
-#define ALLOC2D(ptr,type,num1,num2) \
-  do { \
-    size_t cnt_, num1_=(num1), num2_=(num2); \
-    ALLOC((ptr),type *,num1_); \
-    ALLOC((ptr)[0],type,num1_*num2_); \
-    for (cnt_=1; cnt_<num1_; ++cnt_) \
-      (ptr)[cnt_]=(ptr)[cnt_-1]+num2_; \
-    } while(0)
-#define DEALLOC2D(ptr) \
-  do { if(ptr) DEALLOC((ptr)[0]); DEALLOC(ptr); } while(0)
-
-#define FAPPROX(a,b,eps) \
-  (fabs((a)-(b))<((eps)*fabs(b)))
-
-#endif

libsharp2/sharp_ylmgen.cc

@@ -29,7 +29,6 @@
 #include <cstdlib>
 #include <algorithm>
 #include "libsharp2/sharp_ylmgen.h"
-#include "libsharp2/sharp_utils.h"
 #include "mr_util/error_handling.h"
 
 using namespace std;

libsharp2/sharp_ylmgen.h

@@ -36,7 +36,7 @@ static constexpr double sharp_fbig=0x1p+800,sharp_fsmall=0x1p-800;
 static constexpr double sharp_ftol=0x1p-60;
 static constexpr double sharp_fbighalf=0x1p+400;
 
-typedef struct { double a, b; } sharp_ylmgen_dbl2;
+struct sharp_ylmgen_dbl2 { double a, b; };
 
 class sharp_Ylmgen
   {

mr_util/math_utils.h

+/*
+ *  This file is part of libcxxsupport.
+ *
+ *  libcxxsupport 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.
+ *
+ *  libcxxsupport 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 libcxxsupport; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*! \file math_utils.h
+ *  Various convenience mathematical functions.
+ *
+ *  Copyright (C) 2002-2015 Max-Planck-Society
+ *  \author Martin Reinecke
+ */
+
+#ifndef MRUTIL_MATH_UTILS_H
+#define MRUTIL_MATH_UTILS_H
+
+
+namespace mr {
+
+/*! Returns \e true if | \a a-b | <= \a epsilon * | \a b |, else \e false. */
+template<typename F> inline bool approx (F a, F b, F epsilon=1e-5)
+  {
+  using namespace std;
+  return abs(a-b) <= (epsilon*abs(b));
+  }
+
+}
+
+#endif

test/sharp2_testsuite.cc

@@ -34,15 +34,15 @@ using std::complex;
 #include "libsharp2/sharp_mpi.h"
 #endif
 #include "libsharp2/sharp.h"
-#include "libsharp2/sharp_utils.h"
-#include "libsharp2/sharp_utils.cc"
 #include "libsharp2/sharp_geomhelpers.h"
 #include "libsharp2/sharp_almhelpers.h"
 #include "mr_util/system.h"
 #include "mr_util/error_handling.h"
 #include "mr_util/threading.h"
+#include "mr_util/math_utils.h"
 
 using namespace std;
+using namespace mr;
 
 static void threading_status(void)
   {
@@ -400,72 +400,69 @@ static void check_sign_scale(void)
   sharp_make_triangular_alm_info(lmax,mmax,1,&alms);
   ptrdiff_t nalms = ((mmax+1)*(mmax+2))/2 + (mmax+1)*(lmax-mmax);
 
-  double **map;
-  ALLOC2D(map,double,2,npix);
+  vector<double> bmap(2*npix);
+  vector<double *>map({&bmap[0], &bmap[npix]});
 
-  dcmplx **alm;
-  ALLOC2D(alm,dcmplx,2,nalms);
+  vector<dcmplx> balm(2*nalms);
+  vector<dcmplx *>alm({&balm[0], &balm[nalms]});
   for (int i=0; i<2; ++i)
     for (int j=0; j<nalms; ++j)
       alm[i][j]=dcmplx(1.,1.);
 
-  sharp_execute(SHARP_ALM2MAP,0,&alm[0],&map[0],tinfo,alms,SHARP_DP,
+  sharp_execute(SHARP_ALM2MAP,0,alm.data(),map.data(),tinfo,alms,SHARP_DP,
     NULL,NULL);
-  MR_assert(FAPPROX(map[0][0     ], 3.588246976618616912e+00,1e-12),
+  MR_assert(approx(map[0][0     ], 3.588246976618616912e+00,1e-12),
     "error");
-  MR_assert(FAPPROX(map[0][npix/2], 4.042209792157496651e+01,1e-12),
+  MR_assert(approx(map[0][npix/2], 4.042209792157496651e+01,1e-12),
     "error");
-  MR_assert(FAPPROX(map[0][npix-1],-1.234675107554816442e+01,1e-12),
+  MR_assert(approx(map[0][npix-1],-1.234675107554816442e+01,1e-12),
     "error");
 
-  sharp_execute(SHARP_ALM2MAP,1,&alm[0],&map[0],tinfo,alms,SHARP_DP,
+  sharp_execute(SHARP_ALM2MAP,1,alm.data(),map.data(),tinfo,alms,SHARP_DP,
     NULL,NULL);
-  MR_assert(FAPPROX(map[0][0     ], 2.750897760535633285e+00,1e-12),
+  MR_assert(approx(map[0][0     ], 2.750897760535633285e+00,1e-12),
     "error");
-  MR_assert(FAPPROX(map[0][npix/2], 3.137704477368562905e+01,1e-12),
+  MR_assert(approx(map[0][npix/2], 3.137704477368562905e+01,1e-12),
     "error");
-  MR_assert(FAPPROX(map[0][npix-1],-8.405730859837063917e+01,1e-12),
+  MR_assert(approx(map[0][npix-1],-8.405730859837063917e+01,1e-12),
     "error");
-  MR_assert(FAPPROX(map[1][0     ],-2.398026536095463346e+00,1e-12),
+  MR_assert(approx(map[1][0     ],-2.398026536095463346e+00,1e-12),
     "error");
-  MR_assert(FAPPROX(map[1][npix/2],-4.961140548331700728e+01,1e-12),
+  MR_assert(approx(map[1][npix/2],-4.961140548331700728e+01,1e-12),
     "error");
-  MR_assert(FAPPROX(map[1][npix-1],-1.412765834230440021e+01,1e-12),
+  MR_assert(approx(map[1][npix-1],-1.412765834230440021e+01,1e-12),
     "error");
 
-  sharp_execute(SHARP_ALM2MAP,2,&alm[0],&map[0],tinfo,alms,SHARP_DP,
+  sharp_execute(SHARP_ALM2MAP,2,alm.data(),map.data(),tinfo,alms,SHARP_DP,
     NULL,NULL);
-  MR_assert(FAPPROX(map[0][0     ],-1.398186224727334448e+00,1e-12),
+  MR_assert(approx(map[0][0     ],-1.398186224727334448e+00,1e-12),
     "error");
-  MR_assert(FAPPROX(map[0][npix/2],-2.456676000884031197e+01,1e-12),
+  MR_assert(approx(map[0][npix/2],-2.456676000884031197e+01,1e-12),
     "error");
-  MR_assert(FAPPROX(map[0][npix-1],-1.516249174408820863e+02,1e-12),
+  MR_assert(approx(map[0][npix-1],-1.516249174408820863e+02,1e-12),
     "error");
-  MR_assert(FAPPROX(map[1][0     ],-3.173406200299964119e+00,1e-12),
+  MR_assert(approx(map[1][0     ],-3.173406200299964119e+00,1e-12),
     "error");
-  MR_assert(FAPPROX(map[1][npix/2],-5.831327404513146462e+01,1e-12),
+  MR_assert(approx(map[1][npix/2],-5.831327404513146462e+01,1e-12),
     "error");
-  MR_assert(FAPPROX(map[1][npix-1],-1.863257892248353897e+01,1e-12),
+  MR_assert(approx(map[1][npix-1],-1.863257892248353897e+01,1e-12),
     "error");
 
-  sharp_execute(SHARP_ALM2MAP_DERIV1,1,&alm[0],&map[0],tinfo,alms,
+  sharp_execute(SHARP_ALM2MAP_DERIV1,1,alm.data(),map.data(),tinfo,alms,
     SHARP_DP,NULL,NULL);
-  MR_assert(FAPPROX(map[0][0     ],-6.859393905369091105e-01,1e-11),
+  MR_assert(approx(map[0][0     ],-6.859393905369091105e-01,1e-11),
     "error");
-  MR_assert(FAPPROX(map[0][npix/2],-2.103947835973212364e+02,1e-12),
+  MR_assert(approx(map[0][npix/2],-2.103947835973212364e+02,1e-12),
     "error");
-  MR_assert(FAPPROX(map[0][npix-1],-1.092463246472086439e+03,1e-12),
+  MR_assert(approx(map[0][npix-1],-1.092463246472086439e+03,1e-12),
     "error");
-  MR_assert(FAPPROX(map[1][0     ],-1.411433220713928165e+02,1e-12),
+  MR_assert(approx(map[1][0     ],-1.411433220713928165e+02,1e-12),
     "error");
-  MR_assert(FAPPROX(map[1][npix/2],-1.146122859381925082e+03,1e-12),
+  MR_assert(approx(map[1][npix/2],-1.146122859381925082e+03,1e-12),
     "error");
-  MR_assert(FAPPROX(map[1][npix-1], 7.821618677689795049e+02,1e-12),
+  MR_assert(approx(map[1][npix-1], 7.821618677689795049e+02,1e-12),
     "error");
 
-  DEALLOC2D(map);
-  DEALLOC2D(alm);
-
   sharp_destroy_alm_info(alms);
   sharp_destroy_geom_info(tinfo);
   }
@@ -479,15 +476,18 @@ static void do_sht (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
   int ncomp = (spin==0) ? 1 : 2;
 
   size_t npix = get_npix(ginfo);
-  double **map;
-  ALLOC2D(map,double,ncomp*ntrans,npix);
+  vector<double> bmap(ncomp*ntrans*npix, 0.);
+  vector<double *>map(ncomp*ntrans);
   for (int i=0; i<ncomp*ntrans; ++i)
-    fill(map[i],map[i]+npix,0);
+    map[i]=&bmap[i*npix];
 
-  dcmplx **alm;
-  ALLOC2D(alm,dcmplx,ncomp*ntrans,nalms);
+  vector<dcmplx> balm(ncomp*ntrans*nalms);
+  vector<dcmplx *>alm(ncomp*ntrans);
   for (int i=0; i<ncomp*ntrans; ++i)
+    {
+    alm[i] = &balm[i*nalms];
     random_alm(alm[i],ainfo,spin,i+1);
+    }
 
   double tta2m, ttm2a;
   unsigned long long toa2m, tom2a;
@@ -497,8 +497,8 @@ static void do_sht (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
   for (size_t itrans=0; itrans<ntrans; ++itrans)
     {
 #ifdef USE_MPI
-    sharp_execute_mpi(MPI_COMM_WORLD,SHARP_ALM2MAP,spin,&alm[itrans*ncomp],
-      &map[itrans*ncomp],ginfo,ainfo, SHARP_DP|SHARP_ADD,&tta2m,&toa2m);
+    sharp_execute_mpi(MPI_COMM_WORLD,SHARP_ALM2MAP,spin,alm[itrans*ncomp],
+      map[itrans*ncomp],ginfo,ainfo, SHARP_DP|SHARP_ADD,&tta2m,&toa2m);
 #else
     sharp_execute(SHARP_ALM2MAP,spin,&alm[itrans*ncomp],&map[itrans*ncomp],ginfo,ainfo,
       SHARP_DP,&tta2m,&toa2m);
@@ -506,7 +506,7 @@ static void do_sht (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
     if (t_a2m!=NULL) *t_a2m+=maxTime(tta2m);
     if (op_a2m!=NULL) *op_a2m+=totalops(toa2m);
     }
-  auto sqsum=get_sqsum_and_invert(alm,nalms,ntrans*ncomp);
+  auto sqsum=get_sqsum_and_invert(alm.data(),nalms,ntrans*ncomp);
   if (t_m2a!=NULL) *t_m2a=0;
   if (op_m2a!=NULL) *op_m2a=0;
   for (size_t itrans=0; itrans<ntrans; ++itrans)
@@ -521,10 +521,7 @@ static void do_sht (sharp_geom_info *ginfo, sharp_alm_info *ainfo,
     if (t_m2a!=NULL) *t_m2a+=maxTime(ttm2a);
     if (op_m2a!=NULL) *op_m2a+=totalops(tom2a);
     }
-  get_errors(alm, nalms, ntrans*ncomp, sqsum, err_abs, err_rel);
-
-  DEALLOC2D(map);
-  DEALLOC2D(alm);
+  get_errors(alm.data(), nalms, ntrans*ncomp, sqsum, err_abs, err_rel);
   }
 
 static void check_accuracy (sharp_geom_info *ginfo, sharp_alm_info *ainfo,