sharp2_testsuite.cc 16.1 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
35
 *  \author Martin Reinecke
 */

#include <iostream>
#include <complex>
using std::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
36
#include "mr_util/math_utils.h"
37
#include "mr_util/string_utils.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
78
79
80
81
82
83
84
85
86
87
  exit(1);
  }

typedef complex<double> dcmplx;

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);
  }

Martin Reinecke's avatar
Martin Reinecke committed
88
static void random_alm (dcmplx *alm, sharp_alm_info &helper, int spin, int cnt)
Martin Reinecke's avatar
Martin Reinecke committed
89
90
91
92
93
  {
#pragma omp parallel
{
  int mi;
#pragma omp for schedule (dynamic,100)
Martin Reinecke's avatar
Martin Reinecke committed
94
  for (mi=0;mi<helper.nm; ++mi)
Martin Reinecke's avatar
Martin Reinecke committed
95
    {
Martin Reinecke's avatar
Martin Reinecke committed
96
    int m=helper.mval[mi];
Martin Reinecke's avatar
Martin Reinecke committed
97
    unsigned state=1234567u*(unsigned)cnt+8912u*(unsigned)m; // random seed
Martin Reinecke's avatar
Martin Reinecke committed
98
    for (int l=m;l<=helper.lmax; ++l)
Martin Reinecke's avatar
Martin Reinecke committed
99
100
      {
      if ((l<spin)&&(m<spin))
Martin Reinecke's avatar
Martin Reinecke committed
101
        alm[helper.index(l,mi)] = 0.;
Martin Reinecke's avatar
Martin Reinecke committed
102
103
104
105
      else
        {
        double rv = drand(-1,1,&state);
        double iv = (m==0) ? 0 : drand(-1,1,&state);
Martin Reinecke's avatar
Martin Reinecke committed
106
        alm[helper.index(l,mi)] = dcmplx(rv,iv);
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
127
128
129
130
131
132
        }
      }
    }
} // end of parallel region
  }

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;
  }

Martin Reinecke's avatar
Martin Reinecke committed
133
static ptrdiff_t get_nalms(const sharp_alm_info &ainfo)
Martin Reinecke's avatar
Martin Reinecke committed
134
135
  {
  ptrdiff_t res=0;
Martin Reinecke's avatar
Martin Reinecke committed
136
137
  for (int i=0; i<ainfo.nm; ++i)
    res += ainfo.lmax-ainfo.mval[i]+1;
Martin Reinecke's avatar
Martin Reinecke committed
138
139
140
  return res;
  }

Martin Reinecke's avatar
Martin Reinecke committed
141
static ptrdiff_t get_npix(const sharp_geom_info &ginfo)
Martin Reinecke's avatar
Martin Reinecke committed
142
143
  {
  ptrdiff_t res=0;
Martin Reinecke's avatar
Martin Reinecke committed
144
  for (int i=0; i<ginfo.pair.size(); ++i)
Martin Reinecke's avatar
Martin Reinecke committed
145
    {
Martin Reinecke's avatar
Martin Reinecke committed
146
147
    res += ginfo.pair[i].r1.nph;
    if (ginfo.pair[i].r2.nph>0) res += ginfo.pair[i].r2.nph;
Martin Reinecke's avatar
Martin Reinecke committed
148
149
150
151
    }
  return res;
  }

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
152
static vector<double> get_sqsum_and_invert (dcmplx **alm, ptrdiff_t nalms, int ncomp)
Martin Reinecke's avatar
Martin Reinecke committed
153
  {
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
154
  vector<double> sqsum(ncomp);
Martin Reinecke's avatar
Martin Reinecke committed
155
156
157
158
159
160
161
162
163
164
165
166
  for (int i=0; i<ncomp; ++i)
    {
    sqsum[i]=0;
    for (ptrdiff_t j=0; j<nalms; ++j)
      {
      sqsum[i]+=norm(alm[i][j]);
      alm[i][j]=-alm[i][j];
      }
    }
  return sqsum;
  }

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
167
168
static void get_errors (dcmplx **alm, ptrdiff_t nalms, int ncomp, const vector<double> &sqsum,
  vector<double> &err_abs, vector<double> &err_rel)
Martin Reinecke's avatar
Martin Reinecke committed
169
170
171
  {
  long nalms_tot=nalms;

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
172
173
  err_abs.resize(ncomp);
  err_rel.resize(ncomp);
Martin Reinecke's avatar
Martin Reinecke committed
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
  for (int i=0; i<ncomp; ++i)
    {
    double sum=0, maxdiff=0, sumtot, sqsumtot, maxdifftot;
    for (ptrdiff_t j=0; j<nalms; ++j)
      {
      double sqr=norm(alm[i][j]);
      sum+=sqr;
      if (sqr>maxdiff) maxdiff=sqr;
      }
   maxdiff=sqrt(maxdiff);

    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
190
191
    err_abs[i]=maxdifftot;
    err_rel[i]=sumtot/sqsumtot;
Martin Reinecke's avatar
Martin Reinecke committed
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
    }
  }

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
208
209
static void get_infos (const string &gname, int lmax, int &mmax, int &gpar1,
  int &gpar2, unique_ptr<sharp_geom_info> &ginfo, unique_ptr<sharp_alm_info> &ainfo, int verbose)
Martin Reinecke's avatar
Martin Reinecke committed
210
211
  {
  MR_assert(lmax>=0,"lmax must not be negative");
Martin Reinecke's avatar
Martin Reinecke committed
212
213
  if (mmax<0) mmax=lmax;
  MR_assert(mmax<=lmax,"mmax larger than lmax");
Martin Reinecke's avatar
Martin Reinecke committed
214
215

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

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

Martin Reinecke's avatar
Martin Reinecke committed
220
  if (gname=="healpix")
Martin Reinecke's avatar
Martin Reinecke committed
221
    {
Martin Reinecke's avatar
Martin Reinecke committed
222
223
224
225
    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
226
    }
Martin Reinecke's avatar
Martin Reinecke committed
227
  else if (gname=="gauss")
Martin Reinecke's avatar
Martin Reinecke committed
228
    {
Martin Reinecke's avatar
Martin Reinecke committed
229
230
231
    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
232
    if (verbose)
Martin Reinecke's avatar
Martin Reinecke committed
233
      cout << "Gauss-Legendre grid, nlat=" << gpar1 << ", nlon=" << gpar2 << endl;
Martin Reinecke's avatar
Martin Reinecke committed
234
    }
Martin Reinecke's avatar
Martin Reinecke committed
235
  else if (gname=="fejer1")
Martin Reinecke's avatar
Martin Reinecke committed
236
    {
Martin Reinecke's avatar
Martin Reinecke committed
237
238
239
    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
240
    if (verbose)
Martin Reinecke's avatar
Martin Reinecke committed
241
      cout << "Fejer1 grid, nlat=" << gpar1 << ", nlon=" << gpar2 << endl;
Martin Reinecke's avatar
Martin Reinecke committed
242
    }
Martin Reinecke's avatar
Martin Reinecke committed
243
  else if (gname=="fejer2")
Martin Reinecke's avatar
Martin Reinecke committed
244
    {
Martin Reinecke's avatar
Martin Reinecke committed
245
246
247
    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
248
    if (verbose)
Martin Reinecke's avatar
Martin Reinecke committed
249
      cout << "Fejer2 grid, nlat=" << gpar1 << ", nlon=" << gpar2 << endl;
Martin Reinecke's avatar
Martin Reinecke committed
250
    }
Martin Reinecke's avatar
Martin Reinecke committed
251
  else if (gname=="cc")
Martin Reinecke's avatar
Martin Reinecke committed
252
    {
Martin Reinecke's avatar
Martin Reinecke committed
253
254
255
    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
256
    if (verbose)
Martin Reinecke's avatar
Martin Reinecke committed
257
      cout << "Clenshaw-Curtis grid, nlat=" << gpar1 << ", nlon=" << gpar2 << endl;
Martin Reinecke's avatar
Martin Reinecke committed
258
    }
Martin Reinecke's avatar
Martin Reinecke committed
259
  else if (gname=="smallgauss")
Martin Reinecke's avatar
Martin Reinecke committed
260
    {
Martin Reinecke's avatar
Martin Reinecke committed
261
    int nlat=gpar1, nlon=gpar2;
Martin Reinecke's avatar
Martin Reinecke committed
262
    if (nlat<1) nlat=lmax+1;
Martin Reinecke's avatar
Martin Reinecke committed
263
264
265
    if (nlon<1) nlon=2*mmax+1;
    gpar1=nlat; gpar2=nlon;
    ginfo=sharp_make_gauss_geom_info (nlat, nlon, 0., 1, nlon);
Martin Reinecke's avatar
Martin Reinecke committed
266
267
    ptrdiff_t npix_o=get_npix(*ginfo);
    size_t ofs=0;
Martin Reinecke's avatar
Martin Reinecke committed
268
    for (int i=0; i<ginfo->pair.size(); ++i)
Martin Reinecke's avatar
Martin Reinecke committed
269
      {
Martin Reinecke's avatar
Martin Reinecke committed
270
271
      sharp_ringpair &pair(ginfo->pair[i]);
      int pring=1+2*sharp_get_mlim(lmax,0,pair.r1.sth,pair.r1.cth);
Martin Reinecke's avatar
Martin Reinecke committed
272
273
      if (pring>nlon) pring=nlon;
      pring=good_fft_size(pring);
Martin Reinecke's avatar
Martin Reinecke committed
274
275
276
      pair.r1.nph=pring;
      pair.r1.weight*=nlon*1./pring;
      pair.r1.ofs=ofs;
Martin Reinecke's avatar
Martin Reinecke committed
277
      ofs+=pring;
Martin Reinecke's avatar
Martin Reinecke committed
278
      if (pair.r2.nph>0)
Martin Reinecke's avatar
Martin Reinecke committed
279
        {
Martin Reinecke's avatar
Martin Reinecke committed
280
281
282
        pair.r2.nph=pring;
        pair.r2.weight*=nlon*1./pring;
        pair.r2.ofs=ofs;
Martin Reinecke's avatar
Martin Reinecke committed
283
284
285
286
287
288
        ofs+=pring;
        }
      }
    if (verbose)
      {
      ptrdiff_t npix=get_npix(*ginfo);
Martin Reinecke's avatar
Martin Reinecke committed
289
290
      cout << "Small Gauss grid, nlat=" << nlat << ", npix=" << npix
           << ", savings=" << ((npix_o-npix)*100./npix_o) << "\%" << endl;
Martin Reinecke's avatar
Martin Reinecke committed
291
292
293
294
295
296
297
298
299
300
301
302
303
      }
    }
  else
    MR_fail("unknown grid geometry");
  }

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

  /* flip theta to emulate the "old" Gaussian grid geometry */
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
307
  for (int i=0; i<tinfo->pair.size(); ++i)
Martin Reinecke's avatar
Martin Reinecke committed
308
309
310
311
312
313
314
315
    {
    const double pi=3.141592653589793238462643383279502884197;
    tinfo->pair[i].r1.cth=-tinfo->pair[i].r1.cth;
    tinfo->pair[i].r2.cth=-tinfo->pair[i].r2.cth;
    tinfo->pair[i].r1.theta=pi-tinfo->pair[i].r1.theta;
    tinfo->pair[i].r2.theta=pi-tinfo->pair[i].r2.theta;
    }

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

Martin Reinecke's avatar
Martin Reinecke committed
319
320
  vector<double> bmap(2*npix);
  vector<double *>map({&bmap[0], &bmap[npix]});
Martin Reinecke's avatar
Martin Reinecke committed
321

Martin Reinecke's avatar
Martin Reinecke committed
322
  vector<dcmplx> balm(2*nalms,dcmplx(1.,1.));
Martin Reinecke's avatar
Martin Reinecke committed
323
  vector<dcmplx *>alm({&balm[0], &balm[nalms]});
Martin Reinecke's avatar
Martin Reinecke committed
324

Martin Reinecke's avatar
Martin Reinecke committed
325
  sharp_execute(SHARP_ALM2MAP,0,alm.data(),map.data(),*tinfo,*alms,SHARP_DP,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
326
    nullptr,nullptr);
Martin Reinecke's avatar
Martin Reinecke committed
327
  MR_assert(approx(map[0][0     ], 3.588246976618616912e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
328
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
329
  MR_assert(approx(map[0][npix/2], 4.042209792157496651e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
330
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
331
  MR_assert(approx(map[0][npix-1],-1.234675107554816442e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
332
333
    "error");

Martin Reinecke's avatar
Martin Reinecke committed
334
  sharp_execute(SHARP_ALM2MAP,1,alm.data(),map.data(),*tinfo,*alms,SHARP_DP,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
335
    nullptr,nullptr);
Martin Reinecke's avatar
Martin Reinecke committed
336
  MR_assert(approx(map[0][0     ], 2.750897760535633285e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
337
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
338
  MR_assert(approx(map[0][npix/2], 3.137704477368562905e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
339
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
340
  MR_assert(approx(map[0][npix-1],-8.405730859837063917e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
341
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
342
  MR_assert(approx(map[1][0     ],-2.398026536095463346e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
343
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
344
  MR_assert(approx(map[1][npix/2],-4.961140548331700728e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
345
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
346
  MR_assert(approx(map[1][npix-1],-1.412765834230440021e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
347
348
    "error");

Martin Reinecke's avatar
Martin Reinecke committed
349
  sharp_execute(SHARP_ALM2MAP,2,alm.data(),map.data(),*tinfo,*alms,SHARP_DP,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
350
    nullptr,nullptr);
Martin Reinecke's avatar
Martin Reinecke committed
351
  MR_assert(approx(map[0][0     ],-1.398186224727334448e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
352
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
353
  MR_assert(approx(map[0][npix/2],-2.456676000884031197e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
354
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
355
  MR_assert(approx(map[0][npix-1],-1.516249174408820863e+02,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
356
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
357
  MR_assert(approx(map[1][0     ],-3.173406200299964119e+00,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
358
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
359
  MR_assert(approx(map[1][npix/2],-5.831327404513146462e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
360
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
361
  MR_assert(approx(map[1][npix-1],-1.863257892248353897e+01,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
362
363
    "error");

Martin Reinecke's avatar
Martin Reinecke committed
364
  sharp_execute(SHARP_ALM2MAP_DERIV1,1,alm.data(),map.data(),*tinfo,*alms,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
365
    SHARP_DP,nullptr,nullptr);
Martin Reinecke's avatar
Martin Reinecke committed
366
  MR_assert(approx(map[0][0     ],-6.859393905369091105e-01,1e-11),
Martin Reinecke's avatar
Martin Reinecke committed
367
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
368
  MR_assert(approx(map[0][npix/2],-2.103947835973212364e+02,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
369
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
370
  MR_assert(approx(map[0][npix-1],-1.092463246472086439e+03,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
371
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
372
  MR_assert(approx(map[1][0     ],-1.411433220713928165e+02,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
373
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
374
  MR_assert(approx(map[1][npix/2],-1.146122859381925082e+03,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
375
    "error");
Martin Reinecke's avatar
Martin Reinecke committed
376
  MR_assert(approx(map[1][npix-1], 7.821618677689795049e+02,1e-12),
Martin Reinecke's avatar
Martin Reinecke committed
377
378
379
    "error");
  }

Martin Reinecke's avatar
Martin Reinecke committed
380
static void do_sht (sharp_geom_info &ginfo, sharp_alm_info &ainfo,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
381
  int spin, vector<double> &err_abs, vector<double> &err_rel,
Martin Reinecke's avatar
Martin Reinecke committed
382
383
384
385
386
387
388
  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);
  int ncomp = (spin==0) ? 1 : 2;

  size_t npix = get_npix(ginfo);
Martin Reinecke's avatar
Martin Reinecke committed
389
390
  vector<double> bmap(ncomp*ntrans*npix, 0.);
  vector<double *>map(ncomp*ntrans);
Martin Reinecke's avatar
Martin Reinecke committed
391
  for (int i=0; i<ncomp*ntrans; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
392
    map[i]=&bmap[i*npix];
Martin Reinecke's avatar
Martin Reinecke committed
393

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

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

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
405
406
  if (t_a2m!=nullptr) *t_a2m=0;
  if (op_a2m!=nullptr) *op_a2m=0;
Martin Reinecke's avatar
Martin Reinecke committed
407
408
409
410
  for (size_t itrans=0; itrans<ntrans; ++itrans)
    {
    sharp_execute(SHARP_ALM2MAP,spin,&alm[itrans*ncomp],&map[itrans*ncomp],ginfo,ainfo,
      SHARP_DP,&tta2m,&toa2m);
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
411
412
    if (t_a2m!=nullptr) *t_a2m+=maxTime(tta2m);
    if (op_a2m!=nullptr) *op_a2m+=totalops(toa2m);
Martin Reinecke's avatar
Martin Reinecke committed
413
    }
Martin Reinecke's avatar
Martin Reinecke committed
414
  auto sqsum=get_sqsum_and_invert(alm.data(),nalms,ntrans*ncomp);
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
415
416
  if (t_m2a!=nullptr) *t_m2a=0;
  if (op_m2a!=nullptr) *op_m2a=0;
Martin Reinecke's avatar
Martin Reinecke committed
417
418
419
420
  for (size_t itrans=0; itrans<ntrans; ++itrans)
    {
    sharp_execute(SHARP_MAP2ALM,spin,&alm[itrans*ncomp],&map[itrans*ncomp],ginfo,ainfo,
      SHARP_DP|SHARP_ADD,&ttm2a,&tom2a);
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
421
422
    if (t_m2a!=nullptr) *t_m2a+=maxTime(ttm2a);
    if (op_m2a!=nullptr) *op_m2a+=totalops(tom2a);
Martin Reinecke's avatar
Martin Reinecke committed
423
    }
Martin Reinecke's avatar
Martin Reinecke committed
424
  get_errors(alm.data(), nalms, ntrans*ncomp, sqsum, err_abs, err_rel);
Martin Reinecke's avatar
Martin Reinecke committed
425
426
  }

Martin Reinecke's avatar
Martin Reinecke committed
427
static void check_accuracy (sharp_geom_info &ginfo, sharp_alm_info &ainfo,
Martin Reinecke's avatar
Martin Reinecke committed
428
429
430
  int spin)
  {
  int ncomp = (spin==0) ? 1 : 2;
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
431
  vector<double> err_abs, err_rel;
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
432
433
  do_sht (ginfo, ainfo, spin, err_abs, err_rel, nullptr, nullptr,
    nullptr, nullptr, 1);
Martin Reinecke's avatar
Martin Reinecke committed
434
435
436
437
438
439
  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
440
441
442
443
  unique_ptr<sharp_geom_info> ginfo;
  unique_ptr<sharp_alm_info> ainfo;
  get_infos ("gauss", lmax, mmax, nlat, nlon, ginfo, ainfo, 0);
  check_accuracy(*ginfo,*ainfo,spin);
Martin Reinecke's avatar
Martin Reinecke committed
444
445
446
447
448
449
  }

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

Martin Reinecke's avatar
Martin Reinecke committed
450
  if (mytask==0) cout << "Checking signs and scales.\n";
Martin Reinecke's avatar
Martin Reinecke committed
451
  check_sign_scale();
Martin Reinecke's avatar
Martin Reinecke committed
452
  if (mytask==0) cout << "Passed.\n\n";
Martin Reinecke's avatar
Martin Reinecke committed
453

Martin Reinecke's avatar
Martin Reinecke committed
454
  if (mytask==0) cout << "Testing map analysis accuracy.\n";
Martin Reinecke's avatar
Martin Reinecke committed
455
456
457
458
459
460
461
462
463
464

  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
465
  if (mytask==0) cout << "Passed.\n\n";
Martin Reinecke's avatar
Martin Reinecke committed
466
467
468
469
470
471
472
473
474
475
476
477
478
479
  }

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]");
  int lmax=atoi(argv[3]);
  int mmax=atoi(argv[4]);
  int gpar1=atoi(argv[5]);
  int gpar2=atoi(argv[6]);
  int spin=atoi(argv[7]);
  int ntrans=1;
  if (argc>=9) ntrans=atoi(argv[8]);

Martin Reinecke's avatar
Martin Reinecke committed
480
481
482
  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
483

Martin Reinecke's avatar
Martin Reinecke committed
484
485
486
  unique_ptr<sharp_geom_info> ginfo;
  unique_ptr<sharp_alm_info> ainfo;
  get_infos (argv[2], lmax, mmax, gpar1, gpar2, ginfo, ainfo, 1);
Martin Reinecke's avatar
Martin Reinecke committed
487
488
489
490

  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
491
  vector<double> err_abs, err_rel;
Martin Reinecke's avatar
Martin Reinecke committed
492
493
494
495
496
497
498

  double t_acc=0;
  int nrpt=0;
  while(1)
    {
    ++nrpt;
    double ta2m2, tm2a2;
Martin Reinecke's avatar
Martin Reinecke committed
499
    do_sht (*ginfo, *ainfo, spin, err_abs, err_rel, &ta2m2, &tm2a2,
Martin Reinecke's avatar
Martin Reinecke committed
500
501
502
503
504
505
      &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
506
      if (mytask==0) cout << "Best of " << nrpt << " runs\n";
Martin Reinecke's avatar
Martin Reinecke committed
507
508
509
510
      break;
      }
    }

Martin Reinecke's avatar
Martin Reinecke committed
511
512
513
514
  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
515
516
517

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

Martin Reinecke's avatar
Martin Reinecke committed
520
  double iosize = ntrans*ncomp*(16.*get_nalms(*ainfo) + 8.*get_npix(*ginfo));
Martin Reinecke's avatar
Martin Reinecke committed
521
522
523
524
  iosize = allreduceSumDouble(iosize);

  double tmem=totalMem();
  if (mytask==0)
Martin Reinecke's avatar
Martin Reinecke committed
525
    cout << "\nMemory high water mark: " << tmem/(1<<20) << " MB\n";
Martin Reinecke's avatar
Martin Reinecke committed
526
  if (mytask==0)
Martin Reinecke's avatar
Martin Reinecke committed
527
528
    cout << "Memory overhead: " << (tmem-iosize)/(1<<20) << " MB ("
         << 100.*(1.-iosize/tmem) << "\% of working set)\n";
Martin Reinecke's avatar
Martin Reinecke committed
529
530
531
532
533
534
535
536
537
538
539

  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)
  {
540
541
542
543
544
545
546
  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
547
548
549
  mytask=0; ntasks=1;

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

Martin Reinecke's avatar
Martin Reinecke committed
552
  if (mode=="acctest")
Martin Reinecke's avatar
Martin Reinecke committed
553
    sharp_acctest();
Martin Reinecke's avatar
Martin Reinecke committed
554
  else if (mode=="test")
Martin Reinecke's avatar
Martin Reinecke committed
555
556
557
558
559
560
    sharp_test(argc,argv);
  else
    MR_fail("unknown command");

  return 0;
  }