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