sharp2_testsuite.cc 16.5 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
/*
 *  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_testsuite.c
 *
23
 *  Copyright (C) 2012-2020 Max-Planck-Society
Martin Reinecke's avatar
Martin Reinecke committed
24
25
26
27
28
29
30
31
32
33
34
 *  \author Martin Reinecke
 */

#include <iostream>
#include <complex>
#include "libsharp2/sharp.h"
#include "libsharp2/sharp_geomhelpers.h"
#include "libsharp2/sharp_almhelpers.h"
#include "mr_util/system.h"
#include "mr_util/error_handling.h"
#include "mr_util/threading.h"
Martin Reinecke's avatar
Martin Reinecke committed
35
#include "mr_util/math_utils.h"
36
#include "mr_util/string_utils.h"
37
#include "mr_util/gl_integrator.h"
Martin Reinecke's avatar
Martin Reinecke committed
38
39

using namespace std;
Martin Reinecke's avatar
Martin Reinecke committed
40
using namespace mr;
Martin Reinecke's avatar
Martin Reinecke committed
41
42
43

static void threading_status(void)
  {
44
  cout << "Threading: " << mr::get_default_nthreads() << " threads active." << endl;
Martin Reinecke's avatar
Martin Reinecke committed
45
46
47
48
  }

static void MPI_status(void)
  {
Martin Reinecke's avatar
Martin Reinecke committed
49
  cout << "MPI: not supported by this binary" << endl;
Martin Reinecke's avatar
Martin Reinecke committed
50
51
  }

Martin Reinecke's avatar
Martin Reinecke committed
52
static void sharp_announce (const string &name)
Martin Reinecke's avatar
Martin Reinecke committed
53
  {
Martin Reinecke's avatar
Martin Reinecke committed
54
  size_t m, nlen=name.length();
Martin Reinecke's avatar
Martin Reinecke committed
55
56
57
58
59
60
61
62
63
  cout << "\n+-";
  for (m=0; m<nlen; ++m) cout << "-";
  cout << "-+\n";
  cout << "| " << name << " |\n";
  cout << "+-";
  for (m=0; m<nlen; ++m) cout << "-";
  cout << "-+\n\n";
  cout << "Detected hardware architecture: " << sharp_architecture() << endl;
  cout << "Supported vector length: " << sharp_veclen() << endl;
Martin Reinecke's avatar
Martin Reinecke committed
64
65
  threading_status();
  MPI_status();
Martin Reinecke's avatar
Martin Reinecke committed
66
  cout << endl;
Martin Reinecke's avatar
Martin Reinecke committed
67
68
  }

Martin Reinecke's avatar
Martin Reinecke committed
69
70
static void sharp_module_startup (const string &name, int argc, int argc_expected,
  const string &argv_expected, int verbose)
Martin Reinecke's avatar
Martin Reinecke committed
71
72
73
  {
  if (verbose) sharp_announce (name);
  if (argc==argc_expected) return;
Martin Reinecke's avatar
Martin Reinecke committed
74
  if (verbose) cerr << "Usage: " << name << " " << argv_expected << endl;
Martin Reinecke's avatar
Martin Reinecke committed
75
76
77
  exit(1);
  }

Martin Reinecke's avatar
more    
Martin Reinecke committed
78
using dcmplx = complex<double>;
Martin Reinecke's avatar
Martin Reinecke committed
79
80
81
82
83
84
85
86
87

int ntasks, mytask;

static double drand (double min, double max, unsigned *state)
  {
  *state = (((*state) * 1103515245u) + 12345u) & 0x7fffffffu;
  return min + (max-min)*(*state)/(0x7fffffff+1.0);
  }

88
static void random_alm (dcmplx *alm, sharp_standard_alm_info &helper, size_t spin, size_t cnt)
Martin Reinecke's avatar
Martin Reinecke committed
89
  {
90
// FIXME: re-introduce multi-threading?
91
  for (size_t mi=0;mi<helper.nm(); ++mi)
Martin Reinecke's avatar
Martin Reinecke committed
92
    {
93
    auto m=helper.mval(mi);
94
    unsigned state=1234567u*unsigned(cnt)+8912u*unsigned(m); // random seed
95
    for (auto l=m;l<=helper.lmax(); ++l)
Martin Reinecke's avatar
Martin Reinecke committed
96
97
      {
      if ((l<spin)&&(m<spin))
Martin Reinecke's avatar
Martin Reinecke committed
98
        alm[helper.index(l,mi)] = 0.;
Martin Reinecke's avatar
Martin Reinecke committed
99
100
101
102
      else
        {
        double rv = drand(-1,1,&state);
        double iv = (m==0) ? 0 : drand(-1,1,&state);
Martin Reinecke's avatar
Martin Reinecke committed
103
        alm[helper.index(l,mi)] = dcmplx(rv,iv);
Martin Reinecke's avatar
Martin Reinecke committed
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
        }
      }
    }
  }

static unsigned long long totalops (unsigned long long val)
  {
  return val;
  }

static double maxTime (double val)
  {
  return val;
  }

static double allreduceSumDouble (double val)
  {
  return val;
  }

static double totalMem()
  {
  return mr::getProcessInfo("VmHWM")*1024;
  }

129
static size_t get_nalms(const sharp_alm_info &ainfo)
Martin Reinecke's avatar
Martin Reinecke committed
130
  {
131
  size_t res=0;
132
133
  for (size_t i=0; i<ainfo.nm(); ++i)
    res += ainfo.lmax()-ainfo.mval(i)+1;
Martin Reinecke's avatar
Martin Reinecke committed
134
135
136
  return res;
  }

137
static size_t get_npix(const sharp_geom_info &ginfo)
Martin Reinecke's avatar
Martin Reinecke committed
138
  {
139
140
141
  size_t res=0;
  for (size_t i=0; i<ginfo.nrings(); ++i)
    res += ginfo.nph(i);
Martin Reinecke's avatar
Martin Reinecke committed
142
143
144
  return res;
  }

145
static vector<double> get_sqsum_and_invert (dcmplx **alm, size_t nalms, size_t ncomp)
Martin Reinecke's avatar
Martin Reinecke committed
146
  {
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
147
  vector<double> sqsum(ncomp);
148
  for (size_t i=0; i<ncomp; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
149
150
    {
    sqsum[i]=0;
151
    for (size_t j=0; j<nalms; ++j)
Martin Reinecke's avatar
Martin Reinecke committed
152
153
154
155
156
157
158
159
      {
      sqsum[i]+=norm(alm[i][j]);
      alm[i][j]=-alm[i][j];
      }
    }
  return sqsum;
  }

160
static void get_errors (dcmplx **alm, size_t nalms, size_t ncomp, const vector<double> &sqsum,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
161
  vector<double> &err_abs, vector<double> &err_rel)
Martin Reinecke's avatar
Martin Reinecke committed
162
  {
163
  size_t nalms_tot=nalms;
Martin Reinecke's avatar
Martin Reinecke committed
164

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
165
166
  err_abs.resize(ncomp);
  err_rel.resize(ncomp);
167
  for (size_t i=0; i<ncomp; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
168
169
    {
    double sum=0, maxdiff=0, sumtot, sqsumtot, maxdifftot;
170
    for (size_t j=0; j<nalms; ++j)
Martin Reinecke's avatar
Martin Reinecke committed
171
172
173
174
175
      {
      double sqr=norm(alm[i][j]);
      sum+=sqr;
      if (sqr>maxdiff) maxdiff=sqr;
      }
176
    maxdiff=sqrt(maxdiff);
Martin Reinecke's avatar
Martin Reinecke committed
177
178
179
180
181
182

    sumtot=sum;
    sqsumtot=sqsum[i];
    maxdifftot=maxdiff;
    sumtot=sqrt(sumtot/nalms_tot);
    sqsumtot=sqrt(sqsumtot/nalms_tot);
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
183
184
    err_abs[i]=maxdifftot;
    err_rel[i]=sumtot/sqsumtot;
Martin Reinecke's avatar
Martin Reinecke committed
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
    }
  }

static int good_fft_size(int n)
  {
  if (n<=6) return n;
  int bestfac=2*n;

  for (int f2=1; f2<bestfac; f2*=2)
    for (int f23=f2; f23<bestfac; f23*=3)
      for (int f235=f23; f235<bestfac; f235*=5)
        if (f235>=n) bestfac=f235;

  return bestfac;
  }

Martin Reinecke's avatar
Martin Reinecke committed
201
static void get_infos (const string &gname, int lmax, int &mmax, int &gpar1,
202
  int &gpar2, unique_ptr<sharp_geom_info> &ginfo, unique_ptr<sharp_standard_alm_info> &ainfo, int verbose)
Martin Reinecke's avatar
Martin Reinecke committed
203
204
  {
  MR_assert(lmax>=0,"lmax must not be negative");
Martin Reinecke's avatar
Martin Reinecke committed
205
206
  if (mmax<0) mmax=lmax;
  MR_assert(mmax<=lmax,"mmax larger than lmax");
Martin Reinecke's avatar
Martin Reinecke committed
207
208

  verbose &= (mytask==0);
Martin Reinecke's avatar
Martin Reinecke committed
209
  if (verbose) cout << "lmax: " << lmax << ", mmax: " << mmax << endl;
Martin Reinecke's avatar
Martin Reinecke committed
210

Martin Reinecke's avatar
Martin Reinecke committed
211
  ainfo=sharp_make_triangular_alm_info(lmax,mmax,1);
Martin Reinecke's avatar
Martin Reinecke committed
212

Martin Reinecke's avatar
Martin Reinecke committed
213
  if (gname=="healpix")
Martin Reinecke's avatar
Martin Reinecke committed
214
    {
Martin Reinecke's avatar
Martin Reinecke committed
215
216
217
218
    if (gpar1<1) gpar1=lmax/2;
    if (gpar1==0) ++gpar1;
    ginfo=sharp_make_healpix_geom_info (gpar1, 1);
    if (verbose) cout << "HEALPix grid, nside=" << gpar1 << endl;
Martin Reinecke's avatar
Martin Reinecke committed
219
    }
Martin Reinecke's avatar
Martin Reinecke committed
220
  else if (gname=="gauss")
Martin Reinecke's avatar
Martin Reinecke committed
221
    {
Martin Reinecke's avatar
Martin Reinecke committed
222
223
224
    if (gpar1<1) gpar1=lmax+1;
    if (gpar2<1) gpar2=2*mmax+1;
    ginfo=sharp_make_gauss_geom_info (gpar1, gpar2, 0., 1, gpar2);
Martin Reinecke's avatar
Martin Reinecke committed
225
    if (verbose)
Martin Reinecke's avatar
Martin Reinecke committed
226
      cout << "Gauss-Legendre grid, nlat=" << gpar1 << ", nlon=" << gpar2 << endl;
Martin Reinecke's avatar
Martin Reinecke committed
227
    }
Martin Reinecke's avatar
Martin Reinecke committed
228
  else if (gname=="fejer1")
Martin Reinecke's avatar
Martin Reinecke committed
229
    {
Martin Reinecke's avatar
Martin Reinecke committed
230
231
232
    if (gpar1<1) gpar1=2*lmax+1;
    if (gpar2<1) gpar2=2*mmax+1;
    ginfo=sharp_make_fejer1_geom_info (gpar1, gpar2, 0., 1, gpar2);
Martin Reinecke's avatar
Martin Reinecke committed
233
    if (verbose)
Martin Reinecke's avatar
Martin Reinecke committed
234
      cout << "Fejer1 grid, nlat=" << gpar1 << ", nlon=" << gpar2 << endl;
Martin Reinecke's avatar
Martin Reinecke committed
235
    }
Martin Reinecke's avatar
Martin Reinecke committed
236
  else if (gname=="fejer2")
Martin Reinecke's avatar
Martin Reinecke committed
237
    {
Martin Reinecke's avatar
Martin Reinecke committed
238
239
240
    if (gpar1<1) gpar1=2*lmax+1;
    if (gpar2<1) gpar2=2*mmax+1;
    ginfo=sharp_make_fejer2_geom_info (gpar1, gpar2, 0., 1, gpar2);
Martin Reinecke's avatar
Martin Reinecke committed
241
    if (verbose)
Martin Reinecke's avatar
Martin Reinecke committed
242
      cout << "Fejer2 grid, nlat=" << gpar1 << ", nlon=" << gpar2 << endl;
Martin Reinecke's avatar
Martin Reinecke committed
243
    }
244
245
246
247
248
249
250
251
  else if (gname=="dh")
    {
    if (gpar1<1) gpar1=2*lmax+2;
    if (gpar2<1) gpar2=2*mmax+1;
    ginfo=sharp_make_dh_geom_info (gpar1, gpar2, 0., 1, gpar2);
    if (verbose)
      cout << "Driscoll-Healy grid, nlat=" << gpar1 << ", nlon=" << gpar2 << endl;
    }
Martin Reinecke's avatar
Martin Reinecke committed
252
  else if (gname=="cc")
Martin Reinecke's avatar
Martin Reinecke committed
253
    {
Martin Reinecke's avatar
Martin Reinecke committed
254
255
256
    if (gpar1<1) gpar1=2*lmax+1;
    if (gpar2<1) gpar2=2*mmax+1;
    ginfo=sharp_make_cc_geom_info (gpar1, gpar2, 0., 1, gpar2);
Martin Reinecke's avatar
Martin Reinecke committed
257
    if (verbose)
Martin Reinecke's avatar
Martin Reinecke committed
258
      cout << "Clenshaw-Curtis grid, nlat=" << gpar1 << ", nlon=" << gpar2 << endl;
Martin Reinecke's avatar
Martin Reinecke committed
259
    }
Martin Reinecke's avatar
Martin Reinecke committed
260
  else if (gname=="smallgauss")
Martin Reinecke's avatar
Martin Reinecke committed
261
    {
Martin Reinecke's avatar
Martin Reinecke committed
262
    int nlat=gpar1, nlon=gpar2;
Martin Reinecke's avatar
Martin Reinecke committed
263
    if (nlat<1) nlat=lmax+1;
Martin Reinecke's avatar
Martin Reinecke committed
264
    if (nlon<1) nlon=2*mmax+1;
265
    size_t npix_o = nlat*nlon;
Martin Reinecke's avatar
Martin Reinecke committed
266
    gpar1=nlat; gpar2=nlon;
267
268
269
270
271
272
273
274
275
276
277
    const double pi=3.141592653589793238462643383279502884197;

    vector<size_t> nph(nlat);
    vector<double> phi0_(nlat);
    vector<ptrdiff_t> ofs(nlat);

    GL_Integrator integ(nlat);
    auto theta = integ.coords();
    auto weight = integ.weights();
    size_t ofs_ = 0;
    for (int m=0; m<nlat; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
278
      {
279
280
281
282
283
284
      theta[m] = acos(-theta[m]);
      nph[m] = good_fft_size(min<int>(nlon, 1+2*sharp_get_mlim(lmax,0,sin(theta[m]),cos(theta[m]))));
      phi0_[m] = 0.;
      ofs[m]=ofs_;
      ofs_+=nph[m];
      weight[m]*=2*pi/nph[m];
Martin Reinecke's avatar
Martin Reinecke committed
285
      }
286
287

    ginfo = unique_ptr<sharp_geom_info>(new sharp_standard_geom_info(nlat, nph.data(), ofs.data(), 1, phi0_.data(), theta.data(), weight.data()));
Martin Reinecke's avatar
Martin Reinecke committed
288
289
    if (verbose)
      {
290
      auto npix=get_npix(*ginfo);
Martin Reinecke's avatar
Martin Reinecke committed
291
      cout << "Small Gauss grid, nlat=" << nlat << ", npix=" << npix
292
           << ", savings=" << ((npix_o-npix)*100./npix_o) << "%" << endl;
Martin Reinecke's avatar
Martin Reinecke committed
293
294
295
296
297
298
299
300
      }
    }
  else
    MR_fail("unknown grid geometry");
  }

static void check_sign_scale(void)
  {
301
302
303
304
305
  size_t lmax=50;
  auto mmax=lmax;
  auto nrings=lmax+1;
  auto ppring=2*lmax+2;
  auto npix=ptrdiff_t(nrings*ppring);
Martin Reinecke's avatar
Martin Reinecke committed
306
  auto tinfo = sharp_make_gauss_geom_info (nrings, ppring, 0., 1, ppring);
Martin Reinecke's avatar
Martin Reinecke committed
307

Martin Reinecke's avatar
Martin Reinecke committed
308
  auto alms = sharp_make_triangular_alm_info(lmax,mmax,1);
Martin Reinecke's avatar
Martin Reinecke committed
309
310
  ptrdiff_t nalms = ((mmax+1)*(mmax+2))/2 + (mmax+1)*(lmax-mmax);

Martin Reinecke's avatar
Martin Reinecke committed
311
  vector<double> bmap(2*npix);
Martin Reinecke's avatar
Martin Reinecke committed
312
313
  vector<double *>map1({&bmap[0]});
  vector<double *>map2({&bmap[0], &bmap[npix]});
Martin Reinecke's avatar
Martin Reinecke committed
314

Martin Reinecke's avatar
Martin Reinecke committed
315
  vector<dcmplx> balm(2*nalms,dcmplx(1.,1.));
Martin Reinecke's avatar
Martin Reinecke committed
316
317
  vector<dcmplx *>alm1({&balm[0]});
  vector<dcmplx *>alm2({&balm[0], &balm[nalms]});
Martin Reinecke's avatar
Martin Reinecke committed
318

Martin Reinecke's avatar
Martin Reinecke committed
319
320
321
  /* use mirrored indices to emulate the "old" Gaussian grid geometry */
  /* original indices were 0, npix/2 and npix-1 */
  size_t i0=npix-ppring, i1=npix/2, i2=ppring-1;
Martin Reinecke's avatar
Martin Reinecke committed
322
  sharp_alm2map(&balm[0],&bmap[0],*tinfo,*alms,0,nullptr,nullptr);
Martin Reinecke's avatar
Martin Reinecke committed
323
  MR_assert(approx(map1[0][i0], 3.588246976618616912e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
324
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
325
  MR_assert(approx(map1[0][i1], 4.042209792157496651e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
326
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
327
  MR_assert(approx(map1[0][i2],-1.234675107554816442e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
328
329
    "error");

Martin Reinecke's avatar
Martin Reinecke committed
330
  sharp_alm2map_spin(1, &balm[0], &balm[nalms], &bmap[0], &bmap[npix],*tinfo,*alms,0,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
331
    nullptr,nullptr);
Martin Reinecke's avatar
Martin Reinecke committed
332
  MR_assert(approx(map2[0][i0], 2.750897760535633285e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
333
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
334
  MR_assert(approx(map2[0][i1], 3.137704477368562905e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
335
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
336
  MR_assert(approx(map2[0][i2],-8.405730859837063917e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
337
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
338
  MR_assert(approx(map2[1][i0],-2.398026536095463346e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
339
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
340
  MR_assert(approx(map2[1][i1],-4.961140548331700728e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
341
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
342
  MR_assert(approx(map2[1][i2],-1.412765834230440021e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
343
344
    "error");

Martin Reinecke's avatar
Martin Reinecke committed
345
  sharp_alm2map_spin(2,&balm[0], &balm[nalms],&bmap[0], &bmap[npix],*tinfo,*alms,0,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
346
    nullptr,nullptr);
Martin Reinecke's avatar
Martin Reinecke committed
347
  MR_assert(approx(map2[0][i0],-1.398186224727334448e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
348
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
349
  MR_assert(approx(map2[0][i1],-2.456676000884031197e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
350
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
351
  MR_assert(approx(map2[0][i2],-1.516249174408820863e+02,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
352
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
353
  MR_assert(approx(map2[1][i0],-3.173406200299964119e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
354
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
355
  MR_assert(approx(map2[1][i1],-5.831327404513146462e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
356
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
357
  MR_assert(approx(map2[1][i2],-1.863257892248353897e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
358
359
    "error");

Martin Reinecke's avatar
Martin Reinecke committed
360
  sharp_execute(SHARP_ALM2MAP_DERIV1,1,{&balm[0]},{&bmap[0], &bmap[npix]},*tinfo,*alms,
Martin Reinecke's avatar
Martin Reinecke committed
361
362
    0,nullptr,nullptr);
  MR_assert(approx(map2[0][i0],-6.859393905369091105e-01,1e-11),
Martin Reinecke's avatar
Martin Reinecke committed
363
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
364
  MR_assert(approx(map2[0][i1],-2.103947835973212364e+02,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
365
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
366
  MR_assert(approx(map2[0][i2],-1.092463246472086439e+03,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
367
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
368
  MR_assert(approx(map2[1][i0],-1.411433220713928165e+02,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
369
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
370
  MR_assert(approx(map2[1][i1],-1.146122859381925082e+03,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
371
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
372
  MR_assert(approx(map2[1][i2], 7.821618677689795049e+02,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
373
374
375
    "error");
  }

376
static void do_sht (sharp_geom_info &ginfo, sharp_standard_alm_info &ainfo,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
377
  int spin, vector<double> &err_abs, vector<double> &err_rel,
Martin Reinecke's avatar
Martin Reinecke committed
378
379
380
381
  double *t_a2m, double *t_m2a, unsigned long long *op_a2m,
  unsigned long long *op_m2a, size_t ntrans)
  {
  ptrdiff_t nalms = get_nalms(ainfo);
Martin Reinecke's avatar
Martin Reinecke committed
382
  size_t ncomp = (spin==0) ? 1 : 2;
Martin Reinecke's avatar
Martin Reinecke committed
383
384

  size_t npix = get_npix(ginfo);
Martin Reinecke's avatar
Martin Reinecke committed
385
386
  vector<double> bmap(ncomp*ntrans*npix, 0.);
  vector<double *>map(ncomp*ntrans);
387
  for (size_t i=0; i<ncomp*ntrans; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
388
    map[i]=&bmap[i*npix];
Martin Reinecke's avatar
Martin Reinecke committed
389

Martin Reinecke's avatar
Martin Reinecke committed
390
391
  vector<dcmplx> balm(ncomp*ntrans*nalms);
  vector<dcmplx *>alm(ncomp*ntrans);
392
  for (size_t i=0; i<ncomp*ntrans; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
393
394
    {
    alm[i] = &balm[i*nalms];
Martin Reinecke's avatar
Martin Reinecke committed
395
    random_alm(alm[i],ainfo,spin,i+1);
Martin Reinecke's avatar
Martin Reinecke committed
396
    }
Martin Reinecke's avatar
Martin Reinecke committed
397
398
399
400

  double tta2m, ttm2a;
  unsigned long long toa2m, tom2a;

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
401
402
  if (t_a2m!=nullptr) *t_a2m=0;
  if (op_a2m!=nullptr) *op_a2m=0;
Martin Reinecke's avatar
Martin Reinecke committed
403
404
  for (size_t itrans=0; itrans<ntrans; ++itrans)
    {
Martin Reinecke's avatar
Martin Reinecke committed
405
    vector<any> av, mv;
Martin Reinecke's avatar
Martin Reinecke committed
406
407
408
409
410
411
    for (size_t i=0; i<ncomp; ++i)
      {
      av.push_back(alm[itrans*ncomp+i]);
      mv.push_back(map[itrans*ncomp+i]);
      }
    sharp_execute(SHARP_ALM2MAP,spin,av,mv,ginfo,ainfo, 0,&tta2m,&toa2m);
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
412
413
    if (t_a2m!=nullptr) *t_a2m+=maxTime(tta2m);
    if (op_a2m!=nullptr) *op_a2m+=totalops(toa2m);
Martin Reinecke's avatar
Martin Reinecke committed
414
    }
Martin Reinecke's avatar
Martin Reinecke committed
415
  auto sqsum=get_sqsum_and_invert(alm.data(),nalms,ntrans*ncomp);
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
416
417
  if (t_m2a!=nullptr) *t_m2a=0;
  if (op_m2a!=nullptr) *op_m2a=0;
Martin Reinecke's avatar
Martin Reinecke committed
418
419
  for (size_t itrans=0; itrans<ntrans; ++itrans)
    {
Martin Reinecke's avatar
Martin Reinecke committed
420
    vector<any> av, mv;
Martin Reinecke's avatar
Martin Reinecke committed
421
422
423
424
425
426
427
    for (size_t i=0; i<ncomp; ++i)
      {
      av.push_back(alm[itrans*ncomp+i]);
      mv.push_back(map[itrans*ncomp+i]);
      }
    sharp_execute(SHARP_MAP2ALM,spin,av,mv,ginfo,ainfo,
      SHARP_ADD,&ttm2a,&tom2a);
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
428
429
    if (t_m2a!=nullptr) *t_m2a+=maxTime(ttm2a);
    if (op_m2a!=nullptr) *op_m2a+=totalops(tom2a);
Martin Reinecke's avatar
Martin Reinecke committed
430
    }
Martin Reinecke's avatar
Martin Reinecke committed
431
  get_errors(alm.data(), nalms, ntrans*ncomp, sqsum, err_abs, err_rel);
Martin Reinecke's avatar
Martin Reinecke committed
432
433
  }

434
static void check_accuracy (sharp_geom_info &ginfo, sharp_standard_alm_info &ainfo,
Martin Reinecke's avatar
Martin Reinecke committed
435
436
437
  int spin)
  {
  int ncomp = (spin==0) ? 1 : 2;
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
438
  vector<double> err_abs, err_rel;
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
439
440
  do_sht (ginfo, ainfo, spin, err_abs, err_rel, nullptr, nullptr,
    nullptr, nullptr, 1);
Martin Reinecke's avatar
Martin Reinecke committed
441
442
443
444
445
446
  for (int i=0; i<ncomp; ++i)
    MR_assert((err_rel[i]<1e-10) && (err_abs[i]<1e-10),"error");
  }

static void run(int lmax, int mmax, int nlat, int nlon, int spin)
  {
Martin Reinecke's avatar
Martin Reinecke committed
447
  unique_ptr<sharp_geom_info> ginfo;
448
  unique_ptr<sharp_standard_alm_info> ainfo;
Martin Reinecke's avatar
Martin Reinecke committed
449
450
  get_infos ("gauss", lmax, mmax, nlat, nlon, ginfo, ainfo, 0);
  check_accuracy(*ginfo,*ainfo,spin);
Martin Reinecke's avatar
Martin Reinecke committed
451
452
453
454
455
456
  }

static void sharp_acctest(void)
  {
  if (mytask==0) sharp_module_startup("sharp_acctest",1,1,"",1);

Martin Reinecke's avatar
Martin Reinecke committed
457
  if (mytask==0) cout << "Checking signs and scales.\n";
Martin Reinecke's avatar
Martin Reinecke committed
458
  check_sign_scale();
Martin Reinecke's avatar
Martin Reinecke committed
459
  if (mytask==0) cout << "Passed.\n\n";
Martin Reinecke's avatar
Martin Reinecke committed
460

Martin Reinecke's avatar
Martin Reinecke committed
461
  if (mytask==0) cout << "Testing map analysis accuracy.\n";
Martin Reinecke's avatar
Martin Reinecke committed
462
463
464
465
466
467
468
469
470
471

  run(127, 127, 128, 256, 0);
  run(127, 127, 128, 256, 1);
  run(127, 127, 128, 256, 2);
  run(127, 127, 128, 256, 3);
  run(127, 127, 128, 256, 30);
  run(5, 0, 6, 1, 0);
  run(5, 0, 7, 2, 0);
  run(8, 8, 9, 17, 0);
  run(8, 8, 9, 17, 2);
Martin Reinecke's avatar
Martin Reinecke committed
472
  if (mytask==0) cout << "Passed.\n\n";
Martin Reinecke's avatar
Martin Reinecke committed
473
474
475
476
477
478
  }

static void sharp_test (int argc, const char **argv)
  {
  if (mytask==0) sharp_announce("sharp_test");
  MR_assert(argc>=8,"usage: grid lmax mmax geom1 geom2 spin [ntrans]");
Martin Reinecke's avatar
Martin Reinecke committed
479
480
481
482
483
  auto lmax=stringToData<int>(argv[3]);
  auto mmax=stringToData<int>(argv[4]);
  auto gpar1=stringToData<int>(argv[5]);
  auto gpar2=stringToData<int>(argv[6]);
  auto spin=stringToData<int>(argv[7]);
Martin Reinecke's avatar
Martin Reinecke committed
484
  int ntrans=1;
Martin Reinecke's avatar
Martin Reinecke committed
485
  if (argc>=9) ntrans=stringToData<int>(argv[8]);
Martin Reinecke's avatar
Martin Reinecke committed
486

Martin Reinecke's avatar
Martin Reinecke committed
487
488
489
  if (mytask==0) cout << "Testing map analysis accuracy.\n";
  if (mytask==0) cout << "spin=" << spin << endl;
  if (mytask==0) cout << "ntrans=" << ntrans << endl;
Martin Reinecke's avatar
Martin Reinecke committed
490

Martin Reinecke's avatar
Martin Reinecke committed
491
  unique_ptr<sharp_geom_info> ginfo;
492
  unique_ptr<sharp_standard_alm_info> ainfo;
Martin Reinecke's avatar
Martin Reinecke committed
493
  get_infos (argv[2], lmax, mmax, gpar1, gpar2, ginfo, ainfo, 1);
Martin Reinecke's avatar
Martin Reinecke committed
494
495
496
497

  int ncomp = (spin==0) ? 1 : 2;
  double t_a2m=1e30, t_m2a=1e30;
  unsigned long long op_a2m, op_m2a;
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
498
  vector<double> err_abs, err_rel;
Martin Reinecke's avatar
Martin Reinecke committed
499
500
501
502
503
504
505

  double t_acc=0;
  int nrpt=0;
  while(1)
    {
    ++nrpt;
    double ta2m2, tm2a2;
Martin Reinecke's avatar
Martin Reinecke committed
506
    do_sht (*ginfo, *ainfo, spin, err_abs, err_rel, &ta2m2, &tm2a2,
Martin Reinecke's avatar
Martin Reinecke committed
507
508
509
510
511
512
      &op_a2m, &op_m2a, ntrans);
    if (ta2m2<t_a2m) t_a2m=ta2m2;
    if (tm2a2<t_m2a) t_m2a=tm2a2;
    t_acc+=t_a2m+t_m2a;
    if (t_acc>2.)
      {
Martin Reinecke's avatar
Martin Reinecke committed
513
      if (mytask==0) cout << "Best of " << nrpt << " runs\n";
Martin Reinecke's avatar
Martin Reinecke committed
514
515
516
517
      break;
      }
    }

Martin Reinecke's avatar
Martin Reinecke committed
518
519
520
521
  if (mytask==0) cout << "wall time for alm2map: " << t_a2m << "s\n";
  if (mytask==0) cout << "Performance: " << 1e-9*op_a2m/t_a2m << "GFLOPs/s\n";
  if (mytask==0) cout << "wall time for map2alm: " << t_m2a << "s\n";
  if (mytask==0) cout << "Performance: " << 1e-9*op_m2a/t_m2a << "GFLOPs/s\n";
Martin Reinecke's avatar
Martin Reinecke committed
522
523
524

  if (mytask==0)
    for (int i=0; i<ntrans*ncomp; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
525
      cout << "component " << i << ": rms " << err_rel[i] << ", maxerr " << err_abs[i] << endl;
Martin Reinecke's avatar
Martin Reinecke committed
526

Martin Reinecke's avatar
Martin Reinecke committed
527
  double iosize = ntrans*ncomp*(16.*get_nalms(*ainfo) + 8.*get_npix(*ginfo));
Martin Reinecke's avatar
Martin Reinecke committed
528
529
530
531
  iosize = allreduceSumDouble(iosize);

  double tmem=totalMem();
  if (mytask==0)
Martin Reinecke's avatar
Martin Reinecke committed
532
    cout << "\nMemory high water mark: " << tmem/(1<<20) << " MB\n";
Martin Reinecke's avatar
Martin Reinecke committed
533
  if (mytask==0)
Martin Reinecke's avatar
Martin Reinecke committed
534
    cout << "Memory overhead: " << (tmem-iosize)/(1<<20) << " MB ("
535
         << 100.*(1.-iosize/tmem) << "% of working set)\n";
Martin Reinecke's avatar
Martin Reinecke committed
536
537
538
539
540
541
542
543
544
545
546

  double maxerel=0., maxeabs=0.;
  for (int i=0; i<ncomp; ++i)
    {
    if (maxerel<err_rel[i]) maxerel=err_rel[i];
    if (maxeabs<err_abs[i]) maxeabs=err_abs[i];
    }
  }

int main(int argc, const char **argv)
  {
547
548
549
550
551
552
553
  const char *tmp = getenv("OMP_NUM_THREADS");
  if (tmp!=nullptr)
    {
    string tmp2(tmp);
    if (trim(tmp2)!="")
      mr::set_default_nthreads(stringToData<size_t>(tmp));
    }
Martin Reinecke's avatar
Martin Reinecke committed
554
555
556
  mytask=0; ntasks=1;

  MR_assert(argc>=2,"need at least one command line argument");
Martin Reinecke's avatar
Martin Reinecke committed
557
  auto mode = string(argv[1]);
Martin Reinecke's avatar
Martin Reinecke committed
558

Martin Reinecke's avatar
Martin Reinecke committed
559
  if (mode=="acctest")
Martin Reinecke's avatar
Martin Reinecke committed
560
    sharp_acctest();
Martin Reinecke's avatar
Martin Reinecke committed
561
  else if (mode=="test")
Martin Reinecke's avatar
Martin Reinecke committed
562
563
564
565
566
567
    sharp_test(argc,argv);
  else
    MR_fail("unknown command");

  return 0;
  }