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