NFFT  3.3.0
quadratureS2.c
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id$ */
20 
26 #include "config.h"
27 
28 /* Include standard C headers. */
29 #include <math.h>
30 #include <stdlib.h>
31 #include <stdio.h>
32 #include <string.h>
33 #include <time.h>
34 #ifdef HAVE_COMPLEX_H
35 #include <complex.h>
36 #endif
37 
38 /* Include NFFT 3 utilities headers. */
39 
40 /* Include NFFT 3 library header. */
41 #include "nfft3.h"
42 
43 #include "infft.h"
44 
46 enum boolean {NO = 0, YES = 1};
47 
49 enum testtype {ERROR = 0, TIMING = 1};
50 
52 enum gridtype {GRID_GAUSS_LEGENDRE = 0, GRID_CLENSHAW_CURTIS = 1,
53  GRID_HEALPIX = 2, GRID_EQUIDISTRIBUTION = 3, GRID_EQUIDISTRIBUTION_UNIFORM = 4};
54 
56 enum functiontype {FUNCTION_RANDOM_BANDLIMITED = 0, FUNCTION_F1 = 1,
57  FUNCTION_F2 = 2, FUNCTION_F3 = 3, FUNCTION_F4 = 4, FUNCTION_F5 = 5,
58  FUNCTION_F6 = 6};
59 
61 enum modes {USE_GRID = 0, RANDOM = 1};
62 
71 int main (int argc, char **argv)
72 {
73  int tc;
74  int tc_max;
76  int *NQ;
78  int NQ_max;
80  int *SQ;
82  int SQ_max;
83  int *RQ;
85  int iNQ;
86  int iNQ_max;
87  int testfunction;
88  int N;
90  int use_nfsft;
91  int use_nfft;
92  int use_fpt;
93  int cutoff;
94  double threshold;
96  int gridtype;
97  int repetitions;
98  int mode;
100  double *w;
101  double *x_grid;
102  double *x_compare;
103  double _Complex *f_grid;
104  double _Complex *f_compare;
105  double _Complex *f;
106  double _Complex *f_hat_gen;
108  double _Complex *f_hat;
110  nfsft_plan plan_adjoint;
111  nfsft_plan plan;
112  nfsft_plan plan_gen;
114  double t_avg;
115  double err_infty_avg;
116  double err_2_avg;
118  int i;
119  int k;
120  int n;
121  int d;
123  int m_theta;
125  int m_phi;
127  int m_total;
128  double *theta;
130  double *phi;
132  fftw_plan fplan;
134  //int nside; /**< The size parameter for the HEALPix grid */
135  int d2;
136  int M;
137  double theta_s;
138  double x1,x2,x3,temp;
139  int m_compare;
140  nfsft_plan *plan_adjoint_ptr;
141  nfsft_plan *plan_ptr;
142  double *w_temp;
143  int testmode;
144  ticks t0, t1;
145 
146  /* Read the number of testcases. */
147  fscanf(stdin,"testcases=%d\n",&tc_max);
148  fprintf(stdout,"%d\n",tc_max);
149 
150  /* Process each testcase. */
151  for (tc = 0; tc < tc_max; tc++)
152  {
153  /* Check if the fast transform shall be used. */
154  fscanf(stdin,"nfsft=%d\n",&use_nfsft);
155  fprintf(stdout,"%d\n",use_nfsft);
156  if (use_nfsft != NO)
157  {
158  /* Check if the NFFT shall be used. */
159  fscanf(stdin,"nfft=%d\n",&use_nfft);
160  fprintf(stdout,"%d\n",use_nfsft);
161  if (use_nfft != NO)
162  {
163  /* Read the cut-off parameter. */
164  fscanf(stdin,"cutoff=%d\n",&cutoff);
165  fprintf(stdout,"%d\n",cutoff);
166  }
167  else
168  {
169  /* TODO remove this */
170  /* Initialize unused variable with dummy value. */
171  cutoff = 1;
172  }
173  /* Check if the fast polynomial transform shall be used. */
174  fscanf(stdin,"fpt=%d\n",&use_fpt);
175  fprintf(stdout,"%d\n",use_fpt);
176  if (use_fpt != NO)
177  {
178  /* Read the NFSFT threshold parameter. */
179  fscanf(stdin,"threshold=%lf\n",&threshold);
180  fprintf(stdout,"%lf\n",threshold);
181  }
182  else
183  {
184  /* TODO remove this */
185  /* Initialize unused variable with dummy value. */
186  threshold = 1000.0;
187  }
188  }
189  else
190  {
191  /* TODO remove this */
192  /* Set dummy values. */
193  use_nfft = NO;
194  use_fpt = NO;
195  cutoff = 3;
196  threshold = 1000.0;
197  }
198 
199  /* Read the testmode type. */
200  fscanf(stdin,"testmode=%d\n",&testmode);
201  fprintf(stdout,"%d\n",testmode);
202 
203  if (testmode == ERROR)
204  {
205  /* Read the quadrature grid type. */
206  fscanf(stdin,"gridtype=%d\n",&gridtype);
207  fprintf(stdout,"%d\n",gridtype);
208 
209  /* Read the test function. */
210  fscanf(stdin,"testfunction=%d\n",&testfunction);
211  fprintf(stdout,"%d\n",testfunction);
212 
213  /* Check if random bandlimited function has been chosen. */
214  if (testfunction == FUNCTION_RANDOM_BANDLIMITED)
215  {
216  /* Read the bandwidht. */
217  fscanf(stdin,"bandlimit=%d\n",&N);
218  fprintf(stdout,"%d\n",N);
219  }
220  else
221  {
222  N = 1;
223  }
224 
225  /* Read the number of repetitions. */
226  fscanf(stdin,"repetitions=%d\n",&repetitions);
227  fprintf(stdout,"%d\n",repetitions);
228 
229  fscanf(stdin,"mode=%d\n",&mode);
230  fprintf(stdout,"%d\n",mode);
231 
232  if (mode == RANDOM)
233  {
234  /* Read the bandwidht. */
235  fscanf(stdin,"points=%d\n",&m_compare);
236  fprintf(stdout,"%d\n",m_compare);
237  x_compare = (double*) nfft_malloc(2*m_compare*sizeof(double));
238  d = 0;
239  while (d < m_compare)
240  {
241  x1 = 2.0*(((double)rand())/RAND_MAX) - 1.0;
242  x2 = 2.0*(((double)rand())/RAND_MAX) - 1.0;
243  x3 = 2.0*(((double)rand())/RAND_MAX) - 1.0;
244  temp = sqrt(x1*x1+x2*x2+x3*x3);
245  if (temp <= 1)
246  {
247  x_compare[2*d+1] = acos(x3);
248  if (x_compare[2*d+1] == 0 || x_compare[2*d+1] == KPI)
249  {
250  x_compare[2*d] = 0.0;
251  }
252  else
253  {
254  x_compare[2*d] = atan2(x2/sin(x_compare[2*d+1]),x1/sin(x_compare[2*d+1]));
255  }
256  x_compare[2*d] *= 1.0/(2.0*KPI);
257  x_compare[2*d+1] *= 1.0/(2.0*KPI);
258  d++;
259  }
260  }
261  f_compare = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex));
262  f = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex));
263  }
264  }
265 
266  /* Initialize maximum cut-off degree and grid size parameter. */
267  NQ_max = 0;
268  SQ_max = 0;
269 
270  /* Read the number of cut-off degrees. */
271  fscanf(stdin,"bandwidths=%d\n",&iNQ_max);
272  fprintf(stdout,"%d\n",iNQ_max);
273 
274  /* Allocate memory for the cut-off degrees and grid size parameters. */
275  NQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
276  SQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
277  if (testmode == TIMING)
278  {
279  RQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
280  }
281 
282  /* Read the cut-off degrees and grid size parameters. */
283  for (iNQ = 0; iNQ < iNQ_max; iNQ++)
284  {
285  if (testmode == TIMING)
286  {
287  /* Read cut-off degree and grid size parameter. */
288  fscanf(stdin,"%d %d %d\n",&NQ[iNQ],&SQ[iNQ],&RQ[iNQ]);
289  fprintf(stdout,"%d %d %d\n",NQ[iNQ],SQ[iNQ],RQ[iNQ]);
290  NQ_max = MAX(NQ_max,NQ[iNQ]);
291  SQ_max = MAX(SQ_max,SQ[iNQ]);
292  }
293  else
294  {
295  /* Read cut-off degree and grid size parameter. */
296  fscanf(stdin,"%d %d\n",&NQ[iNQ],&SQ[iNQ]);
297  fprintf(stdout,"%d %d\n",NQ[iNQ],SQ[iNQ]);
298  NQ_max = MAX(NQ_max,NQ[iNQ]);
299  SQ_max = MAX(SQ_max,SQ[iNQ]);
300  }
301  }
302 
303  /* Do precomputation. */
304  //fprintf(stderr,"NFSFT Precomputation\n");
305  //fflush(stderr);
306  nfsft_precompute(NQ_max, threshold,
307  ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U)), 0U);
308 
309  if (testmode == TIMING)
310  {
311  /* Allocate data structures. */
312  f_hat = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(NQ_max)*sizeof(double _Complex));
313  f = (double _Complex*) nfft_malloc(SQ_max*sizeof(double _Complex));
314  x_grid = (double*) nfft_malloc(2*SQ_max*sizeof(double));
315  for (d = 0; d < SQ_max; d++)
316  {
317  f[d] = (((double)rand())/RAND_MAX)-0.5 + _Complex_I*((((double)rand())/RAND_MAX)-0.5);
318  x_grid[2*d] = (((double)rand())/RAND_MAX) - 0.5;
319  x_grid[2*d+1] = (((double)rand())/RAND_MAX) * 0.5;
320  }
321  }
322 
323  //fprintf(stderr,"Entering loop\n");
324  //fflush(stderr);
325  /* Process all cut-off bandwidths. */
326  for (iNQ = 0; iNQ < iNQ_max; iNQ++)
327  {
328  if (testmode == TIMING)
329  {
330  nfsft_init_guru(&plan,NQ[iNQ],SQ[iNQ], NFSFT_NORMALIZED |
331  ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
332  ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
333  PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFTW_MEASURE | FFT_OUT_OF_PLACE,
334  cutoff);
335 
336  plan.f_hat = f_hat;
337  plan.x = x_grid;
338  plan.f = f;
339 
340  nfsft_precompute_x(&plan);
341 
342  t_avg = 0.0;
343 
344  for (i = 0; i < RQ[iNQ]; i++)
345  {
346  t0 = getticks();
347 
348  if (use_nfsft != NO)
349  {
350  /* Execute the adjoint NFSFT transformation. */
351  nfsft_adjoint(&plan);
352  }
353  else
354  {
355  /* Execute the adjoint direct NDSFT transformation. */
356  nfsft_adjoint_direct(&plan);
357  }
358 
359  t1 = getticks();
360  t_avg += nfft_elapsed_seconds(t1,t0);
361  }
362 
363  t_avg = t_avg/((double)RQ[iNQ]);
364 
365  nfsft_finalize(&plan);
366 
367  fprintf(stdout,"%+le\n", t_avg);
368  fprintf(stderr,"%d: %4d %4d %+le\n", tc, NQ[iNQ], SQ[iNQ], t_avg);
369  }
370  else
371  {
372  /* Determine the maximum number of nodes. */
373  switch (gridtype)
374  {
375  case GRID_GAUSS_LEGENDRE:
376  /* Calculate grid dimensions. */
377  m_theta = SQ[iNQ] + 1;
378  m_phi = 2*SQ[iNQ] + 2;
379  m_total = m_theta*m_phi;
380  break;
381  case GRID_CLENSHAW_CURTIS:
382  /* Calculate grid dimensions. */
383  m_theta = 2*SQ[iNQ] + 1;
384  m_phi = 2*SQ[iNQ] + 2;
385  m_total = m_theta*m_phi;
386  break;
387  case GRID_HEALPIX:
388  m_theta = 1;
389  m_phi = 12*SQ[iNQ]*SQ[iNQ];
390  m_total = m_theta * m_phi;
391  //fprintf("HEALPix: SQ = %d, m_theta = %d, m_phi= %d, m");
392  break;
393  case GRID_EQUIDISTRIBUTION:
394  case GRID_EQUIDISTRIBUTION_UNIFORM:
395  m_theta = 2;
396  //fprintf(stderr,"ed: m_theta = %d\n",m_theta);
397  for (k = 1; k < SQ[iNQ]; k++)
398  {
399  m_theta += (int)floor((2*KPI)/acos((cos(KPI/(double)SQ[iNQ])-
400  cos(k*KPI/(double)SQ[iNQ])*cos(k*KPI/(double)SQ[iNQ]))/
401  (sin(k*KPI/(double)SQ[iNQ])*sin(k*KPI/(double)SQ[iNQ]))));
402  //fprintf(stderr,"ed: m_theta = %d\n",m_theta);
403  }
404  //fprintf(stderr,"ed: m_theta final = %d\n",m_theta);
405  m_phi = 1;
406  m_total = m_theta * m_phi;
407  break;
408  }
409 
410  /* Allocate memory for data structures. */
411  w = (double*) nfft_malloc(m_theta*sizeof(double));
412  x_grid = (double*) nfft_malloc(2*m_total*sizeof(double));
413 
414  //fprintf(stderr,"NQ = %d\n",NQ[iNQ]);
415  //fflush(stderr);
416  switch (gridtype)
417  {
418  case GRID_GAUSS_LEGENDRE:
419  //fprintf(stderr,"Generating grid for NQ = %d, SQ = %d\n",NQ[iNQ],SQ[iNQ]);
420  //fflush(stderr);
421 
422  /* Read quadrature weights. */
423  for (k = 0; k < m_theta; k++)
424  {
425  fscanf(stdin,"%le\n",&w[k]);
426  w[k] *= (2.0*KPI)/((double)m_phi);
427  }
428 
429  //fprintf(stderr,"Allocating theta and phi\n");
430  //fflush(stderr);
431  /* Allocate memory to store the grid's angles. */
432  theta = (double*) nfft_malloc(m_theta*sizeof(double));
433  phi = (double*) nfft_malloc(m_phi*sizeof(double));
434 
435  //if (theta == NULL || phi == NULL)
436  //{
437  //fprintf(stderr,"Couldn't allocate theta and phi\n");
438  //fflush(stderr);
439  //}
440 
441 
442  /* Read angles theta. */
443  for (k = 0; k < m_theta; k++)
444  {
445  fscanf(stdin,"%le\n",&theta[k]);
446  }
447 
448  /* Generate the grid angles phi. */
449  for (n = 0; n < m_phi; n++)
450  {
451  phi[n] = n/((double)m_phi);
452  phi[n] -= ((phi[n]>=0.5)?(1.0):(0.0));
453  }
454 
455  //fprintf(stderr,"Generating grid nodes\n");
456  //fflush(stderr);
457 
458  /* Generate the grid's nodes. */
459  d = 0;
460  for (k = 0; k < m_theta; k++)
461  {
462  for (n = 0; n < m_phi; n++)
463  {
464  x_grid[2*d] = phi[n];
465  x_grid[2*d+1] = theta[k];
466  d++;
467  }
468  }
469 
470  //fprintf(stderr,"Freeing theta and phi\n");
471  //fflush(stderr);
472  /* Free the arrays for the grid's angles. */
473  nfft_free(theta);
474  nfft_free(phi);
475 
476  break;
477 
478  case GRID_CLENSHAW_CURTIS:
479 
480  /* Allocate memory to store the grid's angles. */
481  theta = (double*) nfft_malloc(m_theta*sizeof(double));
482  phi = (double*) nfft_malloc(m_phi*sizeof(double));
483 
484  /* Generate the grid angles theta. */
485  for (k = 0; k < m_theta; k++)
486  {
487  theta[k] = k/((double)2*(m_theta-1));
488  }
489 
490  /* Generate the grid angles phi. */
491  for (n = 0; n < m_phi; n++)
492  {
493  phi[n] = n/((double)m_phi);
494  phi[n] -= ((phi[n]>=0.5)?(1.0):(0.0));
495  }
496 
497  /* Generate quadrature weights. */
498  fplan = fftw_plan_r2r_1d(SQ[iNQ]+1, w, w, FFTW_REDFT00, 0U);
499  for (k = 0; k < SQ[iNQ]+1; k++)
500  {
501  w[k] = -2.0/(4*k*k-1);
502  }
503  fftw_execute(fplan);
504  w[0] *= 0.5;
505 
506  for (k = 0; k < SQ[iNQ]+1; k++)
507  {
508  w[k] *= (2.0*KPI)/((double)(m_theta-1)*m_phi);
509  w[m_theta-1-k] = w[k];
510  }
511  fftw_destroy_plan(fplan);
512 
513  /* Generate the grid's nodes. */
514  d = 0;
515  for (k = 0; k < m_theta; k++)
516  {
517  for (n = 0; n < m_phi; n++)
518  {
519  x_grid[2*d] = phi[n];
520  x_grid[2*d+1] = theta[k];
521  d++;
522  }
523  }
524 
525  /* Free the arrays for the grid's angles. */
526  nfft_free(theta);
527  nfft_free(phi);
528 
529  break;
530 
531  case GRID_HEALPIX:
532 
533  d = 0;
534  for (k = 1; k <= SQ[iNQ]-1; k++)
535  {
536  for (n = 0; n <= 4*k-1; n++)
537  {
538  x_grid[2*d+1] = 1 - (k*k)/((double)(3.0*SQ[iNQ]*SQ[iNQ]));
539  x_grid[2*d] = ((n+0.5)/(4*k));
540  x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0);
541  d++;
542  }
543  }
544 
545  d2 = d-1;
546 
547  for (k = SQ[iNQ]; k <= 3*SQ[iNQ]; k++)
548  {
549  for (n = 0; n <= 4*SQ[iNQ]-1; n++)
550  {
551  x_grid[2*d+1] = 2.0/(3*SQ[iNQ])*(2*SQ[iNQ]-k);
552  x_grid[2*d] = (n+((k%2==0)?(0.5):(0.0)))/(4*SQ[iNQ]);
553  x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0);
554  d++;
555  }
556  }
557 
558  for (k = 1; k <= SQ[iNQ]-1; k++)
559  {
560  for (n = 0; n <= 4*k-1; n++)
561  {
562  x_grid[2*d+1] = -x_grid[2*d2+1];
563  x_grid[2*d] = x_grid[2*d2];
564  d++;
565  d2--;
566  }
567  }
568 
569  for (d = 0; d < m_total; d++)
570  {
571  x_grid[2*d+1] = acos(x_grid[2*d+1])/(2.0*KPI);
572  }
573 
574  w[0] = (4.0*KPI)/(m_total);
575  break;
576 
577  case GRID_EQUIDISTRIBUTION:
578  case GRID_EQUIDISTRIBUTION_UNIFORM:
579  /* TODO Compute the weights. */
580 
581  if (gridtype == GRID_EQUIDISTRIBUTION)
582  {
583  w_temp = (double*) nfft_malloc((SQ[iNQ]+1)*sizeof(double));
584  fplan = fftw_plan_r2r_1d(SQ[iNQ]/2+1, w_temp, w_temp, FFTW_REDFT00, 0U);
585  for (k = 0; k < SQ[iNQ]/2+1; k++)
586  {
587  w_temp[k] = -2.0/(4*k*k-1);
588  }
589  fftw_execute(fplan);
590  w_temp[0] *= 0.5;
591 
592  for (k = 0; k < SQ[iNQ]/2+1; k++)
593  {
594  w_temp[k] *= (2.0*KPI)/((double)(SQ[iNQ]));
595  w_temp[SQ[iNQ]-k] = w_temp[k];
596  }
597  fftw_destroy_plan(fplan);
598  }
599 
600  d = 0;
601  x_grid[2*d] = -0.5;
602  x_grid[2*d+1] = 0.0;
603  if (gridtype == GRID_EQUIDISTRIBUTION)
604  {
605  w[d] = w_temp[0];
606  }
607  else
608  {
609  w[d] = (4.0*KPI)/(m_total);
610  }
611  d = 1;
612  x_grid[2*d] = -0.5;
613  x_grid[2*d+1] = 0.5;
614  if (gridtype == GRID_EQUIDISTRIBUTION)
615  {
616  w[d] = w_temp[SQ[iNQ]];
617  }
618  else
619  {
620  w[d] = (4.0*KPI)/(m_total);
621  }
622  d = 2;
623 
624  for (k = 1; k < SQ[iNQ]; k++)
625  {
626  theta_s = (double)k*KPI/(double)SQ[iNQ];
627  M = (int)floor((2.0*KPI)/acos((cos(KPI/(double)SQ[iNQ])-
628  cos(theta_s)*cos(theta_s))/(sin(theta_s)*sin(theta_s))));
629 
630  for (n = 0; n < M; n++)
631  {
632  x_grid[2*d] = (n + 0.5)/M;
633  x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0);
634  x_grid[2*d+1] = theta_s/(2.0*KPI);
635  if (gridtype == GRID_EQUIDISTRIBUTION)
636  {
637  w[d] = w_temp[k]/((double)(M));
638  }
639  else
640  {
641  w[d] = (4.0*KPI)/(m_total);
642  }
643  d++;
644  }
645  }
646 
647  if (gridtype == GRID_EQUIDISTRIBUTION)
648  {
649  nfft_free(w_temp);
650  }
651  break;
652 
653  default:
654  break;
655  }
656 
657  /* Allocate memory for grid values. */
658  f_grid = (double _Complex*) nfft_malloc(m_total*sizeof(double _Complex));
659 
660  if (mode == RANDOM)
661  {
662  }
663  else
664  {
665  m_compare = m_total;
666  f_compare = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex));
667  x_compare = x_grid;
668  f = f_grid;
669  }
670 
671  //fprintf(stderr,"Generating test function\n");
672  //fflush(stderr);
673  switch (testfunction)
674  {
675  case FUNCTION_RANDOM_BANDLIMITED:
676  f_hat_gen = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(N)*sizeof(double _Complex));
677  //fprintf(stderr,"Generating random test function\n");
678  //fflush(stderr);
679  /* Generate random function samples by sampling a bandlimited
680  * function. */
681  nfsft_init_guru(&plan_gen,N,m_total, NFSFT_NORMALIZED |
682  ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
683  ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
684  ((N>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
685  FFT_OUT_OF_PLACE, cutoff);
686 
687  plan_gen.f_hat = f_hat_gen;
688  plan_gen.x = x_grid;
689  plan_gen.f = f_grid;
690 
691  nfsft_precompute_x(&plan_gen);
692 
693  for (k = 0; k < plan_gen.N_total; k++)
694  {
695  f_hat_gen[k] = 0.0;
696  }
697 
698  for (k = 0; k <= N; k++)
699  {
700  for (n = -k; n <= k; n++)
701  {
702  f_hat_gen[NFSFT_INDEX(k,n,&plan_gen)] =
703  (((double)rand())/RAND_MAX)-0.5 + _Complex_I*((((double)rand())/RAND_MAX)-0.5);
704  }
705  }
706 
707  if (use_nfsft != NO)
708  {
709  /* Execute the NFSFT transformation. */
710  nfsft_trafo(&plan_gen);
711  }
712  else
713  {
714  /* Execute the direct NDSFT transformation. */
715  nfsft_trafo_direct(&plan_gen);
716  }
717 
718  nfsft_finalize(&plan_gen);
719 
720  if (mode == RANDOM)
721  {
722  nfsft_init_guru(&plan_gen,N,m_compare, NFSFT_NORMALIZED |
723  ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
724  ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
725  ((N>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
726  FFT_OUT_OF_PLACE, cutoff);
727 
728  plan_gen.f_hat = f_hat_gen;
729  plan_gen.x = x_compare;
730  plan_gen.f = f_compare;
731 
732  nfsft_precompute_x(&plan_gen);
733 
734  if (use_nfsft != NO)
735  {
736  /* Execute the NFSFT transformation. */
737  nfsft_trafo(&plan_gen);
738  }
739  else
740  {
741  /* Execute the direct NDSFT transformation. */
742  nfsft_trafo_direct(&plan_gen);
743  }
744 
745  nfsft_finalize(&plan_gen);
746  }
747  else
748  {
749  memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
750  }
751 
752  nfft_free(f_hat_gen);
753 
754  break;
755 
756  case FUNCTION_F1:
757  for (d = 0; d < m_total; d++)
758  {
759  x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI);
760  x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI);
761  x3 = cos(x_grid[2*d+1]*2.0*KPI);
762  f_grid[d] = x1*x2*x3;
763  }
764  if (mode == RANDOM)
765  {
766  for (d = 0; d < m_compare; d++)
767  {
768  x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI);
769  x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI);
770  x3 = cos(x_compare[2*d+1]*2.0*KPI);
771  f_compare[d] = x1*x2*x3;
772  }
773  }
774  else
775  {
776  memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
777  }
778  break;
779  case FUNCTION_F2:
780  for (d = 0; d < m_total; d++)
781  {
782  x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI);
783  x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI);
784  x3 = cos(x_grid[2*d+1]*2.0*KPI);
785  f_grid[d] = 0.1*exp(x1+x2+x3);
786  }
787  if (mode == RANDOM)
788  {
789  for (d = 0; d < m_compare; d++)
790  {
791  x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI);
792  x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI);
793  x3 = cos(x_compare[2*d+1]*2.0*KPI);
794  f_compare[d] = 0.1*exp(x1+x2+x3);
795  }
796  }
797  else
798  {
799  memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
800  }
801  break;
802  case FUNCTION_F3:
803  for (d = 0; d < m_total; d++)
804  {
805  x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI);
806  x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI);
807  x3 = cos(x_grid[2*d+1]*2.0*KPI);
808  temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
809  f_grid[d] = 0.1*temp;
810  }
811  if (mode == RANDOM)
812  {
813  for (d = 0; d < m_compare; d++)
814  {
815  x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI);
816  x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI);
817  x3 = cos(x_compare[2*d+1]*2.0*KPI);
818  temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
819  f_compare[d] = 0.1*temp;
820  }
821  }
822  else
823  {
824  memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
825  }
826  break;
827  case FUNCTION_F4:
828  for (d = 0; d < m_total; d++)
829  {
830  x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI);
831  x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI);
832  x3 = cos(x_grid[2*d+1]*2.0*KPI);
833  temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
834  f_grid[d] = 1.0/(temp);
835  }
836  if (mode == RANDOM)
837  {
838  for (d = 0; d < m_compare; d++)
839  {
840  x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI);
841  x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI);
842  x3 = cos(x_compare[2*d+1]*2.0*KPI);
843  temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
844  f_compare[d] = 1.0/(temp);
845  }
846  }
847  else
848  {
849  memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
850  }
851  break;
852  case FUNCTION_F5:
853  for (d = 0; d < m_total; d++)
854  {
855  x1 = sin(x_grid[2*d+1]*2.0*KPI)*cos(x_grid[2*d]*2.0*KPI);
856  x2 = sin(x_grid[2*d+1]*2.0*KPI)*sin(x_grid[2*d]*2.0*KPI);
857  x3 = cos(x_grid[2*d+1]*2.0*KPI);
858  temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
859  f_grid[d] = 0.1*sin(1+temp)*sin(1+temp);
860  }
861  if (mode == RANDOM)
862  {
863  for (d = 0; d < m_compare; d++)
864  {
865  x1 = sin(x_compare[2*d+1]*2.0*KPI)*cos(x_compare[2*d]*2.0*KPI);
866  x2 = sin(x_compare[2*d+1]*2.0*KPI)*sin(x_compare[2*d]*2.0*KPI);
867  x3 = cos(x_compare[2*d+1]*2.0*KPI);
868  temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
869  f_compare[d] = 0.1*sin(1+temp)*sin(1+temp);
870  }
871  }
872  else
873  {
874  memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
875  }
876  break;
877  case FUNCTION_F6:
878  for (d = 0; d < m_total; d++)
879  {
880  if (x_grid[2*d+1] <= 0.25)
881  {
882  f_grid[d] = 1.0;
883  }
884  else
885  {
886  f_grid[d] = 1.0/(sqrt(1+3*cos(2.0*KPI*x_grid[2*d+1])*cos(2.0*KPI*x_grid[2*d+1])));
887  }
888  }
889  if (mode == RANDOM)
890  {
891  for (d = 0; d < m_compare; d++)
892  {
893  if (x_compare[2*d+1] <= 0.25)
894  {
895  f_compare[d] = 1.0;
896  }
897  else
898  {
899  f_compare[d] = 1.0/(sqrt(1+3*cos(2.0*KPI*x_compare[2*d+1])*cos(2.0*KPI*x_compare[2*d+1])));
900  }
901  }
902  }
903  else
904  {
905  memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
906  }
907  break;
908  default:
909  //fprintf(stderr,"Generating one function\n");
910  //fflush(stderr);
911  for (d = 0; d < m_total; d++)
912  {
913  f_grid[d] = 1.0;
914  }
915  if (mode == RANDOM)
916  {
917  for (d = 0; d < m_compare; d++)
918  {
919  f_compare[d] = 1.0;
920  }
921  }
922  else
923  {
924  memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
925  }
926  break;
927  }
928 
929  //fprintf(stderr,"Initializing trafo\n");
930  //fflush(stderr);
931  /* Init transform plan. */
932  nfsft_init_guru(&plan_adjoint,NQ[iNQ],m_total, NFSFT_NORMALIZED |
933  ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
934  ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
935  ((NQ[iNQ]>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
936  FFT_OUT_OF_PLACE, cutoff);
937 
938  plan_adjoint_ptr = &plan_adjoint;
939 
940  if (mode == RANDOM)
941  {
942  nfsft_init_guru(&plan,NQ[iNQ],m_compare, NFSFT_NORMALIZED |
943  ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
944  ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
945  ((NQ[iNQ]>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
946  FFT_OUT_OF_PLACE, cutoff);
947  plan_ptr = &plan;
948  }
949  else
950  {
951  plan_ptr = &plan_adjoint;
952  }
953 
954  f_hat = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(NQ[iNQ])*sizeof(double _Complex));
955 
956  plan_adjoint_ptr->f_hat = f_hat;
957  plan_adjoint_ptr->x = x_grid;
958  plan_adjoint_ptr->f = f_grid;
959 
960  plan_ptr->f_hat = f_hat;
961  plan_ptr->x = x_compare;
962  plan_ptr->f = f;
963 
964  //fprintf(stderr,"Precomputing for x\n");
965  //fflush(stderr);
966  nfsft_precompute_x(plan_adjoint_ptr);
967  if (plan_adjoint_ptr != plan_ptr)
968  {
969  nfsft_precompute_x(plan_ptr);
970  }
971 
972  /* Initialize cumulative time variable. */
973  t_avg = 0.0;
974  err_infty_avg = 0.0;
975  err_2_avg = 0.0;
976 
977  /* Cycle through all runs. */
978  for (i = 0; i < 1/*repetitions*/; i++)
979  {
980  //fprintf(stderr,"Copying original values\n");
981  //fflush(stderr);
982  /* Copy exact funtion values to working array. */
983  //memcpy(f,f_grid,m_total*sizeof(double _Complex));
984 
985  /* Initialize time measurement. */
986  t0 = getticks();
987 
988  //fprintf(stderr,"Multiplying with quadrature weights\n");
989  //fflush(stderr);
990  /* Multiplication with the quadrature weights. */
991  /*fprintf(stderr,"\n");*/
992  d = 0;
993  for (k = 0; k < m_theta; k++)
994  {
995  for (n = 0; n < m_phi; n++)
996  {
997  /*fprintf(stderr,"f_ref[%d] = %le + I*%le,\t f[%d] = %le + I*%le, \t w[%d] = %le\n",
998  d,creal(f_ref[d]),cimag(f_ref[d]),d,creal(f[d]),cimag(f[d]),k,
999  w[k]);*/
1000  f_grid[d] *= w[k];
1001  d++;
1002  }
1003  }
1004 
1005  t1 = getticks();
1006  t_avg += nfft_elapsed_seconds(t1,t0);
1007 
1008  nfft_free(w);
1009 
1010  t0 = getticks();
1011 
1012  /*fprintf(stderr,"\n");
1013  d = 0;
1014  for (d = 0; d < grid_total; d++)
1015  {
1016  fprintf(stderr,"f[%d] = %le + I*%le, theta[%d] = %le, phi[%d] = %le\n",
1017  d,creal(f[d]),cimag(f[d]),d,x[2*d+1],d,x[2*d]);
1018  }*/
1019 
1020  //fprintf(stderr,"Executing adjoint\n");
1021  //fflush(stderr);
1022  /* Check if the fast NFSFT algorithm shall be tested. */
1023  if (use_nfsft != NO)
1024  {
1025  /* Execute the adjoint NFSFT transformation. */
1026  nfsft_adjoint(plan_adjoint_ptr);
1027  }
1028  else
1029  {
1030  /* Execute the adjoint direct NDSFT transformation. */
1031  nfsft_adjoint_direct(plan_adjoint_ptr);
1032  }
1033 
1034  /* Multiplication with the Fourier-Legendre coefficients. */
1035  /*for (k = 0; k <= m[im]; k++)
1036  {
1037  for (n = -k; n <= k; n++)
1038  {
1039  fprintf(stderr,"f_hat[%d,%d] = %le\t + I*%le\n",k,n,
1040  creal(f_hat[NFSFT_INDEX(k,n,&plan_adjoint)]),
1041  cimag(f_hat[NFSFT_INDEX(k,n,&plan_adjoint)]));
1042  }
1043  }*/
1044 
1045  //fprintf(stderr,"Executing trafo\n");
1046  //fflush(stderr);
1047  if (use_nfsft != NO)
1048  {
1049  /* Execute the NFSFT transformation. */
1050  nfsft_trafo(plan_ptr);
1051  }
1052  else
1053  {
1054  /* Execute the direct NDSFT transformation. */
1055  nfsft_trafo_direct(plan_ptr);
1056  }
1057 
1058  t1 = getticks();
1059  t_avg += nfft_elapsed_seconds(t1,t0);
1060 
1061  //fprintf(stderr,"Finalizing\n");
1062  //fflush(stderr);
1063  /* Finalize the NFSFT plans */
1064  nfsft_finalize(plan_adjoint_ptr);
1065  if (plan_ptr != plan_adjoint_ptr)
1066  {
1067  nfsft_finalize(plan_ptr);
1068  }
1069 
1070  /* Free data arrays. */
1071  nfft_free(f_hat);
1072  nfft_free(x_grid);
1073 
1074  err_infty_avg += X(error_l_infty_complex)(f, f_compare, m_compare);
1075  err_2_avg += X(error_l_2_complex)(f, f_compare, m_compare);
1076 
1077  nfft_free(f_grid);
1078 
1079  if (mode == RANDOM)
1080  {
1081  }
1082  else
1083  {
1084  nfft_free(f_compare);
1085  }
1086 
1087  /*for (d = 0; d < m_total; d++)
1088  {
1089  fprintf(stderr,"f_ref[%d] = %le + I*%le,\t f[%d] = %le + I*%le\n",
1090  d,creal(f_ref[d]),cimag(f_ref[d]),d,creal(f[d]),cimag(f[d]));
1091  }*/
1092  }
1093 
1094  //fprintf(stderr,"Calculating the error\n");
1095  //fflush(stderr);
1096  /* Calculate average time needed. */
1097  t_avg = t_avg/((double)repetitions);
1098 
1099  /* Calculate the average error. */
1100  err_infty_avg = err_infty_avg/((double)repetitions);
1101 
1102  /* Calculate the average error. */
1103  err_2_avg = err_2_avg/((double)repetitions);
1104 
1105  /* Print out the error measurements. */
1106  fprintf(stdout,"%+le %+le %+le\n", t_avg, err_infty_avg, err_2_avg);
1107  fprintf(stderr,"%d: %4d %4d %+le %+le %+le\n", tc, NQ[iNQ], SQ[iNQ],
1108  t_avg, err_infty_avg, err_2_avg);
1109  }
1110  } /* for (im = 0; im < im_max; im++) - Process all cut-off
1111  * bandwidths.*/
1112  fprintf(stderr,"\n");
1113 
1114  /* Delete precomputed data. */
1115  nfsft_forget();
1116 
1117  /* Free memory for cut-off bandwidths and grid size parameters. */
1118  nfft_free(NQ);
1119  nfft_free(SQ);
1120  if (testmode == TIMING)
1121  {
1122  nfft_free(RQ);
1123  }
1124 
1125  if (mode == RANDOM)
1126  {
1127  nfft_free(x_compare);
1128  nfft_free(f_compare);
1129  nfft_free(f);
1130  }
1131 
1132  if (testmode == TIMING)
1133  {
1134  /* Allocate data structures. */
1135  nfft_free(f_hat);
1136  nfft_free(f);
1137  nfft_free(x_grid);
1138  }
1139 
1140  } /* for (tc = 0; tc < tc_max; tc++) - Process each testcase. */
1141 
1142  /* Return exit code for successful run. */
1143  return EXIT_SUCCESS;
1144 }
1145 /* \} */
double * x
the nodes for ,
Definition: nfft3.h:576
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:576
R nfft_elapsed_seconds(ticks t1, ticks t0)
Return number of elapsed seconds between two time points.
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
Definition: nfft3.h:576
modes
TODO Add comment here.
Definition: quadratureS2.c:61
testtype
Enumeration for test modes.
Definition: quadratureS2.c:49
void nfft_free(void *p)
#define X(name)
Include header for C99 complex datatype.
Definition: fastsum.h:53
functiontype
Enumeration for test functions.
Definition: quadratureS2.c:56
void * nfft_malloc(size_t n)
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:576
fftw_complex * f
Samples.
Definition: nfft3.h:576
double threshold
The threshold /f$/f$.
gridtype
Enumeration for quadrature grid types.
Definition: quadratureS2.c:52