Added sharp_make_mmajor_real_packed_alm_info

libsharp/sharp.c

@@ -859,6 +859,8 @@ static void sharp_build_job_common (sharp_job *job, sharp_jobtype type,
   job->flags = flags;
   if ((job->flags&SHARP_NVMAX)==0)
     job->flags|=sharp_nv_oracle (type, spin, ntrans);
+  if (alm_info->flags&SHARP_REAL_HARMONICS)
+    job->flags|=SHARP_REAL_HARMONICS;
   job->time = 0.;
   job->opcnt = 0;
   job->ntrans = ntrans;

libsharp/sharp_almhelpers.c

@@ -68,3 +68,27 @@ void sharp_make_rectangular_alm_info (int lmax, int mmax, int stride,
     }
   *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 = RALLOC(sharp_alm_info,1);
+  info->lmax = lmax;
+  info->nm = nm;
+  info->mval = RALLOC(int,nm);
+  info->mvstart = RALLOC(ptrdiff_t,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;
+  }

libsharp/sharp_almhelpers.h

@@ -50,6 +50,14 @@ void sharp_make_triangular_alm_info (int lmax, int mmax, int stride,
 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);
+
 #ifdef __cplusplus
 }
 #endif

libsharp/sharp_lowlevel.h

@@ -86,10 +86,21 @@ typedef struct
   } sharp_alm_info;
 
 /*! alm_info flags */
-typedef enum { SHARP_PACKED = 1
+typedef enum { SHARP_PACKED = 1,
-               /*< m=0-coefficients are packed so that the (zero) imaginary part is
+               /*!< m=0-coefficients are packed so that the (zero) imaginary part is
-                   not present. mvstart is in units of *real* float/double for all
+                    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 */
+                    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).
+                */
              } sharp_almflags;
 
 
@@ -172,17 +183,10 @@ typedef enum { SHARP_DP              = 1<<4,
                SHARP_ADD             = 1<<5,
                /*!< results are added to the output arrays, instead of
                     overwriting them */
-               SHARP_REAL_HARMONICS  = 1<<6,
+
-               /*!< Use the real spherical harmonic convention. For
+               /* NOTE: SHARP_REAL_HARMONICS, 1<<6, is also available in sharp_jobflags,
-                    m==0, the alm are treated exactly the same as in
+                  but its use here is deprecated in favor of having it in the sharp_alm_info */
-                    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).
-                */
                SHARP_NO_FFT          = 1<<7,
 
                SHARP_USE_WEIGHTS     = 1<<20,    /* internal use only */