Subversion Repositories shark

Rev

Rev 3 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 1
/*
2
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
3
 *
4
 * This program is free software; you can redistribute it and/or modify
5
 * it under the terms of the GNU General Public License as published by
6
 * the Free Software Foundation; either version 2 of the License, or
7
 * (at your option) any later version.
8
 *
9
 * This program is distributed in the hope that it will be useful,
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12
 * GNU General Public License for more details.
13
 *
14
 * You should have received a copy of the GNU General Public License
15
 * along with this program; if not, write to the Free Software
16
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
17
 *
18
 */
19
 
20
/*
21
 * fftw_test.c : test program for complex-complex transforms
22
 */
23
 
24
/* $Id: rfftw_te.c,v 1.1.1.1 2002-03-29 14:12:59 pj Exp $ */
25
#include <fftw-int.h>
26
#include <stdlib.h>
27
#include <stdio.h>
28
#include <string.h>
29
#include <math.h>
30
#include <time.h>
31
 
32
#include <rfftw.h>
33
 
34
#include "test_main.h"
35
 
36
char fftw_prefix[] = "rfftw";
37
 
38
/*************************************************
39
 * Speed tests
40
 *************************************************/
41
 
42
void zero_arr(int n, fftw_real * a)
43
{
44
     int i;
45
     for (i = 0; i < n; ++i)
46
          a[i] = 0.0;
47
}
48
 
49
void test_speed_aux(int n, fftw_direction dir, int flags, int specific)
50
{
51
     fftw_real *in, *out;
52
     fftw_plan plan;
53
     double t;
54
     fftw_time begin, end;
55
 
56
     in = (fftw_real *) fftw_malloc(n * howmany_fields
57
                                    * sizeof(fftw_real));
58
     out = (fftw_real *) fftw_malloc(n * howmany_fields
59
                                     * sizeof(fftw_real));
60
 
61
     if (specific) {
62
          begin = fftw_get_time();
63
          plan = rfftw_create_plan_specific(n, dir,
64
                                        speed_flag | flags | wisdom_flag,
65
                                            in, howmany_fields,
66
                                            out, howmany_fields);
67
          end = fftw_get_time();
68
     } else {
69
          begin = fftw_get_time();
70
          plan = rfftw_create_plan(n, dir, speed_flag | flags | wisdom_flag);
71
          end = fftw_get_time();
72
     }
73
     CHECK(plan != NULL, "can't create plan");
74
 
75
     t = fftw_time_to_sec(fftw_time_diff(end, begin));
76
     WHEN_VERBOSE(2, printf("time for planner: %f s\n", t));
77
 
78
     WHEN_VERBOSE(2, rfftw_print_plan(plan));
79
 
80
     FFTW_TIME_FFT(rfftw(plan, howmany_fields,
81
                         in, howmany_fields, 1, out, howmany_fields, 1),
82
                   in, n * howmany_fields, t);
83
 
84
     rfftw_destroy_plan(plan);
85
 
86
     WHEN_VERBOSE(1, printf("time for one fft: %s", smart_sprint_time(t)));
87
     WHEN_VERBOSE(1, printf(" (%s/point)\n", smart_sprint_time(t / n)));
88
     WHEN_VERBOSE(1, printf("\"mflops\" = 5/2 (n log2 n) / (t in microseconds)"
89
                        " = %f\n", 0.5 * howmany_fields * mflops(t, n)));
90
 
91
     fftw_free(in);
92
     fftw_free(out);
93
 
94
     WHEN_VERBOSE(1, printf("\n"));
95
}
96
 
97
void test_speed_nd_aux(struct size sz,
98
                       fftw_direction dir, int flags, int specific)
99
{
100
     fftw_real *in;
101
     fftwnd_plan plan;
102
     double t;
103
     fftw_time begin, end;
104
     int i, N;
105
 
106
     /* only bench in-place multi-dim transforms */
107
     flags |= FFTW_IN_PLACE;   
108
 
109
     N = 1;
110
     for (i = 0; i < sz.rank - 1; ++i)
111
          N *= sz.narray[i];
112
 
113
     N *= (sz.narray[i] + 2);
114
 
115
     in = (fftw_real *) fftw_malloc(N * howmany_fields * sizeof(fftw_real));
116
 
117
     if (specific) {
118
          begin = fftw_get_time();
119
          plan = rfftwnd_create_plan_specific(sz.rank, sz.narray, dir,
120
                                        speed_flag | flags | wisdom_flag,
121
                                              in, howmany_fields, 0, 1);
122
     } else {
123
          begin = fftw_get_time();
124
          plan = rfftwnd_create_plan(sz.rank, sz.narray,
125
                                  dir, speed_flag | flags | wisdom_flag);
126
     }
127
     end = fftw_get_time();
128
     CHECK(plan != NULL, "can't create plan");
129
 
130
     t = fftw_time_to_sec(fftw_time_diff(end, begin));
131
     WHEN_VERBOSE(2, printf("time for planner: %f s\n", t));
132
 
133
     WHEN_VERBOSE(2, printf("\n"));
134
     WHEN_VERBOSE(2, (rfftwnd_print_plan(plan)));
135
     WHEN_VERBOSE(2, printf("\n"));
136
 
137
     if (dir == FFTW_REAL_TO_COMPLEX) {
138
          FFTW_TIME_FFT(rfftwnd_real_to_complex(plan, howmany_fields,
139
                                                in, howmany_fields, 1,
140
                                                0, 0, 0),
141
                        in, N * howmany_fields, t);
142
     } else {
143
          FFTW_TIME_FFT(rfftwnd_complex_to_real(plan, howmany_fields,
144
                                                (fftw_complex *) in,
145
                                                howmany_fields, 1,
146
                                                0, 0, 0),
147
                        in, N * howmany_fields, t);
148
     }
149
 
150
     rfftwnd_destroy_plan(plan);
151
 
152
     WHEN_VERBOSE(1, printf("time for one fft: %s", smart_sprint_time(t)));
153
     WHEN_VERBOSE(1, printf(" (%s/point)\n", smart_sprint_time(t / N)));
154
     WHEN_VERBOSE(1, printf("\"mflops\" = 5/2 (N log2 N) / (t in microseconds)"
155
                        " = %f\n", 0.5 * howmany_fields * mflops(t, N)));
156
 
157
     fftw_free(in);
158
 
159
     WHEN_VERBOSE(1, printf("\n"));
160
}
161
 
162
/*************************************************
163
 * correctness tests
164
 *************************************************/
165
 
166
void fill_random(fftw_real * a, int n, int stride)
167
{
168
     int i;
169
 
170
     /* generate random inputs */
171
     for (i = 0; i < n; ++i)
172
          a[i * stride] = DRAND();
173
}
174
 
175
double compute_error(fftw_real * A, int astride,
176
                     fftw_real * B, int bstride, int n)
177
{
178
     /* compute the relative error */
179
     double error = 0.0;
180
     int i;
181
 
182
     for (i = 0; i < n; ++i) {
183
          double a;
184
          double mag;
185
          a = fabs(A[i * astride] - B[i * bstride]);
186
          mag = 0.5 * (fabs(A[i * astride]) + fabs(B[i * bstride])) + TOLERANCE;
187
 
188
          a /= mag;
189
          if (a > error)
190
               error = a;
191
 
192
#ifdef HAVE_ISNAN
193
          CHECK(!isnan(a), "NaN in answer");
194
#endif
195
     }
196
     return error;
197
}
198
 
199
void array_compare(fftw_real * A, fftw_real * B, int n)
200
{
201
     CHECK(compute_error(A, 1, B, 1, n) < TOLERANCE,
202
           "failure in RFFTW verification");
203
}
204
 
205
void test_out_of_place(int n, int istride, int ostride,
206
                       int howmany, fftw_direction dir,
207
                       fftw_plan validated_plan, int specific)
208
{
209
     fftw_complex *in2, *out2;
210
     fftw_real *in1, *out1, *out3;
211
     fftw_plan plan;
212
     int i, j;
213
     int flags = measure_flag | wisdom_flag;
214
 
215
     if (coinflip())
216
          flags |= FFTW_THREADSAFE;
217
 
218
     in1 = (fftw_real *) fftw_malloc(istride * n * sizeof(fftw_real) * howmany);
219
     in2 = (fftw_complex *) fftw_malloc(n * sizeof(fftw_complex));
220
     out1 = (fftw_real *) fftw_malloc(ostride * n * sizeof(fftw_real) * howmany);
221
     out2 = (fftw_complex *) fftw_malloc(n * sizeof(fftw_complex));
222
     out3 = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
223
 
224
     if (!specific)
225
          plan = rfftw_create_plan(n, dir, flags);
226
     else
227
          plan = rfftw_create_plan_specific(n, dir, flags,
228
                                            in1, istride, out1, ostride);
229
     CHECK(plan != NULL, "can't create plan");
230
 
231
     /* generate random inputs */
232
     fill_random(in1, n, istride);
233
     for (j = 1; j < howmany; ++j)
234
          for (i = 0; i < n; ++i)
235
               in1[(j * n + i) * istride] = in1[i * istride];
236
 
237
     /* copy random inputs to complex array for comparison with fftw: */
238
     if (dir == FFTW_REAL_TO_COMPLEX)
239
          for (i = 0; i < n; ++i) {
240
               c_re(in2[i]) = in1[i * istride];
241
               c_im(in2[i]) = 0.0;
242
     } else {
243
          int n2 = (n + 1) / 2;
244
          c_re(in2[0]) = in1[0];
245
          c_im(in2[0]) = 0.0;
246
          for (i = 1; i < n2; ++i) {
247
               c_re(in2[i]) = in1[i * istride];
248
               c_im(in2[i]) = in1[(n - i) * istride];
249
          }
250
          if (n2 * 2 == n) {
251
               c_re(in2[n2]) = in1[n2 * istride];
252
               c_im(in2[n2]) = 0.0;
253
               ++i;
254
          }
255
          for (; i < n; ++i) {
256
               c_re(in2[i]) = c_re(in2[n - i]);
257
               c_im(in2[i]) = -c_im(in2[n - i]);
258
          }
259
     }
260
 
261
     /*
262
      * fill in other positions of the array, to make sure that
263
      * rfftw doesn't overwrite them
264
      */
265
     for (j = 1; j < istride; ++j)
266
          for (i = 0; i < n * howmany; ++i)
267
               in1[i * istride + j] = i * istride + j;
268
 
269
     for (j = 1; j < ostride; ++j)
270
          for (i = 0; i < n * howmany; ++i)
271
               out1[i * ostride + j] = -i * ostride + j;
272
 
273
     WHEN_VERBOSE(2, rfftw_print_plan(plan));
274
 
275
     /* fft-ize */
276
     if (howmany != 1 || istride != 1 || ostride != 1 || coinflip())
277
          rfftw(plan, howmany, in1, istride, n * istride, out1, ostride,
278
                n * ostride);
279
     else
280
          rfftw_one(plan, in1, out1);
281
 
282
     rfftw_destroy_plan(plan);
283
 
284
     /* check for overwriting */
285
     for (j = 1; j < istride; ++j)
286
          for (i = 0; i < n * howmany; ++i)
287
               CHECK(in1[i * istride + j] == i * istride + j,
288
                     "input has been overwritten");
289
     for (j = 1; j < ostride; ++j)
290
          for (i = 0; i < n * howmany; ++i)
291
               CHECK(out1[i * ostride + j] == -i * ostride + j,
292
                     "output has been overwritten");
293
 
294
     fftw(validated_plan, 1, in2, 1, n, out2, 1, n);
295
 
296
     if (dir == FFTW_REAL_TO_COMPLEX) {
297
          int n2 = (n + 1) / 2;
298
          out3[0] = c_re(out2[0]);
299
          for (i = 1; i < n2; ++i) {
300
               out3[i] = c_re(out2[i]);
301
               out3[n - i] = c_im(out2[i]);
302
          }
303
          if (n2 * 2 == n)
304
               out3[n2] = c_re(out2[n2]);
305
     } else {
306
          for (i = 0; i < n; ++i)
307
               out3[i] = c_re(out2[i]);
308
     }
309
 
310
     for (j = 0; j < howmany; ++j)
311
          CHECK(compute_error(out1 + j * n * ostride, ostride, out3, 1, n)
312
                < TOLERANCE,
313
                "test_out_of_place: wrong answer");
314
     WHEN_VERBOSE(2, printf("OK\n"));
315
 
316
     fftw_free(in1);
317
     fftw_free(in2);
318
     fftw_free(out1);
319
     fftw_free(out2);
320
     fftw_free(out3);
321
}
322
 
323
void test_in_place(int n, int istride,
324
                   int howmany, fftw_direction dir,
325
                   fftw_plan validated_plan, int specific)
326
{
327
     fftw_complex *in2, *out2;
328
     fftw_real *in1, *out1, *out3;
329
     fftw_plan plan;
330
     int i, j;
331
     int ostride = istride;
332
     int flags = measure_flag | wisdom_flag | FFTW_IN_PLACE;
333
 
334
     if (coinflip())
335
          flags |= FFTW_THREADSAFE;
336
 
337
     in1 = (fftw_real *) fftw_malloc(istride * n * sizeof(fftw_real) * howmany);
338
     in2 = (fftw_complex *) fftw_malloc(n * sizeof(fftw_complex));
339
     out1 = in1;
340
     out2 = (fftw_complex *) fftw_malloc(n * sizeof(fftw_complex));
341
     out3 = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
342
 
343
     if (!specific)
344
          plan = rfftw_create_plan(n, dir, flags);
345
     else
346
          plan = rfftw_create_plan_specific(n, dir, flags,
347
                                            in1, istride, out1, ostride);
348
     CHECK(plan != NULL, "can't create plan");
349
 
350
     /* generate random inputs */
351
     fill_random(in1, n, istride);
352
     for (j = 1; j < howmany; ++j)
353
          for (i = 0; i < n; ++i)
354
               in1[(j * n + i) * istride] = in1[i * istride];
355
 
356
     /* copy random inputs to complex array for comparison with fftw: */
357
     if (dir == FFTW_REAL_TO_COMPLEX)
358
          for (i = 0; i < n; ++i) {
359
               c_re(in2[i]) = in1[i * istride];
360
               c_im(in2[i]) = 0.0;
361
     } else {
362
          int n2 = (n + 1) / 2;
363
          c_re(in2[0]) = in1[0];
364
          c_im(in2[0]) = 0.0;
365
          for (i = 1; i < n2; ++i) {
366
               c_re(in2[i]) = in1[i * istride];
367
               c_im(in2[i]) = in1[(n - i) * istride];
368
          }
369
          if (n2 * 2 == n) {
370
               c_re(in2[n2]) = in1[n2 * istride];
371
               c_im(in2[n2]) = 0.0;
372
               ++i;
373
          }
374
          for (; i < n; ++i) {
375
               c_re(in2[i]) = c_re(in2[n - i]);
376
               c_im(in2[i]) = -c_im(in2[n - i]);
377
          }
378
     }
379
 
380
     /*
381
      * fill in other positions of the array, to make sure that
382
      * rfftw doesn't overwrite them
383
      */
384
     for (j = 1; j < istride; ++j)
385
          for (i = 0; i < n * howmany; ++i)
386
               in1[i * istride + j] = i * istride + j;
387
 
388
     WHEN_VERBOSE(2, rfftw_print_plan(plan));
389
 
390
     /* fft-ize */
391
     if (howmany != 1 || istride != 1 || coinflip())
392
          rfftw(plan, howmany, in1, istride, n * istride, 0, 0, 0);
393
     else
394
          rfftw_one(plan, in1, NULL);
395
 
396
     rfftw_destroy_plan(plan);
397
 
398
     /* check for overwriting */
399
     for (j = 1; j < ostride; ++j)
400
          for (i = 0; i < n * howmany; ++i)
401
               CHECK(out1[i * ostride + j] == i * ostride + j,
402
                     "output has been overwritten");
403
 
404
     fftw(validated_plan, 1, in2, 1, n, out2, 1, n);
405
 
406
     if (dir == FFTW_REAL_TO_COMPLEX) {
407
          int n2 = (n + 1) / 2;
408
          out3[0] = c_re(out2[0]);
409
          for (i = 1; i < n2; ++i) {
410
               out3[i] = c_re(out2[i]);
411
               out3[n - i] = c_im(out2[i]);
412
          }
413
          if (n2 * 2 == n)
414
               out3[n2] = c_re(out2[n2]);
415
     } else {
416
          for (i = 0; i < n; ++i)
417
               out3[i] = c_re(out2[i]);
418
     }
419
 
420
     for (j = 0; j < howmany; ++j)
421
          CHECK(compute_error(out1 + j * n * ostride, ostride, out3, 1, n)
422
                < TOLERANCE,
423
                "test_in_place: wrong answer");
424
     WHEN_VERBOSE(2, printf("OK\n"));
425
 
426
     fftw_free(in1);
427
     fftw_free(in2);
428
     fftw_free(out2);
429
     fftw_free(out3);
430
}
431
 
432
void test_out_of_place_both(int n, int istride, int ostride,
433
                            int howmany,
434
                            fftw_plan validated_plan_forward,
435
                            fftw_plan validated_plan_backward)
436
{
437
     int specific;
438
 
439
     for (specific = 0; specific <= 1; ++specific) {
440
          WHEN_VERBOSE(2,
441
               printf("TEST CORRECTNESS (out of place, FFTW_FORWARD, %s)"
442
                   " n = %d  istride = %d  ostride = %d  howmany = %d\n",
443
                      SPECIFICP(specific),
444
                      n, istride, ostride, howmany));
445
          test_out_of_place(n, istride, ostride, howmany, FFTW_FORWARD,
446
                            validated_plan_forward, specific);
447
 
448
          WHEN_VERBOSE(2,
449
              printf("TEST CORRECTNESS (out of place, FFTW_BACKWARD, %s)"
450
                   " n = %d  istride = %d  ostride = %d  howmany = %d\n",
451
                     SPECIFICP(specific),
452
                     n, istride, ostride, howmany));
453
          test_out_of_place(n, istride, ostride, howmany, FFTW_BACKWARD,
454
                            validated_plan_backward, specific);
455
     }
456
}
457
 
458
void test_in_place_both(int n, int istride, int howmany,
459
                        fftw_plan validated_plan_forward,
460
                        fftw_plan validated_plan_backward)
461
{
462
     int specific;
463
 
464
     for (specific = 0; specific <= 1; ++specific) {
465
          WHEN_VERBOSE(2,
466
                  printf("TEST CORRECTNESS (in place, FFTW_FORWARD, %s) "
467
                         "n = %d  istride = %d  howmany = %d\n",
468
                         SPECIFICP(specific),
469
                         n, istride, howmany));
470
          test_in_place(n, istride, howmany, FFTW_FORWARD,
471
                        validated_plan_forward, specific);
472
 
473
          WHEN_VERBOSE(2,
474
                 printf("TEST CORRECTNESS (in place, FFTW_BACKWARD, %s) "
475
                        "n = %d  istride = %d  howmany = %d\n",
476
                        SPECIFICP(specific),
477
                        n, istride, howmany));
478
          test_in_place(n, istride, howmany, FFTW_BACKWARD,
479
                        validated_plan_backward, specific);
480
     }
481
}
482
 
483
void test_correctness(int n)
484
{
485
     int istride, ostride, howmany;
486
     fftw_plan validated_plan_forward, validated_plan_backward;
487
 
488
     WHEN_VERBOSE(1,
489
                  printf("Testing correctness for n = %d...", n);
490
                  fflush(stdout));
491
 
492
     /* produce a *good* plan (validated by Ergun's test procedure) */
493
     validated_plan_forward =
494
         fftw_create_plan(n, FFTW_FORWARD, measure_flag | wisdom_flag);
495
     validated_plan_backward =
496
         fftw_create_plan(n, FFTW_BACKWARD, measure_flag | wisdom_flag);
497
     CHECK(validated_plan_forward != NULL, "can't create plan");
498
     CHECK(validated_plan_backward != NULL, "can't create plan");
499
 
500
     for (istride = 1; istride <= MAX_STRIDE; ++istride)
501
          for (ostride = 1; ostride <= MAX_STRIDE; ++ostride)
502
               for (howmany = 1; howmany <= MAX_HOWMANY; ++howmany)
503
                    test_out_of_place_both(n, istride, ostride, howmany,
504
                                           validated_plan_forward,
505
                                           validated_plan_backward);
506
 
507
     for (istride = 1; istride <= MAX_STRIDE; ++istride)
508
          for (howmany = 1; howmany <= MAX_HOWMANY; ++howmany)
509
               test_in_place_both(n, istride, howmany,
510
                                  validated_plan_forward,
511
                                  validated_plan_backward);
512
 
513
     fftw_destroy_plan(validated_plan_forward);
514
     fftw_destroy_plan(validated_plan_backward);
515
 
516
     if (!(wisdom_flag & FFTW_USE_WISDOM) && chk_mem_leak)
517
          fftw_check_memory_leaks();
518
 
519
     WHEN_VERBOSE(1, printf("OK\n"));
520
}
521
 
522
/*************************************************
523
 * multi-dimensional correctness tests
524
 *************************************************/
525
 
526
void testnd_out_of_place(int rank, int *n, fftwnd_plan validated_plan)
527
{
528
     int istride, ostride;
529
     int N, dim, i, j, k;
530
     int nc, nhc, nr;
531
     fftw_real *in1, *out3;
532
     fftw_complex *in2, *out1, *out2;
533
     fftwnd_plan p, ip;
534
     int flags = measure_flag | wisdom_flag;
535
 
536
     if (coinflip())
537
          flags |= FFTW_THREADSAFE;
538
 
539
     N = nc = nr = nhc = 1;
540
     for (dim = 0; dim < rank; ++dim)
541
          N *= n[dim];
542
     if (rank > 0) {
543
          nr = n[rank - 1];
544
          nc = N / nr;
545
          nhc = nr / 2 + 1;
546
     }
547
     in1 = (fftw_real *) fftw_malloc(N * MAX_STRIDE * sizeof(fftw_real));
548
     out3 = (fftw_real *) fftw_malloc(N * MAX_STRIDE * sizeof(fftw_real));
549
     out1 = (fftw_complex *) fftw_malloc(nhc * nc * MAX_STRIDE
550
                                         * sizeof(fftw_complex));
551
     in2 = (fftw_complex *) fftw_malloc(N * sizeof(fftw_complex));
552
     out2 = (fftw_complex *) fftw_malloc(N * sizeof(fftw_complex));
553
 
554
     p = rfftwnd_create_plan(rank, n, FFTW_REAL_TO_COMPLEX, flags);
555
     ip = rfftwnd_create_plan(rank, n, FFTW_COMPLEX_TO_REAL, flags);
556
     CHECK(p != NULL && ip != NULL, "can't create plan");
557
 
558
     for (istride = 1; istride <= MAX_STRIDE; ++istride) {
559
          /* generate random inputs */
560
          for (i = 0; i < nc; ++i)
561
               for (j = 0; j < nr; ++j) {
562
                    c_re(in2[i * nr + j]) = DRAND();
563
                    c_im(in2[i * nr + j]) = 0.0;
564
                    for (k = 0; k < istride; ++k)
565
                         in1[(i * nr + j) * istride + k]
566
                             = c_re(in2[i * nr + j]);
567
               }
568
          for (i = 0; i < N * istride; ++i)
569
               out3[i] = 0.0;
570
 
571
          fftwnd(validated_plan, 1, in2, 1, 1, out2, 1, 1);
572
 
573
          for (ostride = 1; ostride <= MAX_STRIDE; ++ostride) {
574
               int howmany = (istride < ostride) ? istride : ostride;
575
 
576
               WHEN_VERBOSE(2, printf("\n    testing stride %d/%d...",
577
                                      istride, ostride));
578
 
579
               if (howmany != 1 || istride != 1 || ostride != 1 || coinflip())
580
                    rfftwnd_real_to_complex(p, howmany, in1, istride, 1,
581
                                            out1, ostride, 1);
582
               else
583
                    rfftwnd_one_real_to_complex(p, in1, out1);
584
 
585
               for (i = 0; i < nc; ++i)
586
                    for (k = 0; k < howmany; ++k)
587
                         CHECK(compute_error_complex(out1 + i * nhc * ostride + k,
588
                                                     ostride,
589
                                                     out2 + i * nr, 1,
590
                                                     nhc) < TOLERANCE,
591
                               "out-of-place (r2c): wrong answer");
592
 
593
               if (howmany != 1 || istride != 1 || ostride != 1 || coinflip())
594
                    rfftwnd_complex_to_real(ip, howmany, out1, ostride, 1,
595
                                            out3, istride, 1);
596
               else
597
                    rfftwnd_one_complex_to_real(ip, out1, out3);
598
 
599
               for (i = 0; i < N * istride; ++i)
600
                    out3[i] *= 1.0 / N;
601
 
602
               if (istride == howmany)
603
                    CHECK(compute_error(out3, 1, in1, 1, N * istride)
604
                        < TOLERANCE, "out-of-place (c2r): wrong answer");
605
               for (i = 0; i < nc; ++i)
606
                    for (k = 0; k < howmany; ++k)
607
                         CHECK(compute_error(out3 + i * nr * istride + k,
608
                                             istride,
609
                                         (fftw_real *) (in2 + i * nr), 2,
610
                                             nr) < TOLERANCE,
611
                           "out-of-place (c2r): wrong answer (check 2)");
612
          }
613
     }
614
 
615
     rfftwnd_destroy_plan(p);
616
     rfftwnd_destroy_plan(ip);
617
 
618
     fftw_free(out3);
619
     fftw_free(out2);
620
     fftw_free(in2);
621
     fftw_free(out1);
622
     fftw_free(in1);
623
}
624
 
625
void testnd_in_place(int rank, int *n, fftwnd_plan validated_plan,
626
                     int alternate_api, int specific)
627
{
628
     int istride, ostride, howmany;
629
     int N, dim, i, j, k;
630
     int nc, nhc, nr;
631
     fftw_real *in1, *out3;
632
     fftw_complex *in2, *out1, *out2;
633
     fftwnd_plan p, ip;
634
     int flags = measure_flag | wisdom_flag | FFTW_IN_PLACE;
635
 
636
     if (coinflip())
637
          flags |= FFTW_THREADSAFE;
638
 
639
     N = nc = nr = nhc = 1;
640
     for (dim = 0; dim < rank; ++dim)
641
          N *= n[dim];
642
     if (rank > 0) {
643
          nr = n[rank - 1];
644
          nc = N / nr;
645
          nhc = nr / 2 + 1;
646
     }
647
     in1 = (fftw_real *) fftw_malloc(2 * nhc * nc * MAX_STRIDE * sizeof(fftw_real));
648
     out3 = in1;
649
     out1 = (fftw_complex *) in1;
650
     in2 = (fftw_complex *) fftw_malloc(N * sizeof(fftw_complex));
651
     out2 = (fftw_complex *) fftw_malloc(N * sizeof(fftw_complex));
652
 
653
     if (alternate_api && specific && (rank == 2 || rank == 3)) {
654
          if (rank == 2) {
655
               p = rfftw2d_create_plan_specific(n[0], n[1],
656
                                             FFTW_REAL_TO_COMPLEX, flags,
657
                                                in1, MAX_STRIDE, 0, 0);
658
               ip = rfftw2d_create_plan_specific(n[0], n[1],
659
                                             FFTW_COMPLEX_TO_REAL, flags,
660
                                                 in1, MAX_STRIDE, 0, 0);
661
          } else {
662
               p = rfftw3d_create_plan_specific(n[0], n[1], n[2],
663
                                             FFTW_REAL_TO_COMPLEX, flags,
664
                                                in1, MAX_STRIDE, 0, 0);
665
               ip = rfftw3d_create_plan_specific(n[0], n[1], n[2],
666
                                             FFTW_COMPLEX_TO_REAL, flags,
667
                                                 in1, MAX_STRIDE, 0, 0);
668
          }
669
     } else if (specific) {
670
          p = rfftwnd_create_plan_specific(rank, n, FFTW_REAL_TO_COMPLEX,
671
                                           flags,
672
                                       in1, MAX_STRIDE, in1, MAX_STRIDE);
673
          ip = rfftwnd_create_plan_specific(rank, n, FFTW_COMPLEX_TO_REAL,
674
                                            flags,
675
                                       in1, MAX_STRIDE, in1, MAX_STRIDE);
676
     } else if (alternate_api && (rank == 2 || rank == 3)) {
677
          if (rank == 2) {
678
               p = rfftw2d_create_plan(n[0], n[1], FFTW_REAL_TO_COMPLEX,
679
                                       flags);
680
               ip = rfftw2d_create_plan(n[0], n[1], FFTW_COMPLEX_TO_REAL,
681
                                        flags);
682
          } else {
683
               p = rfftw3d_create_plan(n[0], n[1], n[2], FFTW_REAL_TO_COMPLEX,
684
                                       flags);
685
               ip = rfftw3d_create_plan(n[0], n[1], n[2], FFTW_COMPLEX_TO_REAL,
686
                                        flags);
687
          }
688
     } else {
689
          p = rfftwnd_create_plan(rank, n, FFTW_REAL_TO_COMPLEX, flags);
690
          ip = rfftwnd_create_plan(rank, n, FFTW_COMPLEX_TO_REAL, flags);
691
     }
692
 
693
     CHECK(p != NULL && ip != NULL, "can't create plan");
694
 
695
     for (i = 0; i < nc * nhc * 2 * MAX_STRIDE; ++i)
696
          out3[i] = 0;
697
 
698
     for (istride = 1; istride <= MAX_STRIDE; ++istride) {
699
          /* generate random inputs */
700
          for (i = 0; i < nc; ++i)
701
               for (j = 0; j < nr; ++j) {
702
                    c_re(in2[i * nr + j]) = DRAND();
703
                    c_im(in2[i * nr + j]) = 0.0;
704
                    for (k = 0; k < istride; ++k)
705
                         in1[(i * nhc * 2 + j) * istride + k]
706
                             = c_re(in2[i * nr + j]);
707
               }
708
 
709
          fftwnd(validated_plan, 1, in2, 1, 1, out2, 1, 1);
710
 
711
          howmany = ostride = istride;
712
 
713
          WHEN_VERBOSE(2, printf("\n    testing in-place stride %d...",
714
                                 istride));
715
 
716
          if (howmany != 1 || istride != 1 || ostride != 1 || coinflip())
717
               rfftwnd_real_to_complex(p, howmany, in1, istride, 1,
718
                                       out1, ostride, 1);
719
          else
720
               rfftwnd_one_real_to_complex(p, in1, NULL);
721
 
722
          for (i = 0; i < nc; ++i)
723
               for (k = 0; k < howmany; ++k)
724
                    CHECK(compute_error_complex(out1 + i * nhc * ostride + k,
725
                                                ostride,
726
                                                out2 + i * nr, 1,
727
                                                nhc) < TOLERANCE,
728
                          "in-place (r2c): wrong answer");
729
 
730
          if (howmany != 1 || istride != 1 || ostride != 1 || coinflip())
731
               rfftwnd_complex_to_real(ip, howmany, out1, ostride, 1,
732
                                       out3, istride, 1);
733
          else
734
               rfftwnd_one_complex_to_real(ip, out1, NULL);
735
 
736
          for (i = 0; i < nc * nhc * 2 * istride; ++i)
737
               out3[i] *= 1.0 / N;
738
 
739
          for (i = 0; i < nc; ++i)
740
               for (k = 0; k < howmany; ++k)
741
                    CHECK(compute_error(out3 + i * nhc * 2 * istride + k,
742
                                        istride,
743
                                        (fftw_real *) (in2 + i * nr), 2,
744
                                        nr) < TOLERANCE,
745
                          "in-place (c2r): wrong answer (check 2)");
746
     }
747
 
748
     rfftwnd_destroy_plan(p);
749
     rfftwnd_destroy_plan(ip);
750
 
751
     fftw_free(out2);
752
     fftw_free(in2);
753
     fftw_free(in1);
754
}
755
 
756
void testnd_correctness(struct size sz, fftw_direction dir,
757
                        int alt_api, int specific, int force_buf)
758
{
759
     fftwnd_plan validated_plan;
760
 
761
     if (dir != FFTW_FORWARD)
762
          return;
763
     if (force_buf)
764
          return;
765
 
766
     validated_plan = fftwnd_create_plan(sz.rank, sz.narray,
767
                                         dir, measure_flag | wisdom_flag);
768
     CHECK(validated_plan != NULL, "can't create plan");
769
 
770
     testnd_out_of_place(sz.rank, sz.narray, validated_plan);
771
     testnd_in_place(sz.rank, sz.narray,
772
                     validated_plan, alt_api, specific);
773
 
774
     fftwnd_destroy_plan(validated_plan);
775
}
776
 
777
/*************************************************
778
 * planner tests
779
 *************************************************/
780
 
781
void test_planner(int rank)
782
{
783
     /*
784
      * create and destroy many plans, at random.  Check the
785
      * garbage-collecting allocator of twiddle factors
786
      */
787
     int i, dim;
788
     int r, s;
789
     fftw_plan p[PLANNER_TEST_SIZE];
790
     fftwnd_plan pnd[PLANNER_TEST_SIZE];
791
     int *narr, maxdim;
792
 
793
     chk_mem_leak = 0;
794
     verbose--;
795
 
796
     please_wait();
797
     if (rank < 1)
798
          rank = 1;
799
 
800
     narr = (int *) fftw_malloc(rank * sizeof(int));
801
 
802
     maxdim = (int) pow(8192.0, 1.0/rank);
803
 
804
     for (i = 0; i < PLANNER_TEST_SIZE; ++i) {
805
          p[i] = (fftw_plan) 0;
806
          pnd[i] = (fftwnd_plan) 0;
807
     }
808
 
809
     for (i = 0; i < PLANNER_TEST_SIZE * PLANNER_TEST_SIZE; ++i) {
810
          r = rand();
811
          if (r < 0)
812
               r = -r;
813
          r = r % PLANNER_TEST_SIZE;
814
 
815
          for (dim = 0; dim < rank; ++dim) {
816
               do {
817
                    s = rand();
818
                    if (s < 0)
819
                         s = -s;
820
                    s = s % maxdim + 1;
821
               } while (s == 0);
822
               narr[dim] = s;
823
          }
824
 
825
          if (rank == 1) {
826
               if (p[r])
827
                    rfftw_destroy_plan(p[r]);
828
 
829
               p[r] = rfftw_create_plan(narr[0], random_dir(), measure_flag |
830
                                        wisdom_flag);
831
               if (paranoid && narr[0] < 200)
832
                    test_correctness(narr[0]);
833
          }
834
          if (pnd[r])
835
               rfftwnd_destroy_plan(pnd[r]);
836
 
837
          pnd[r] = rfftwnd_create_plan(rank, narr,
838
                                       random_dir(), measure_flag |
839
                                       wisdom_flag);
840
 
841
          if (i % (PLANNER_TEST_SIZE * PLANNER_TEST_SIZE / 20) == 0) {
842
               WHEN_VERBOSE(0, printf("test planner: so far so good\n"));
843
               WHEN_VERBOSE(0, printf("test planner: iteration %d out of %d\n",
844
                              i, PLANNER_TEST_SIZE * PLANNER_TEST_SIZE));
845
          }
846
     }
847
 
848
     for (i = 0; i < PLANNER_TEST_SIZE; ++i) {
849
          if (p[i])
850
               rfftw_destroy_plan(p[i]);
851
          if (pnd[i])
852
               rfftwnd_destroy_plan(pnd[i]);
853
     }
854
 
855
     fftw_free(narr);
856
     verbose++;
857
     chk_mem_leak = 1;
858
}
859
 
860
 
861
/*************************************************
862
 * Ergun's test for real->complex transforms
863
 *************************************************/
864
static void rfill_random(fftw_real *a, int n)
865
{
866
     int i;
867
 
868
     for (i = 0; i < n; ++i) {
869
          a[i] = DRAND();
870
     }
871
}
872
 
873
static void rarray_copy(fftw_real *out, fftw_real *in, int n)
874
{
875
     int i;
876
 
877
     for (i = 0; i < n; ++i)
878
          out[i] = in[i];
879
}
880
 
881
/* C = A + B */
882
void rarray_add(fftw_real *C, fftw_real *A, fftw_real *B, int n)
883
{
884
     int i;
885
 
886
     for (i = 0; i < n; ++i) {
887
          C[i] = A[i] + B[i];
888
     }
889
}
890
 
891
/* C = A - B */
892
void rarray_sub(fftw_real *C, fftw_real *A, fftw_real *B, int n)
893
{
894
     int i;
895
 
896
     for (i = 0; i < n; ++i) {
897
          C[i] = A[i] - B[i];
898
     }
899
}
900
 
901
/* B = rotate left A */
902
void rarray_rol(fftw_real *B, fftw_real *A,
903
               int n, int n_before, int n_after)
904
{
905
     int i, ib, ia;
906
 
907
     for (ib = 0; ib < n_before; ++ib) {
908
          for (i = 0; i < n - 1; ++i)
909
               for (ia = 0; ia < n_after; ++ia)
910
                    B[(ib * n + i) * n_after + ia] =
911
                        A[(ib * n + i + 1) * n_after + ia];
912
 
913
          for (ia = 0; ia < n_after; ++ia)
914
               B[(ib * n + n - 1) * n_after + ia] = A[ib * n * n_after + ia];
915
     }
916
}
917
 
918
/* A = alpha * B  (out of place) */
919
void rarray_scale(fftw_real *A, fftw_real *B, fftw_real alpha, int n)
920
{
921
     int i;
922
 
923
     for (i = 0; i < n; ++i) {
924
          A[i] = B[i] * alpha;
925
     }
926
}
927
 
928
void rarray_compare(fftw_real *A, fftw_real *B, int n)
929
{
930
     double d = compute_error(A, 1, B, 1, n);
931
     if (d > TOLERANCE) {
932
          fflush(stdout);
933
          fprintf(stderr, "Found relative error %e\n", d);
934
          fftw_die("failure in Ergun's verification procedure\n");
935
     }
936
}
937
 
938
/*
939
 * guaranteed out-of-place transform.  Does the necessary
940
 * copying if the plan is in-place.
941
 */
942
static void rfftw_out_of_place(fftw_plan plan, int n,
943
                              fftw_real *in, fftw_real *out)
944
{
945
     if (plan->flags & FFTW_IN_PLACE) {
946
          rarray_copy(out, in, n);
947
          rfftw(plan, 1, out, 1, n, (fftw_real *)0, 1, n);
948
     } else {
949
          rfftw(plan, 1, in, 1, n, out, 1, n);
950
     }
951
}
952
 
953
/*
954
 * This is a real (as opposed to complex) variation of the FFT tester
955
 * described in
956
 *
957
 * Funda Ergün. Testing multivariate linear functions: Overcoming the
958
 * generator bottleneck. In Proceedings of the Twenty-Seventh Annual
959
 * ACM Symposium on the Theory of Computing, pages 407-416, Las Vegas,
960
 * Nevada, 29 May--1 June 1995.
961
 */
962
void test_ergun(int n, fftw_direction dir, fftw_plan plan)
963
{
964
     fftw_real *inA, *inB, *inC, *outA, *outB, *outC;
965
     fftw_real *inA1, *inB1;
966
     fftw_real *tmp;
967
     int i;
968
     int rounds = 20;
969
     FFTW_TRIG_REAL twopin = FFTW_K2PI / (FFTW_TRIG_REAL) n;
970
 
971
     inA = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
972
     inB = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
973
     inA1 = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
974
     inB1 = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
975
     inC = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
976
     outA = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
977
     outB = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
978
     outC = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
979
     tmp = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
980
 
981
     WHEN_VERBOSE(2,
982
                  printf("Validating plan, n = %d, dir = %s\n", n,
983
                         dir == FFTW_REAL_TO_COMPLEX ?
984
                         "REAL_TO_COMPLEX" : "COMPLEX_TO_REAL"));
985
 
986
     /* test 1: check linearity */
987
     for (i = 0; i < rounds; ++i) {
988
          fftw_real alpha, beta;
989
          alpha = DRAND();
990
          beta = DRAND();
991
          rfill_random(inA, n);
992
          rfill_random(inB, n);
993
          rarray_scale(inA1, inA, alpha, n);
994
          rarray_scale(inB1, inB, beta, n);
995
          rarray_add(inC, inA1, inB1, n);
996
          rfftw_out_of_place(plan, n, inA, outA);
997
          rfftw_out_of_place(plan, n, inB, outB);
998
          rarray_scale(outA, outA, alpha, n);
999
          rarray_scale(outB, outB, beta, n);
1000
          rarray_add(tmp, outA, outB, n);
1001
          rfftw_out_of_place(plan, n, inC, outC);
1002
          rarray_compare(outC, tmp, n);
1003
     }
1004
 
1005
     /* test 2: check that the unit impulse is transformed properly */
1006
     for (i = 0; i < n; ++i) {
1007
          /* impulse */
1008
          inA[i] = 0.0;
1009
 
1010
          /* transform of the impulse */
1011
          if (2 * i <= n)
1012
               outA[i] = 1.0;
1013
          else
1014
               outA[i] = 0.0;
1015
     }
1016
     inA[0] = 1.0;
1017
 
1018
     if (dir == FFTW_REAL_TO_COMPLEX) {
1019
          for (i = 0; i < rounds; ++i) {
1020
               rfill_random(inB, n);
1021
               rarray_sub(inC, inA, inB, n);
1022
               rfftw_out_of_place(plan, n, inB, outB);
1023
               rfftw_out_of_place(plan, n, inC, outC);
1024
               rarray_add(tmp, outB, outC, n);
1025
               rarray_compare(tmp, outA, n);
1026
          }
1027
     } else {
1028
          for (i = 0; i < rounds; ++i) {
1029
               rfill_random(outB, n);
1030
               rarray_sub(outC, outA, outB, n);
1031
               rfftw_out_of_place(plan, n, outB, inB);
1032
               rfftw_out_of_place(plan, n, outC, inC);
1033
               rarray_add(tmp, inB, inC, n);
1034
               rarray_scale(tmp, tmp, 1.0 / ((double) n), n);
1035
               rarray_compare(tmp, inA, n);
1036
          }
1037
     }
1038
 
1039
     /* test 3: check the time-shift property */
1040
     /* the paper performs more tests, but this code should be fine too */
1041
     if (dir == FFTW_REAL_TO_COMPLEX) {
1042
          for (i = 0; i < rounds; ++i) {
1043
               int j;
1044
 
1045
               rfill_random(inA, n);
1046
               rarray_rol(inB, inA, n, 1, 1);
1047
               rfftw_out_of_place(plan, n, inA, outA);
1048
               rfftw_out_of_place(plan, n, inB, outB);
1049
               tmp[0] = outB[0];
1050
               for (j = 1; 2 * j < n; ++j) {
1051
                    FFTW_TRIG_REAL s = dir * FFTW_TRIG_SIN(j * twopin);
1052
                    FFTW_TRIG_REAL c = FFTW_TRIG_COS(j * twopin);
1053
                    tmp[j] = outB[j] * c - outB[n - j] * s;
1054
                    tmp[n - j] = outB[j] * s + outB[n - j] * c;
1055
               }
1056
               if (2 * j == n)
1057
                    tmp[j] = -outB[j];
1058
 
1059
               rarray_compare(tmp, outA, n);
1060
          }
1061
     } else {
1062
          for (i = 0; i < rounds; ++i) {
1063
               int j;
1064
 
1065
               rfill_random(inA, n);
1066
               inB[0] = inA[0];
1067
               for (j = 1; 2 * j < n; ++j) {
1068
                    FFTW_TRIG_REAL s = dir * FFTW_TRIG_SIN(j * twopin);
1069
                    FFTW_TRIG_REAL c = FFTW_TRIG_COS(j * twopin);
1070
                    inB[j] = inA[j] * c - inA[n - j] * s;
1071
                    inB[n - j] = inA[j] * s + inA[n - j] * c;
1072
               }
1073
               if (2 * j == n)
1074
                    inB[j] = -inA[j];
1075
 
1076
               rfftw_out_of_place(plan, n, inA, outA);
1077
               rfftw_out_of_place(plan, n, inB, outB);         
1078
               rarray_rol(tmp, outA, n, 1, 1);
1079
               rarray_compare(tmp, outB, n);
1080
          }
1081
     }
1082
 
1083
     WHEN_VERBOSE(2, printf("Validation done\n"));
1084
 
1085
     fftw_free(tmp);
1086
     fftw_free(outC);
1087
     fftw_free(outB);
1088
     fftw_free(outA);
1089
     fftw_free(inC);
1090
     fftw_free(inB1);
1091
     fftw_free(inA1);
1092
     fftw_free(inB);
1093
     fftw_free(inA);
1094
}
1095
 
1096
static void rfftw_plan_hook_function(fftw_plan p)
1097
{
1098
     WHEN_VERBOSE(3, printf("Validating tentative plan\n"));
1099
     WHEN_VERBOSE(3, fftw_print_plan(p));
1100
     test_ergun(p->n, p->dir, p);
1101
}
1102
 
1103
/*************************************************
1104
 * test initialization
1105
 *************************************************/
1106
 
1107
void test_init(int *argc, char **argv)
1108
{
1109
}
1110
 
1111
void test_finish(void)
1112
{
1113
}
1114
 
1115
extern void (*rfftw_plan_hook)(fftw_plan plan);
1116
void enter_paranoid_mode(void)
1117
{
1118
     rfftw_plan_hook = rfftw_plan_hook_function;
1119
}
1120
 
1121
int get_option(int argc, char **argv, char *argval, int argval_maxlen)
1122
{
1123
     return default_get_option(argc,argv,argval,argval_maxlen);
1124
}