hp_generate_coords: parallelized loop with OpenMP.

src/mpim_grid/grid_healpix.cc

@@ -22,21 +22,28 @@ hp_get_indexing(const char *indexingName)
   return indexing;
 }
 
-inline size_t
-hp_convert_index(HpIndexing indexing, size_t ipix, size_t nside)
+static int64_t
+hp_xy_to_xy(int64_t ipix, int nside)
 {
-  // clang-format off
-  return (indexing == HpIndexing::Ring) ? healpixl_ring_to_xy(ipix, nside) :
-         (indexing == HpIndexing::Nest) ? healpixl_nested_to_xy(ipix, nside) : ipix;
-  // clang-format on
+  (void)nside;
+  return ipix;
 }
 
 void
-hp_generate_coords(HpIndexing indexing, size_t nside, size_t nvals, double *xvals, double *yvals, bool withBounds, double *xbounds, double *ybounds)
+hp_generate_coords(HpIndexing indexing, int nside, int64_t nvals, double *xvals, double *yvals, bool withBounds, double *xbounds, double *ybounds)
 {
-  for (size_t ipix = 0; ipix < nvals; ipix++)
+  int64_t (*convert_index)(int64_t ipix, int nside);
+  // clang-format off
+  convert_index = (indexing == HpIndexing::Ring) ? &healpixl_ring_to_xy :
+                  (indexing == HpIndexing::Nest) ? &healpixl_nested_to_xy : &hp_xy_to_xy;
+  // clang-format on
+
+#ifdef _OPENMP
+#pragma omp parallel for if (nvals > 99999) default(shared) schedule(static)
+#endif
+  for (int64_t ipix = 0; ipix < nvals; ipix++)
     {
-      const auto index = hp_convert_index(indexing, ipix, nside);
+      const auto index = convert_index(ipix, nside);
 
       healpixl_to_radec(index, nside, 0.5, 0.5, &xvals[ipix], &yvals[ipix]);
 

src/mpim_grid/grid_healpix.h

 #ifndef GRID_HEALPIX_H
 #define GRID_HEALPIX_H
 
+#include <cinttypes>
+
 enum class HpIndexing
   {
     Undef,
@@ -10,6 +12,6 @@ enum class HpIndexing
   };
 
 HpIndexing hp_get_indexing(const char *indexingName);
-void hp_generate_coords(HpIndexing indexing, size_t nside, size_t nvals, double *xvals, double *yvals, bool withBounds, double *xbounds, double *ybounds);
+void hp_generate_coords(HpIndexing indexing, int nside, int64_t nvals, double *xvals, double *yvals, bool withBounds, double *xbounds, double *ybounds);
 
 #endif