fewer pointers

libsharp2/sharp.cc

@@ -112,38 +112,33 @@ struct ringhelper
     }
   };
 
-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)
-  {
-  sharp_alm_info *info = new sharp_alm_info;
-  info->lmax = lmax;
-  info->nm = nm;
-  info->mval.resize(nm);
-  info->mvstart.resize(nm);
-  info->stride = stride;
-  info->flags = flags;
+sharp_alm_info::sharp_alm_info (int lmax, int nm, int stride, const int *mval,
+  const ptrdiff_t *mstart, int flags)
+  : lmax(lmax), nm(nm), mval(nm), mvstart(nm), stride(stride), flags(flags)
+  {
   for (int mi=0; mi<nm; ++mi)
     {
-    info->mval[mi] = mval[mi];
-    info->mvstart[mi] = mstart[mi];
+    this->mval[mi] = mval[mi];
+    this->mvstart[mi] = mstart[mi];
     }
-  *alm_info = info;
   }
 
-void sharp_make_alm_info (int lmax, int mmax, int stride,
-  const ptrdiff_t *mstart, sharp_alm_info **alm_info)
+sharp_alm_info::sharp_alm_info (int lmax, int mmax, int stride,
+  const ptrdiff_t *mstart)
+  : lmax(lmax), nm(mmax+1), mval(mmax+1), mvstart(mmax+1), stride(stride), flags(0)
   {
-  vector<int> mval(mmax+1);
   for (int i=0; i<=mmax; ++i)
+    {
     mval[i]=i;
-  sharp_make_general_alm_info (lmax, mmax+1, stride, mval.data(), mstart, 0, alm_info);
+    this->mvstart[i] = mstart[i];
+    }
   }
 
-ptrdiff_t sharp_alm_index (const sharp_alm_info *self, int l, int mi)
+ptrdiff_t sharp_alm_info::index (int l, int mi)
   {
-  MR_assert(!(self->flags & SHARP_PACKED),
+  MR_assert(!(flags & SHARP_PACKED),
     "sharp_alm_index not applicable with SHARP_PACKED alms");
-  return self->mvstart[mi]+self->stride*l;
+  return mvstart[mi]+stride*l;
   }
 
 ptrdiff_t sharp_alm_count(const sharp_alm_info *self)
@@ -159,14 +154,9 @@ ptrdiff_t sharp_alm_count(const sharp_alm_info *self)
   return result;
   }
 
-void sharp_destroy_alm_info (sharp_alm_info *info)
-  {
-  delete info;
-  }
-
-void sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
+unique_ptr<sharp_geom_info> sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
   const int *stride, const double *phi0, const double *theta,
-  const double *wgt, sharp_geom_info **geom_info)
+  const double *wgt)
   {
   sharp_geom_info *info = new sharp_geom_info;
   vector<sharp_ringinfo> infos(nrings);
@@ -175,7 +165,6 @@ void sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
   info->pair.resize(nrings);
   int npairs=0;
   info->nphmax=0;
-  *geom_info = info;
 
   for (int m=0; m<nrings; ++m)
     {
@@ -220,22 +209,7 @@ void sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
         (a.r1.cth>b.r1.cth));
     return a.r1.nph<b.r1.nph;
     });
-  }
-
-ptrdiff_t sharp_map_size(const sharp_geom_info *info)
-  {
-  ptrdiff_t result = 0;
-  for (int m=0; m<info->pair.size(); ++m)
-    {
-      result+=info->pair[m].r1.nph;
-      result+=(info->pair[m].r2.nph>=0) ? (info->pair[m].r2.nph) : 0;
-    }
-  return result;
-  }
-
-void sharp_destroy_geom_info (sharp_geom_info *geom_info)
-  {
-  delete geom_info;
+  return unique_ptr<sharp_geom_info>(info);
   }
 
 /* This currently requires all m values from 0 to nm-1 to be present.
@@ -806,49 +780,49 @@ MRUTIL_NOINLINE static void phase2map (sharp_job *job, int mmax, int llim, int u
     }
   }
 
-MRUTIL_NOINLINE static void sharp_execute_job (sharp_job *job)
+MRUTIL_NOINLINE static void sharp_execute_job (sharp_job &job)
   {
   mr::timers::SimpleTimer timer;
-  job->opcnt=0;
-  int lmax = job->ainfo->lmax,
-      mmax=sharp_get_mmax(job->ainfo->mval);
+  job.opcnt=0;
+  int lmax = job.ainfo->lmax,
+      mmax=sharp_get_mmax(job.ainfo->mval);
 
-  job->norm_l = (job->type==SHARP_ALM2MAP_DERIV1) ?
+  job.norm_l = (job.type==SHARP_ALM2MAP_DERIV1) ?
      sharp_Ylmgen::get_d1norm (lmax) :
-     sharp_Ylmgen::get_norm (lmax, job->spin);
+     sharp_Ylmgen::get_norm (lmax, job.spin);
 
 /* clear output arrays if requested */
-  init_output (job);
+  init_output (&job);
 
   int nchunks, chunksize;
-  get_chunk_info(job->ginfo->pair.size(),sharp_veclen()*sharp_max_nvec(job->spin),
+  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, phasebuffer);
+  alloc_phase (&job,mmax+1,chunksize, phasebuffer);
   std::atomic<size_t> opcnt = 0;
 
 /* chunk loop */
   for (int chunk=0; chunk<nchunks; ++chunk)
     {
-    int llim=chunk*chunksize, ulim=min<int>(llim+chunksize,job->ginfo->pair.size());
+    int llim=chunk*chunksize, ulim=min<int>(llim+chunksize,job.ginfo->pair.size());
     vector<int> ispair(ulim-llim);
     vector<int> mlim(ulim-llim);
     vector<double> cth(ulim-llim), sth(ulim-llim);
     for (int i=0; i<ulim-llim; ++i)
       {
-      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]);
+      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]);
       }
 
 /* map->phase where necessary */
-    map2phase (job, mmax, llim, ulim);
+    map2phase (&job, mmax, llim, ulim);
 
-    mr::execDynamic(job->ainfo->nm, 0, 1, [&](mr::Scheduler &sched)
+    mr::execDynamic(job.ainfo->nm, 0, 1, [&](mr::Scheduler &sched)
       {
-      sharp_job ljob = *job;
+      sharp_job ljob = job;
       ljob.opcnt=0;
       sharp_Ylmgen generator(lmax,mmax,ljob.spin);
       aligned_array<dcmplx> almbuffer;
@@ -869,47 +843,47 @@ MRUTIL_NOINLINE static void sharp_execute_job (sharp_job *job)
       }); /* end of parallel region */
 
 /* phase->map where necessary */
-    phase2map (job, mmax, llim, ulim);
+    phase2map (&job, mmax, llim, ulim);
     } /* end of chunk loop */
 
-  job->opcnt = opcnt;
-  job->time=timer();
+  job.opcnt = opcnt;
+  job.time=timer();
   }
 
-static void sharp_build_job_common (sharp_job *job, sharp_jobtype type,
-  int spin, void *alm, void *map, const sharp_geom_info *geom_info,
-  const sharp_alm_info *alm_info, int flags)
+static void sharp_build_job_common (sharp_job &job, sharp_jobtype type,
+  int spin, void *alm, void *map, const sharp_geom_info &geom_info,
+  const sharp_alm_info &alm_info, int flags)
   {
   if (type==SHARP_ALM2MAP_DERIV1) spin=1;
   if (type==SHARP_MAP2ALM) flags|=SHARP_USE_WEIGHTS;
   if (type==SHARP_Yt) type=SHARP_MAP2ALM;
   if (type==SHARP_WY) { type=SHARP_ALM2MAP; flags|=SHARP_USE_WEIGHTS; }
 
-  MR_assert((spin>=0)&&(spin<=alm_info->lmax), "bad spin");
-  job->type = type;
-  job->spin = spin;
-  job->nmaps = (type==SHARP_ALM2MAP_DERIV1) ? 2 : ((spin>0) ? 2 : 1);
-  job->nalm = (type==SHARP_ALM2MAP_DERIV1) ? 1 : ((spin>0) ? 2 : 1);
-  job->ginfo = geom_info;
-  job->ainfo = alm_info;
-  job->flags = flags;
-  if (alm_info->flags&SHARP_REAL_HARMONICS)
-    job->flags|=SHARP_REAL_HARMONICS;
-  job->time = 0.;
-  job->opcnt = 0;
-  job->alm=(void **)alm;
-  job->map=(void **)map;
+  MR_assert((spin>=0)&&(spin<=alm_info.lmax), "bad spin");
+  job.type = type;
+  job.spin = spin;
+  job.nmaps = (type==SHARP_ALM2MAP_DERIV1) ? 2 : ((spin>0) ? 2 : 1);
+  job.nalm = (type==SHARP_ALM2MAP_DERIV1) ? 1 : ((spin>0) ? 2 : 1);
+  job.ginfo = &geom_info;
+  job.ainfo = &alm_info;
+  job.flags = flags;
+  if (alm_info.flags&SHARP_REAL_HARMONICS)
+    job.flags|=SHARP_REAL_HARMONICS;
+  job.time = 0.;
+  job.opcnt = 0;
+  job.alm=(void **)alm;
+  job.map=(void **)map;
   }
 
 void sharp_execute (sharp_jobtype type, int spin, void *alm, void *map,
-  const sharp_geom_info *geom_info, const sharp_alm_info *alm_info,
+  const sharp_geom_info &geom_info, const sharp_alm_info &alm_info,
   int flags, double *time, unsigned long long *opcnt)
   {
   sharp_job job;
-  sharp_build_job_common (&job, type, spin, alm, map, geom_info, alm_info,
+  sharp_build_job_common (job, type, spin, alm, map, geom_info, alm_info,
     flags);
 
-  sharp_execute_job (&job);
+  sharp_execute_job (job);
   if (time!=NULL) *time = job.time;
   if (opcnt!=NULL) *opcnt = job.opcnt;
   }

libsharp2/sharp.h

@@ -30,6 +30,7 @@
 
 #include <cstddef>
 #include <vector>
+#include <memory>
 
 /*! \internal
     Helper type containing information about a single ring. */
@@ -76,27 +77,6 @@ struct sharp_alm_info
   std::vector<ptrdiff_t> mvstart;
   /*! Stride between a_lm and a_(l+1),m */
   ptrdiff_t stride;
-  };
-
-/*! alm_info flags */
-enum sharp_almflags { SHARP_PACKED = 1,
-               /*!< m=0-coefficients are packed so that the (zero) imaginary part is
-                    not present. mvstart is in units of *real* float/double for all
-                    m; stride is in units of reals for m=0 and complex for m!=0 */
-               SHARP_REAL_HARMONICS  = 1<<6
-               /*!< Use the real spherical harmonic convention. For
-                    m==0, the alm are treated exactly the same as in
-                    the complex case.  For m!=0, alm[i] represent a
-                    pair (+abs(m), -abs(m)) instead of (real, imag),
-                    and the coefficients are scaled by a factor of
-                    sqrt(2) relative to the complex case.  In other
-                    words, (sqrt(.5) * alm[i]) recovers the
-                    corresponding complex coefficient (when accessed
-                    as complex).
-                */
-             };
-
-
 
 /*! Creates an a_lm data structure from the following parameters:
     \param lmax maximum \a l quantum number (>=0)
@@ -105,10 +85,9 @@ enum sharp_almflags { SHARP_PACKED = 1,
       differing by 1.
     \param mstart the index of the (hypothetical) coefficient with the
       quantum numbers 0,\a m. Must have \a mmax+1 entries.
-    \param alm_info will hold a pointer to the newly created data structure
  */
-void sharp_make_alm_info (int lmax, int mmax, int stride,
-  const ptrdiff_t *mstart, sharp_alm_info **alm_info);
+  sharp_alm_info(int lmax, int mmax, int stride, const ptrdiff_t *mstart);
+
 /*! Creates an a_lm data structure which from the following parameters:
     \param lmax maximum \a l quantum number (\a >=0)
     \param nm number of different \a m (\a 0<=nm<=lmax+1)
@@ -118,21 +97,33 @@ void sharp_make_alm_info (int lmax, int mmax, int stride,
     \param mvstart array with \a nm entries containing the (hypothetical)
       indices of the coefficients with the quantum numbers 0,\a mval[i]
     \param flags a combination of sharp_almflags (pass 0 unless you know you need this)
-    \param alm_info will hold a pointer to the newly created data structure
  */
-void sharp_make_general_alm_info (int lmax, int nm, int stride, const int *mval,
-  const ptrdiff_t *mvstart, int flags, sharp_alm_info **alm_info);
+  sharp_alm_info (int lmax, int nm, int stride, const int *mval,
+    const ptrdiff_t *mvstart, int flags);
 /*! Returns the index of the coefficient with quantum numbers \a l,
     \a mval[mi].
     \note for a \a sharp_alm_info generated by sharp_make_alm_info() this is
     the index for the coefficient with the quantum numbers \a l, \a mi. */
-ptrdiff_t sharp_alm_index (const sharp_alm_info *self, int l, int mi);
-/*! Returns the number of alm coefficients described by \a self. If the SHARP_PACKED
-    flag is set, this is number of "real" coeffecients (for m < 0 and m >= 0),
-    otherwise it is the number of complex coefficients (with m>=0). */
-ptrdiff_t sharp_alm_count(const sharp_alm_info *self);
-/*! Deallocates the a_lm info object. */
-void sharp_destroy_alm_info (sharp_alm_info *info);
+  ptrdiff_t index (int l, int mi);
+  };
+
+/*! alm_info flags */
+enum sharp_almflags { SHARP_PACKED = 1,
+               /*!< m=0-coefficients are packed so that the (zero) imaginary part is
+                    not present. mvstart is in units of *real* float/double for all
+                    m; stride is in units of reals for m=0 and complex for m!=0 */
+               SHARP_REAL_HARMONICS  = 1<<6
+               /*!< Use the real spherical harmonic convention. For
+                    m==0, the alm are treated exactly the same as in
+                    the complex case.  For m!=0, alm[i] represent a
+                    pair (+abs(m), -abs(m)) instead of (real, imag),
+                    and the coefficients are scaled by a factor of
+                    sqrt(2) relative to the complex case.  In other
+                    words, (sqrt(.5) * alm[i]) recovers the
+                    corresponding complex coefficient (when accessed
+                    as complex).
+                */
+             };
 
 /*! \} */
 
@@ -152,17 +143,9 @@ void sharp_destroy_alm_info (sharp_alm_info *info);
       Pass NULL to use 1.0 as weight for all rings.
     \param geom_info will hold a pointer to the newly created data structure
  */
-void sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
+std::unique_ptr<sharp_geom_info> sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
   const int *stride, const double *phi0, const double *theta,
-  const double *wgt, sharp_geom_info **geom_info);
-
-/*! Counts the number of grid points needed for (the local part of) a map described
-    by \a info.
- */
-ptrdiff_t sharp_map_size(const sharp_geom_info *info);
-
-/*! Deallocates the geometry information in \a info. */
-void sharp_destroy_geom_info (sharp_geom_info *info);
+  const double *wgt);
 
 /*! \} */
 
@@ -220,7 +203,7 @@ enum sharp_jobflags { SHARP_DP              = 1<<4,
   \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, void *alm, void *map,
-  const sharp_geom_info *geom_info, const sharp_alm_info *alm_info,
+  const sharp_geom_info &geom_info, const sharp_alm_info &alm_info,
   int flags, double *time, unsigned long long *opcnt);
 
 void sharp_set_chunksize_min(int new_chunksize_min);

libsharp2/sharp_almhelpers.cc

@@ -27,63 +27,22 @@
 
 #include "libsharp2/sharp_almhelpers.h"
 
-void sharp_make_triangular_alm_info (int lmax, int mmax, int stride,
-  sharp_alm_info **alm_info)
+using namespace std;
+
+unique_ptr<sharp_alm_info> sharp_make_triangular_alm_info (int lmax, int mmax, int stride)
   {
-  sharp_alm_info *info = new sharp_alm_info;
-  info->lmax = lmax;
-  info->nm = mmax+1;
-  info->mval.resize(mmax+1);
-  info->mvstart.resize(mmax+1);
-  info->stride = stride;
-  info->flags = 0;
+  vector<ptrdiff_t> mvstart(mmax+1);
   ptrdiff_t tval = 2*lmax+1;
   for (ptrdiff_t m=0; m<=mmax; ++m)
-    {
-    info->mval[m] = m;
-    info->mvstart[m] = stride*((m*(tval-m))>>1);
-    }
-  *alm_info = info;
+    mvstart[m] = stride*((m*(tval-m))>>1);
+  return make_unique<sharp_alm_info>(lmax, mmax, stride, mvstart.data());
   }
 
-void sharp_make_rectangular_alm_info (int lmax, int mmax, int stride,
-  sharp_alm_info **alm_info)
+unique_ptr<sharp_alm_info> sharp_make_rectangular_alm_info (int lmax, int mmax, int stride)
   {
-  sharp_alm_info *info = new sharp_alm_info;
-  info->lmax = lmax;
-  info->nm = mmax+1;
-  info->mval.resize(mmax+1);
-  info->mvstart.resize(mmax+1);
-  info->stride = stride;
-  info->flags = 0;
+  vector<ptrdiff_t> mvstart(mmax+1);
+  ptrdiff_t tval = 2*lmax+1;
   for (ptrdiff_t m=0; m<=mmax; ++m)
-    {
-    info->mval[m] = m;
-    info->mvstart[m] = stride*m*(lmax+1);
-    }
-  *alm_info = info;
-  }
-
-void sharp_make_mmajor_real_packed_alm_info (int lmax, int stride,
-  int nm, const int *ms, sharp_alm_info **alm_info)
-  {
-  ptrdiff_t idx;
-  int f;
-  sharp_alm_info *info = new sharp_alm_info;
-  info->lmax = lmax;
-  info->nm = nm;
-  info->mval.resize(nm);
-  info->mvstart.resize(nm);
-  info->stride = stride;
-  info->flags = SHARP_PACKED | SHARP_REAL_HARMONICS;
-  idx = 0;  /* tracks the number of 'consumed' elements so far; need to correct by m */
-  for (int im=0; im!=nm; ++im)
-    {
-    int m=(ms==NULL)?im:ms[im];
-    f = (m==0) ? 1 : 2;
-    info->mval[im] = m;
-    info->mvstart[im] = stride * (idx - f * m);
-    idx += f * (lmax + 1 - m);
-    }
-  *alm_info = info;
+    mvstart[m] = stride*m*(lmax+1);
+  return make_unique<sharp_alm_info>(lmax, mmax, stride, mvstart.data());
   }

libsharp2/sharp_almhelpers.h

@@ -28,26 +28,18 @@
 #ifndef SHARP2_ALMHELPERS_H
 #define SHARP2_ALMHELPERS_H
 
+
+#include <memory>
 #include "libsharp2/sharp.h"
 
 /*! Initialises an a_lm data structure according to the scheme used by
     Healpix_cxx.
     \ingroup almgroup */
-void sharp_make_triangular_alm_info (int lmax, int mmax, int stride,
-  sharp_alm_info **alm_info);
+std::unique_ptr<sharp_alm_info> sharp_make_triangular_alm_info (int lmax, int mmax, int stride);
 
 /*! Initialises an a_lm data structure according to the scheme used by
     Fortran Healpix
     \ingroup almgroup */
-void sharp_make_rectangular_alm_info (int lmax, int mmax, int stride,
-  sharp_alm_info **alm_info);
-
-/*! Initialises alm_info for mmajor, real, packed spherical harmonics.
-    Pass \a mmax + 1 to nm and NULL to \a ms in order to use everything;
-    otherwise you can pick a subset of m to process (should only be used
-    for MPI parallelization).
-    \ingroup almgroup */
-void sharp_make_mmajor_real_packed_alm_info (int lmax, int stride,
-  int nm, const int *ms, sharp_alm_info **alm_info);
+std::unique_ptr<sharp_alm_info> sharp_make_rectangular_alm_info (int lmax, int mmax, int stride);
 
 #endif

libsharp2/sharp_geomhelpers.cc

@@ -34,8 +34,8 @@
 
 using namespace std;
 
-void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
-  const int *rings, const double *weight, sharp_geom_info **geom_info)
+unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
+  const int *rings, const double *weight)
   {
   const double pi=3.141592653589793238462643383279502884197;
   ptrdiff_t npix=(ptrdiff_t)nside*nside*12;
@@ -84,18 +84,17 @@ void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
     curofs+=nph[m];
     }
 
-  sharp_make_geom_info (nrings, nph.data(), ofs.data(), stride_.data(), phi0.data(), theta.data(), weight_.data(),
+  return sharp_make_geom_info(nrings, nph.data(), ofs.data(), stride_.data(), phi0.data(), theta.data(), weight_.data());