SHARP_PACKED alm_info flag

libsharp/sharp.c

@@ -124,7 +124,7 @@ static int ringpair_compare (const void *xa, const void *xb)
   }
 
 void sharp_make_general_alm_info (int lmax, int nm, int stride, const int *mval,
-  const ptrdiff_t *mstart, sharp_alm_info **alm_info)
+  const ptrdiff_t *mstart, int flags, sharp_alm_info **alm_info)
   {
   sharp_alm_info *info = RALLOC(sharp_alm_info,1);
   info->lmax = lmax;
@@ -132,6 +132,7 @@ void sharp_make_general_alm_info (int lmax, int nm, int stride, const int *mval,
   info->mval = RALLOC(int,nm);
   info->mvstart = RALLOC(ptrdiff_t,nm);
   info->stride = stride;
+  info->flags = flags;
   for (int mi=0; mi<nm; ++mi)
     {
     info->mval[mi] = mval[mi];
@@ -146,12 +147,16 @@ void sharp_make_alm_info (int lmax, int mmax, int stride,
   int *mval=RALLOC(int,mmax+1);
   for (int i=0; i<=mmax; ++i)
     mval[i]=i;
-  sharp_make_general_alm_info (lmax, mmax+1, stride, mval, mstart, alm_info);
+  sharp_make_general_alm_info (lmax, mmax+1, stride, mval, mstart, 0, alm_info);
   DEALLOC(mval);
   }
 
 ptrdiff_t sharp_alm_index (const sharp_alm_info *self, int l, int mi)
-  { return self->mvstart[mi]+self->stride*l; }
+  {
+  UTIL_ASSERT(!(self->flags & SHARP_PACKED),
+              "sharp_alm_index not applicable with SHARP_PACKED alms");
+  return self->mvstart[mi]+self->stride*l; 
+  }
 
 void sharp_destroy_alm_info (sharp_alm_info *info)
   {
@@ -358,17 +363,40 @@ static void fill_map (const sharp_geom_info *ginfo, void *map, double value,
     }
   }
 
-static void fill_alm (const sharp_alm_info *ainfo, void *alm, dcmplx value,
+static void clear_alm (const sharp_alm_info *ainfo, void *alm, int flags)
-  int flags)
   {
-  if (flags&SHARP_DP)
+#define CLEARLOOP(real_t,body)             \
-    for (int mi=0;mi<ainfo->nm;++mi)
+      {                                    \
-      for (int l=ainfo->mval[mi];l<=ainfo->lmax;++l)
+        real_t *talm = (real_t *)alm;      \
-        ((dcmplx *)alm)[sharp_alm_index(ainfo,l,mi)] = value;
+          for (int l=m;l<=ainfo->lmax;++l) \
-  else
+            body                           \
-    for (int mi=0;mi<ainfo->nm;++mi)
+      }
-      for (int l=ainfo->mval[mi];l<=ainfo->lmax;++l)
+
-        ((fcmplx *)alm)[sharp_alm_index(ainfo,l,mi)] = (fcmplx)value;
+  for (int mi=0;mi<ainfo->nm;++mi)
+    {
+      int m=ainfo->mval[mi];
+      ptrdiff_t mvstart = ainfo->mvstart[mi];
+      ptrdiff_t stride = ainfo->stride;
+      if (!(ainfo->flags&SHARP_PACKED))
+        mvstart*=2;
+      if ((ainfo->flags&SHARP_PACKED)&&(m==0))
+        {
+        if (flags&SHARP_DP)
+          CLEARLOOP(double, talm[mvstart+l*stride] = 0.;)
+        else
+          CLEARLOOP(float, talm[mvstart+l*stride] = 0.;)
+        }
+      else
+        {
+        stride*=2;
+        if (flags&SHARP_DP)
+          CLEARLOOP(double, talm[mvstart+l*stride] = talm[mvstart+l*stride+1] = 0.;)
+        else
+          CLEARLOOP(float, talm[mvstart+l*stride] = talm[mvstart+l*stride+1] = 0.;)
+        }
+
+#undef CLEARLOOP
+    }
   }
 
 static void init_output (sharp_job *job)
@@ -376,7 +404,7 @@ static void init_output (sharp_job *job)
   if (job->flags&SHARP_ADD) return;
   if (job->type == SHARP_MAP2ALM)
     for (int i=0; i<job->ntrans*job->nalm; ++i)
-      fill_alm (job->ainfo,job->alm[i],0.,job->flags);
+      clear_alm (job->ainfo,job->alm[i],job->flags);
   else
     for (int i=0; i<job->ntrans*job->nmaps; ++i)
       fill_map (job->ginfo,job->map[i],0.,job->flags);
@@ -418,11 +446,11 @@ static void dealloc_almtmp (sharp_job *job)
 static void alm2almtmp (sharp_job *job, int lmax, int mi)
   {
 
-#define COPY_LOOP(source_t, expr_of_x)                      \
+#define COPY_LOOP(real_t, source_t, expr_of_x)                      \
   for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)            \
     for (int i=0; i<job->ntrans*job->nalm; ++i)             \
       {                                                     \
-        source_t x = ((source_t *)job->alm[i])[ofs+l*stride]; \
+        source_t x = *(source_t *)(((real_t *)job->alm[i])+ofs+l*stride); \
         job->almtmp[job->ntrans*job->nalm*l+i] = expr_of_x; \
       }
 
@@ -430,19 +458,44 @@ static void alm2almtmp (sharp_job *job, int lmax, int mi)
     {
     ptrdiff_t ofs=job->ainfo->mvstart[mi];
     int stride=job->ainfo->stride;
+    int packed_m0=(job->ainfo->flags&SHARP_PACKED)&&(job->ainfo->mval[mi]==0);
+    if (!(job->ainfo->flags&SHARP_PACKED))
+      ofs *= 2;
+    if (!packed_m0)
+      stride *= 2;
     if (job->spin==0)
       {
-      if (job->flags&SHARP_DP)
+      if (packed_m0)
-        COPY_LOOP(dcmplx, x)
+        {
+        if (job->flags&SHARP_DP)
+          COPY_LOOP(double, double, x)
+        else
+          COPY_LOOP(float, float, x)
+        }
       else
-        COPY_LOOP(fcmplx, x)
+        {
+        if (job->flags&SHARP_DP)
+          COPY_LOOP(double, dcmplx, x)
+        else
+          COPY_LOOP(float, fcmplx, x)
+        }
       }
     else
       {
-      if (job->flags&SHARP_DP)
+      if (packed_m0)
-        COPY_LOOP(dcmplx, x*job->norm_l[l])
+        {
+        if (job->flags&SHARP_DP)
+          COPY_LOOP(double, double, x*job->norm_l[l])
+        else
+          COPY_LOOP(float, float, x*job->norm_l[l])
+        }
       else
-        COPY_LOOP(fcmplx, x*job->norm_l[l])
+        {
+        if (job->flags&SHARP_DP)
+          COPY_LOOP(double, dcmplx, x*job->norm_l[l])
+        else
+          COPY_LOOP(float, fcmplx, x*job->norm_l[l])
+        }
       }
     }
   else
@@ -455,30 +508,55 @@ static void alm2almtmp (sharp_job *job, int lmax, int mi)
 static void almtmp2alm (sharp_job *job, int lmax, int mi)
   {
 
-#define COPY_LOOP(target_t, expr_of_x)                       \
+#define COPY_LOOP(real_t, target_t, expr_of_x)               \
   for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)             \
     for (int i=0; i<job->ntrans*job->nalm; ++i)              \
       {                                                      \
         dcmplx x = job->almtmp[job->ntrans*job->nalm*l+i];   \
-        ((target_t *)job->alm[i])[ofs+l*stride] += expr_of_x; \
+        *(target_t *)(((real_t *)job->alm[i])+ofs+l*stride) += expr_of_x; \
       }
 
   if (job->type != SHARP_MAP2ALM) return;
   ptrdiff_t ofs=job->ainfo->mvstart[mi];
   int stride=job->ainfo->stride;
+  int packed_m0=(job->ainfo->flags&SHARP_PACKED)&&(job->ainfo->mval[mi]==0);
+  if (!(job->ainfo->flags&SHARP_PACKED))
+    ofs *= 2;
+  if (!packed_m0)
+    stride *= 2;
   if (job->spin==0)
     {
-    if (job->flags&SHARP_DP)
+    if (packed_m0)
-      COPY_LOOP(dcmplx, x)
+      {
+      if (job->flags&SHARP_DP)
+        COPY_LOOP(double, double, creal(x))
+      else
+        COPY_LOOP(float, float, crealf(x))
+      }
     else
-      COPY_LOOP(fcmplx, (fcmplx)x)
+      {
+      if (job->flags&SHARP_DP)
+        COPY_LOOP(double, dcmplx, x)
+      else
+        COPY_LOOP(float, fcmplx, (fcmplx)x)
+      }
     }
   else
     {
-    if (job->flags&SHARP_DP)
+    if (packed_m0)
-      COPY_LOOP(dcmplx, x * job->norm_l[l])
+      {
+      if (job->flags&SHARP_DP)
+        COPY_LOOP(double, double, creal(x) * job->norm_l[l])
+      else
+        COPY_LOOP(float, fcmplx, (float)(creal(x) * job->norm_l[l]))
+      }
     else
-      COPY_LOOP(fcmplx, (fcmplx)(x * job->norm_l[l]))
+      {
+      if (job->flags&SHARP_DP)
+        COPY_LOOP(double, dcmplx, x * job->norm_l[l])
+      else
+        COPY_LOOP(float, fcmplx, (fcmplx)(x * job->norm_l[l]))
+      }
     }
 
 #undef COPY_LOOP

libsharp/sharp_almhelpers.c

@@ -41,6 +41,7 @@ void sharp_make_triangular_alm_info (int lmax, int mmax, int stride,
   info->mval = RALLOC(int,mmax+1);
   info->mvstart = RALLOC(ptrdiff_t,mmax+1);
   info->stride = stride;
+  info->flags = 0;
   int tval = 2*lmax+1;
   for (ptrdiff_t m=0; m<=mmax; ++m)
     {
@@ -59,6 +60,7 @@ void sharp_make_rectangular_alm_info (int lmax, int mmax, int stride,
   info->mval = RALLOC(int,mmax+1);
   info->mvstart = RALLOC(ptrdiff_t,mmax+1);
   info->stride = stride;
+  info->flags = 0;
   for (ptrdiff_t m=0; m<=mmax; ++m)
     {
     info->mval[m] = m;

libsharp/sharp_lowlevel.h

@@ -76,6 +76,8 @@ typedef struct
   int nm;
   /*! Array with \a nm entries containing the individual m values */
   int *mval;
+  /*! Combination of flags from sharp_almflags */
+  int flags;
   /*! Array with \a nm entries containing the (hypothetical) indices of
       the coefficients with quantum numbers 0,\a mval[i] */
   ptrdiff_t *mvstart;
@@ -83,6 +85,15 @@ typedef struct
   ptrdiff_t stride;
   } sharp_alm_info;
 
+/*! alm_info flags */
+typedef enum { 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_almflags;
+
+
+
 /*! Creates an a_lm data structure from the following parameters:
     \param lmax maximum \a l quantum number (>=0)
     \param mmax maximum \a m quantum number (0<= \a mmax <= \a lmax)
@@ -102,10 +113,11 @@ void sharp_make_alm_info (int lmax, int mmax, int stride,
     \param mval array with \a nm entries containing the individual m values
     \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, sharp_alm_info **alm_info);
+  const ptrdiff_t *mvstart, int flags, sharp_alm_info **alm_info);
 /*! 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