sharp.h 8.82 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
/*
 *  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 */

/*! \file sharp.h
 *  Portable interface for the spherical transform library.
 *
 *  Copyright (C) 2012-2019 Max-Planck-Society
 *  \author Martin Reinecke \author Dag Sverre Seljebotn
 */

#ifndef SHARP_SHARP_H
#define SHARP_SHARP_H

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
31
32
#include <cstddef>
#include <vector>
Martin Reinecke's avatar
Martin Reinecke committed
33
#include <memory>
Martin Reinecke's avatar
Martin Reinecke committed
34
35
36

/*! \internal
    Helper type containing information about a single ring. */
Martin Reinecke's avatar
Martin Reinecke committed
37
struct sharp_ringinfo
Martin Reinecke's avatar
Martin Reinecke committed
38
39
40
41
  {
  double theta, phi0, weight, cth, sth;
  ptrdiff_t ofs;
  int nph, stride;
Martin Reinecke's avatar
Martin Reinecke committed
42
  };
Martin Reinecke's avatar
Martin Reinecke committed
43
44
45
46

/*! \internal
    Helper type containing information about a pair of rings with colatitudes
    symmetric around the equator. */
Martin Reinecke's avatar
Martin Reinecke committed
47
struct sharp_ringpair
Martin Reinecke's avatar
Martin Reinecke committed
48
49
  {
  sharp_ringinfo r1,r2;
Martin Reinecke's avatar
Martin Reinecke committed
50
  };
Martin Reinecke's avatar
Martin Reinecke committed
51
52
53

/*! \internal
    Type holding all required information about a map geometry. */
Martin Reinecke's avatar
Martin Reinecke committed
54
struct sharp_geom_info
Martin Reinecke's avatar
Martin Reinecke committed
55
  {
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
56
57
  std::vector<sharp_ringpair> pair;
  int nphmax;
Martin Reinecke's avatar
Martin Reinecke committed
58
  };
Martin Reinecke's avatar
Martin Reinecke committed
59
60
61
62
63
64

/*! \defgroup almgroup Helpers for dealing with a_lm */
/*! \{ */

/*! \internal
    Helper type for index calculation in a_lm arrays. */
Martin Reinecke's avatar
Martin Reinecke committed
65
struct sharp_alm_info
Martin Reinecke's avatar
Martin Reinecke committed
66
67
68
69
70
71
  {
  /*! Maximum \a l index of the array */
  int lmax;
  /*! Number of different \a m values in this object */
  int nm;
  /*! Array with \a nm entries containing the individual m values */
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
72
  std::vector<int> mval;
Martin Reinecke's avatar
Martin Reinecke committed
73
74
75
76
  /*! 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] */
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
77
  std::vector<ptrdiff_t> mvstart;
Martin Reinecke's avatar
Martin Reinecke committed
78
79
80
81
82
83
84
85
86
87
88
  /*! Stride between a_lm and a_(l+1),m */
  ptrdiff_t stride;

/*! 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)
    \param stride the stride between entries with identical \a m, and \a l
      differing by 1.
    \param mstart the index of the (hypothetical) coefficient with the
      quantum numbers 0,\a m. Must have \a mmax+1 entries.
 */
Martin Reinecke's avatar
Martin Reinecke committed
89
90
  sharp_alm_info(int lmax, int mmax, int stride, const ptrdiff_t *mstart);

Martin Reinecke's avatar
Martin Reinecke committed
91
92
93
94
95
96
97
98
99
100
/*! 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)
    \param stride the stride between entries with identical \a m, and \a l
      differing by 1.
    \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)
 */
Martin Reinecke's avatar
Martin Reinecke committed
101
102
  sharp_alm_info (int lmax, int nm, int stride, const int *mval,
    const ptrdiff_t *mvstart, int flags);
Martin Reinecke's avatar
Martin Reinecke committed
103
104
105
106
/*! 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. */
Martin Reinecke's avatar
Martin Reinecke committed
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
  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).
                */
             };
Martin Reinecke's avatar
Martin Reinecke committed
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

/*! \} */

/*! \defgroup geominfogroup Functions for dealing with geometry information */
/*! \{ */

/*! Creates a geometry information from a set of ring descriptions.
    All arrays passed to this function must have \a nrings elements.
    \param nrings the number of rings in the map
    \param nph the number of pixels in each ring
    \param ofs the index of the first pixel in each ring in the map array
    \param stride the stride between consecutive pixels
    \param phi0 the azimuth (in radians) of the first pixel in each ring
    \param theta the colatitude (in radians) of each ring
    \param wgt the pixel weight to be used for the ring in map2alm
      and adjoint map2alm transforms.
      Pass NULL to use 1.0 as weight for all rings.
    \param geom_info will hold a pointer to the newly created data structure
 */
Martin Reinecke's avatar
Martin Reinecke committed
146
std::unique_ptr<sharp_geom_info> sharp_make_geom_info (int nrings, const int *nph, const ptrdiff_t *ofs,
Martin Reinecke's avatar
Martin Reinecke committed
147
  const int *stride, const double *phi0, const double *theta,
Martin Reinecke's avatar
Martin Reinecke committed
148
  const double *wgt);
Martin Reinecke's avatar
Martin Reinecke committed
149
150
151
152
153
154
155

/*! \} */

/*! \defgroup lowlevelgroup Low-level libsharp2 SHT interface */
/*! \{ */

/*! Enumeration of SHARP job types. */
Martin Reinecke's avatar
Martin Reinecke committed
156
enum sharp_jobtype { SHARP_YtW=0,               /*!< analysis */
Martin Reinecke's avatar
Martin Reinecke committed
157
158
159
160
161
162
               SHARP_MAP2ALM=SHARP_YtW,   /*!< analysis */
               SHARP_Y=1,                 /*!< synthesis */
               SHARP_ALM2MAP=SHARP_Y,     /*!< synthesis */
               SHARP_Yt=2,                /*!< adjoint synthesis */
               SHARP_WY=3,                /*!< adjoint analysis */
               SHARP_ALM2MAP_DERIV1=4     /*!< synthesis of first derivatives */
Martin Reinecke's avatar
Martin Reinecke committed
163
             };
Martin Reinecke's avatar
Martin Reinecke committed
164
165

/*! Job flags */
Martin Reinecke's avatar
Martin Reinecke committed
166
enum sharp_jobflags { SHARP_DP              = 1<<4,
Martin Reinecke's avatar
Martin Reinecke committed
167
168
169
170
171
172
173
174
175
176
177
               /*!< map and a_lm are in double precision */
               SHARP_ADD             = 1<<5,
               /*!< results are added to the output arrays, instead of
                    overwriting them */

               /* NOTE: SHARP_REAL_HARMONICS, 1<<6, is also available in sharp_jobflags,
                  but its use here is deprecated in favor of having it in the sharp_alm_info */

               SHARP_NO_FFT          = 1<<7,

               SHARP_USE_WEIGHTS     = 1<<20,    /* internal use only */
Martin Reinecke's avatar
Martin Reinecke committed
178
             };
Martin Reinecke's avatar
Martin Reinecke committed
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205

/*! Performs a libsharp2 SHT job. The interface deliberately does not use
  the C99 "complex" data type, in order to be callable from C89 and C++.
  \param type the type of SHT
  \param spin the spin of the quantities to be transformed
  \param alm contains pointers to the a_lm coefficients. If \a spin==0,
    alm[0] points to the a_lm of the SHT. If \a spin>0, alm[0] and alm[1]
    point to the two a_lm sets of the SHT. The exact data type of \a alm
    depends on whether the SHARP_DP flag is set.
  \param map contains pointers to the maps. If \a spin==0,
    map[0] points to the map of the SHT. If \a spin>0, or \a type is
    SHARP_ALM2MAP_DERIV1, map[0] and map[1] point to the two maps of the SHT.
    The exact data type of \a map depends on whether the SHARP_DP flag is set.
  \param geom_info A \c sharp_geom_info object compatible with the provided
    \a map arrays.
  \param alm_info A \c sharp_alm_info object compatible with the provided
    \a alm arrays. All \c m values from 0 to some \c mmax<=lmax must be present
    exactly once.
  \param flags See sharp_jobflags. In particular, if SHARP_DP is set, then
    \a alm is expected to have the type "complex double **" and \a map is
    expected to have the type "double **"; otherwise, the expected
    types are "complex float **" and "float **", respectively.
  \param time If not NULL, the wall clock time required for this SHT
    (in seconds) will be written here.
  \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,
Martin Reinecke's avatar
Martin Reinecke committed
206
  const sharp_geom_info &geom_info, const sharp_alm_info &alm_info,
Martin Reinecke's avatar
Martin Reinecke committed
207
208
209
210
211
212
213
214
215
216
217
218
  int flags, double *time, unsigned long long *opcnt);

void sharp_set_chunksize_min(int new_chunksize_min);
void sharp_set_nchunks_max(int new_nchunks_max);

/*! \} */

int sharp_get_mlim (int lmax, int spin, double sth, double cth);
int sharp_veclen(void);
const char *sharp_architecture(void);

#endif