← Back to branch summary

~diresu/blender/blender-command-port

« back to all changes in this revision

Viewing changes to extern/fftw/libbench2/verify-r2r.c

  • Committer: theeth
  • Date: 2008-10-14 16:52:04 UTC
  • Revision ID: vcs-imports@canonical.com-20081014165204-r32w2gm6s0osvdhn
copy back trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/fftw/libbench2/verify-r2r.c
 
1
/*
 
2
 * Copyright (c) 2003, 2006 Matteo Frigo
 
3
 * Copyright (c) 2003, 2006 Massachusetts Institute of Technology
 
4
 *
 
5
 * This program is free software; you can redistribute it and/or modify
 
6
 * it under the terms of the GNU General Public License as published by
 
7
 * the Free Software Foundation; either version 2 of the License, or
 
8
 * (at your option) any later version.
 
9
 *
 
10
 * This program is distributed in the hope that it will be useful,
 
11
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
 * GNU General Public License for more details.
 
14
 *
 
15
 * You should have received a copy of the GNU General Public License
 
16
 * along with this program; if not, write to the Free Software
 
17
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
18
 *
 
19
 */
 
20
 
 
21
/* Lots of ugly duplication from verify-lib.c, plus lots of ugliness in
 
22
   general for all of the r2r variants...oh well, for now */
 
23
 
 
24
#include "verify.h"
 
25
#include <math.h>
 
26
#include <stdlib.h>
 
27
#include <stdio.h>
 
28
 
 
29
typedef struct {
 
30
     bench_problem *p;
 
31
     bench_tensor *probsz;
 
32
     bench_tensor *totalsz;
 
33
     bench_tensor *pckdsz;
 
34
     bench_tensor *pckdvecsz;
 
35
} info;
 
36
 
 
37
/*
 
38
 * Utility functions:
 
39
 */
 
40
 
 
41
static double dabs(double x) { return (x < 0.0) ? -x : x; }
 
42
static double dmin(double x, double y) { return (x < y) ? x : y; }
 
43
 
 
44
static double raerror(R *a, R *b, int n)
 
45
{
 
46
     if (n > 0) {
 
47
          /* compute the relative Linf error */
 
48
          double e = 0.0, mag = 0.0;
 
49
          int i;
 
50
 
 
51
          for (i = 0; i < n; ++i) {
 
52
               e = dmax(e, dabs(a[i] - b[i]));
 
53
               mag = dmax(mag, dmin(dabs(a[i]), dabs(b[i])));
 
54
          }
 
55
          if (dabs(mag) < 1e-14 && dabs(e) < 1e-14)
 
56
               e = 0.0;
 
57
          else
 
58
               e /= mag;
 
59
 
 
60
#ifdef HAVE_ISNAN
 
61
          BENCH_ASSERT(!isnan(e));
 
62
#endif
 
63
          return e;
 
64
     } else
 
65
          return 0.0;
 
66
}
 
67
 
 
68
#define by2pi(m, n) ((K2PI * (m)) / (n))
 
69
 
 
70
/*
 
71
 * Improve accuracy by reducing x to range [0..1/8]
 
72
 * before multiplication by 2 * PI.
 
73
 */
 
74
 
 
75
static trigreal bench_sincos(trigreal m, trigreal n, int sinp)
 
76
{
 
77
     /* waiting for C to get tail recursion... */
 
78
     trigreal half_n = n * 0.5;
 
79
     trigreal quarter_n = half_n * 0.5;
 
80
     trigreal eighth_n = quarter_n * 0.5;
 
81
     trigreal sgn = 1.0;
 
82
 
 
83
     if (sinp) goto sin;
 
84
 cos:
 
85
     if (m < 0) { m = -m; /* goto cos; */ }
 
86
     if (m > half_n) { m = n - m; goto cos; }
 
87
     if (m > eighth_n) { m = quarter_n - m; goto sin; }
 
88
     return sgn * COS(by2pi(m, n));
 
89
 
 
90
 msin:
 
91
     sgn = -sgn;
 
92
 sin:
 
93
     if (m < 0) { m = -m; goto msin; }
 
94
     if (m > half_n) { m = n - m; goto msin; }
 
95
     if (m > eighth_n) { m = quarter_n - m; goto cos; }
 
96
     return sgn * SIN(by2pi(m, n));
 
97
}
 
98
 
 
99
static trigreal cos2pi(int m, int n)
 
100
{
 
101
     return bench_sincos((trigreal)m, (trigreal)n, 0);
 
102
}
 
103
 
 
104
static trigreal sin2pi(int m, int n)
 
105
{
 
106
     return bench_sincos((trigreal)m, (trigreal)n, 1);
 
107
}
 
108
 
 
109
static trigreal cos00(int i, int j, int n)
 
110
{
 
111
     return cos2pi(i * j, n);
 
112
}
 
113
 
 
114
static trigreal cos01(int i, int j, int n)
 
115
{
 
116
     return cos00(i, 2*j + 1, 2*n);
 
117
}
 
118
 
 
119
static trigreal cos10(int i, int j, int n)
 
120
{
 
121
     return cos00(2*i + 1, j, 2*n);
 
122
}
 
123
 
 
124
static trigreal cos11(int i, int j, int n)
 
125
{
 
126
     return cos00(2*i + 1, 2*j + 1, 4*n);
 
127
}
 
128
 
 
129
static trigreal sin00(int i, int j, int n)
 
130
{
 
131
     return sin2pi(i * j, n);
 
132
}
 
133
 
 
134
static trigreal sin01(int i, int j, int n)
 
135
{
 
136
     return sin00(i, 2*j + 1, 2*n);
 
137
}
 
138
 
 
139
static trigreal sin10(int i, int j, int n)
 
140
{
 
141
     return sin00(2*i + 1, j, 2*n);
 
142
}
 
143
 
 
144
static trigreal sin11(int i, int j, int n)
 
145
{
 
146
     return sin00(2*i + 1, 2*j + 1, 4*n);
 
147
}
 
148
 
 
149
static trigreal realhalf(int i, int j, int n)
 
150
{
 
151
     UNUSED(i);
 
152
     if (j <= n - j)
 
153
          return 1.0;
 
154
     else
 
155
          return 0.0;
 
156
}
 
157
 
 
158
static trigreal coshalf(int i, int j, int n)
 
159
{
 
160
     if (j <= n - j)
 
161
          return cos00(i, j, n);
 
162
     else
 
163
          return cos00(i, n - j, n);
 
164
}
 
165
 
 
166
static trigreal unity(int i, int j, int n)
 
167
{
 
168
     UNUSED(i);
 
169
     UNUSED(j);
 
170
     UNUSED(n);
 
171
     return 1.0;
 
172
}
 
173
 
 
174
typedef trigreal (*trigfun)(int, int, int);
 
175
 
 
176
static void rarand(R *a, int n)
 
177
{
 
178
     int i;
 
179
 
 
180
     /* generate random inputs */
 
181
     for (i = 0; i < n; ++i) {
 
182
          a[i] = mydrand();
 
183
     }
 
184
}
 
185
 
 
186
/* C = A + B */
 
187
static void raadd(R *c, R *a, R *b, int n)
 
188
{
 
189
     int i;
 
190
 
 
191
     for (i = 0; i < n; ++i) {
 
192
          c[i] = a[i] + b[i];
 
193
     }
 
194
}
 
195
 
 
196
/* C = A - B */
 
197
static void rasub(R *c, R *a, R *b, int n)
 
198
{
 
199
     int i;
 
200
 
 
201
     for (i = 0; i < n; ++i) {
 
202
          c[i] = a[i] - b[i];
 
203
     }
 
204
}
 
205
 
 
206
/* B = rotate left A + rotate right A */
 
207
static void rarolr(R *b, R *a, int n, int nb, int na, 
 
208
                   r2r_kind_t k)
 
209
{
 
210
     int isL0 = 0, isL1 = 0, isR0 = 0, isR1 = 0;
 
211
     int i, ib, ia;
 
212
 
 
213
     for (ib = 0; ib < nb; ++ib) {
 
214
          for (i = 0; i < n - 1; ++i)
 
215
               for (ia = 0; ia < na; ++ia)
 
216
                    b[(ib * n + i) * na + ia] =
 
217
                         a[(ib * n + i + 1) * na + ia];
 
218
 
 
219
          /* ugly switch to do boundary conditions for various r2r types */
 
220
          switch (k) {
 
221
               /* periodic boundaries */
 
222
              case R2R_DHT:
 
223
              case R2R_R2HC:
 
224
                   for (ia = 0; ia < na; ++ia) {
 
225
                        b[(ib * n + n - 1) * na + ia] = 
 
226
                             a[(ib * n + 0) * na + ia];
 
227
                        b[(ib * n + 0) * na + ia] += 
 
228
                             a[(ib * n + n - 1) * na + ia];
 
229
                   }
 
230
                   break;
 
231
                   
 
232
              case R2R_HC2R: /* ugh (hermitian halfcomplex boundaries) */
 
233
                   if (n > 2) {
 
234
                        if (n % 2 == 0)
 
235
                             for (ia = 0; ia < na; ++ia) {
 
236
                                  b[(ib * n + n - 1) * na + ia] = 0.0;
 
237
                                  b[(ib * n + 0) * na + ia] += 
 
238
                                       a[(ib * n + 1) * na + ia];
 
239
                                  b[(ib * n + n/2) * na + ia] += 
 
240
                                       + a[(ib * n + n/2 - 1) * na + ia]
 
241
                                       - a[(ib * n + n/2 + 1) * na + ia];
 
242
                                  b[(ib * n + n/2 + 1) * na + ia] += 
 
243
                                       - a[(ib * n + n/2) * na + ia];
 
244
                             }
 
245
                        else 
 
246
                             for (ia = 0; ia < na; ++ia) {
 
247
                                  b[(ib * n + n - 1) * na + ia] = 0.0;
 
248
                                  b[(ib * n + 0) * na + ia] += 
 
249
                                       a[(ib * n + 1) * na + ia];
 
250
                                  b[(ib * n + n/2) * na + ia] += 
 
251
                                       + a[(ib * n + n/2) * na + ia]
 
252
                                       - a[(ib * n + n/2 + 1) * na + ia];
 
253
                                  b[(ib * n + n/2 + 1) * na + ia] += 
 
254
                                       - a[(ib * n + n/2 + 1) * na + ia]
 
255
                                       - a[(ib * n + n/2) * na + ia];
 
256
                             }
 
257
                   } else /* n <= 2 */ {
 
258
                        for (ia = 0; ia < na; ++ia) {
 
259
                             b[(ib * n + n - 1) * na + ia] =
 
260
                                  a[(ib * n + 0) * na + ia];
 
261
                             b[(ib * n + 0) * na + ia] += 
 
262
                                  a[(ib * n + n - 1) * na + ia];
 
263
                        }
 
264
                   }
 
265
                   break;
 
266
                   
 
267
              /* various even/odd boundary conditions */
 
268
              case R2R_REDFT00:
 
269
                   isL1 = isR1 = 1;
 
270
                   goto mirrors;
 
271
              case R2R_REDFT01:
 
272
                   isL1 = 1;
 
273
                   goto mirrors;
 
274
              case R2R_REDFT10:
 
275
                   isL0 = isR0 = 1;
 
276
                   goto mirrors;
 
277
              case R2R_REDFT11:
 
278
                   isL0 = 1;
 
279
                   isR0 = -1;
 
280
                   goto mirrors;
 
281
              case R2R_RODFT00:
 
282
                   goto mirrors;
 
283
              case R2R_RODFT01:
 
284
                   isR1 = 1;
 
285
                   goto mirrors;
 
286
              case R2R_RODFT10:
 
287
                   isL0 = isR0 = -1;
 
288
                   goto mirrors;
 
289
              case R2R_RODFT11:
 
290
                   isL0 = -1;
 
291
                   isR0 = 1;
 
292
                   goto mirrors;
 
293
 
 
294
          mirrors:
 
295
                   
 
296
                   for (ia = 0; ia < na; ++ia)
 
297
                        b[(ib * n + n - 1) * na + ia] = 
 
298
                             isR0 * a[(ib * n + n - 1) * na + ia]
 
299
                             + (n > 1 ? isR1 * a[(ib * n + n - 2) * na + ia]
 
300
                                : 0);
 
301
                   
 
302
                   for (ia = 0; ia < na; ++ia)
 
303
                        b[(ib * n) * na + ia] += 
 
304
                             isL0 * a[(ib * n) * na + ia]
 
305
                             + (n > 1 ? isL1 * a[(ib * n + 1) * na + ia] : 0);
 
306
                   
 
307
          }
 
308
 
 
309
          for (i = 1; i < n; ++i)
 
310
               for (ia = 0; ia < na; ++ia)
 
311
                    b[(ib * n + i) * na + ia] +=
 
312
                         a[(ib * n + i - 1) * na + ia];
 
313
     }
 
314
}
 
315
 
 
316
static void raphase_shift(R *b, R *a, int n, int nb, int na,
 
317
                         int n0, int k0, trigfun t)
 
318
{
 
319
     int j, jb, ja;
 
320
 
 
321
     for (jb = 0; jb < nb; ++jb)
 
322
          for (j = 0; j < n; ++j) {
 
323
               trigreal c = 2.0 * t(1, j + k0, n0);
 
324
 
 
325
               for (ja = 0; ja < na; ++ja) {
 
326
                    int k = (jb * n + j) * na + ja;
 
327
                    b[k] = a[k] * c;
 
328
               }
 
329
          }
 
330
}
 
331
 
 
332
/* A = alpha * A  (real, in place) */
 
333
static void rascale(R *a, R alpha, int n)
 
334
{
 
335
     int i;
 
336
 
 
337
     for (i = 0; i < n; ++i) {
 
338
          a[i] *= alpha;
 
339
     }
 
340
}
 
341
 
 
342
/*
 
343
 * compute rdft:
 
344
 */
 
345
 
 
346
/* copy real A into real B, using output stride of A and input stride of B */
 
347
typedef struct {
 
348
     dotens2_closure k;
 
349
     R *ra;
 
350
     R *rb;
 
351
} cpyr_closure;
 
352
 
 
353
static void cpyr0(dotens2_closure *k_, 
 
354
                  int indxa, int ondxa, int indxb, int ondxb)
 
355
{
 
356
     cpyr_closure *k = (cpyr_closure *)k_;
 
357
     k->rb[indxb] = k->ra[ondxa];
 
358
     UNUSED(indxa); UNUSED(ondxb);
 
359
}
 
360
 
 
361
static void cpyr(R *ra, bench_tensor *sza, R *rb, bench_tensor *szb)
 
362
{
 
363
     cpyr_closure k;
 
364
     k.k.apply = cpyr0;
 
365
     k.ra = ra; k.rb = rb;
 
366
     bench_dotens2(sza, szb, &k.k);
 
367
}
 
368
 
 
369
static void dofft(info *nfo, R *in, R *out)
 
370
{
 
371
     cpyr(in, nfo->pckdsz, (R *) nfo->p->in, nfo->totalsz);
 
372
     doit(1, nfo->p);
 
373
     cpyr((R *) nfo->p->out, nfo->totalsz, out, nfo->pckdsz);
 
374
}
 
375
 
 
376
static double racmp(R *a, R *b, int n, const char *test, double tol)
 
377
{
 
378
     double d = raerror(a, b, n);
 
379
     if (d > tol) {
 
380
          fprintf(stderr, "Found relative error %e (%s)\n", d, test);
 
381
          {
 
382
               int i;
 
383
               for (i = 0; i < n; ++i)
 
384
                    fprintf(stderr, "%8d %16.12f   %16.12f\n", i, 
 
385
                            (double) a[i],
 
386
                            (double) b[i]);
 
387
          }
 
388
          exit(EXIT_FAILURE);
 
389
     }
 
390
     return d;
 
391
}
 
392
 
 
393
/***********************************************************************/
 
394
 
 
395
typedef struct {
 
396
     int n; /* physical size */
 
397
     int n0; /* "logical" transform size */
 
398
     int i0, k0; /* shifts of input/output */
 
399
     trigfun ti, ts;  /* impulse/shift trig functions */
 
400
} dim_stuff;
 
401
 
 
402
static void impulse_response(int rnk, dim_stuff *d, R impulse_amp,
 
403
                             R *A, int N)
 
404
{
 
405
     if (rnk == 0)
 
406
          A[0] = impulse_amp;
 
407
     else {
 
408
          int i;
 
409
          N /= d->n;
 
410
          for (i = 0; i < d->n; ++i) {
 
411
               impulse_response(rnk - 1, d + 1,
 
412
                                impulse_amp * d->ti(d->i0, d->k0 + i, d->n0),
 
413
                                A + i * N, N);
 
414
          }
 
415
     }
 
416
}
 
417
 
 
418
/***************************************************************************/
 
419
 
 
420
/*
 
421
 * Implementation of the FFT tester described in
 
422
 *
 
423
 * Funda Erg�n. Testing multivariate linear functions: Overcoming the
 
424
 * generator bottleneck. In Proceedings of the Twenty-Seventh Annual
 
425
 * ACM Symposium on the Theory of Computing, pages 407-416, Las Vegas,
 
426
 * Nevada, 29 May--1 June 1995.
 
427
 *
 
428
 * Also: F. Ergun, S. R. Kumar, and D. Sivakumar, "Self-testing without
 
429
 * the generator bottleneck," SIAM J. on Computing 29 (5), 1630-51 (2000).
 
430
 */
 
431
 
 
432
static double rlinear(int n, info *nfo, R *inA, R *inB, R *inC, R *outA,
 
433
                      R *outB, R *outC, R *tmp, int rounds, double tol)
 
434
{
 
435
     double e = 0.0;
 
436
     int j;
 
437
 
 
438
     for (j = 0; j < rounds; ++j) {
 
439
          R alpha, beta;
 
440
          alpha = mydrand();
 
441
          beta = mydrand();
 
442
          rarand(inA, n);
 
443
          rarand(inB, n);
 
444
          dofft(nfo, inA, outA);
 
445
          dofft(nfo, inB, outB);
 
446
 
 
447
          rascale(outA, alpha, n);
 
448
          rascale(outB, beta, n);
 
449
          raadd(tmp, outA, outB, n);
 
450
          rascale(inA, alpha, n);
 
451
          rascale(inB, beta, n);
 
452
          raadd(inC, inA, inB, n);
 
453
          dofft(nfo, inC, outC);
 
454
 
 
455
          e = dmax(e, racmp(outC, tmp, n, "linear", tol));
 
456
     }
 
457
     return e;
 
458
}
 
459
 
 
460
static double rimpulse(dim_stuff *d, R impulse_amp,
 
461
                       int n, int vecn, info *nfo, 
 
462
                       R *inA, R *inB, R *inC,
 
463
                       R *outA, R *outB, R *outC,
 
464
                       R *tmp, int rounds, double tol)
 
465
{
 
466
     double e = 0.0;
 
467
     int N = n * vecn;
 
468
     int i;
 
469
     int j;
 
470
 
 
471
     /* test 2: check that the unit impulse is transformed properly */
 
472
 
 
473
     for (i = 0; i < N; ++i) {
 
474
          /* pls */
 
475
          inA[i] = 0.0;
 
476
     }
 
477
     for (i = 0; i < vecn; ++i) {
 
478
          inA[i * n] = (i+1) / (double)(vecn+1);
 
479
     
 
480
          /* transform of the pls */
 
481
          impulse_response(nfo->probsz->rnk, d, impulse_amp * inA[i * n],
 
482
                           outA + i * n, n);
 
483
     }
 
484
 
 
485
     dofft(nfo, inA, tmp);
 
486
     e = dmax(e, racmp(tmp, outA, N, "impulse 1", tol));
 
487
 
 
488
     for (j = 0; j < rounds; ++j) {
 
489
          rarand(inB, N);
 
490
          rasub(inC, inA, inB, N);
 
491
          dofft(nfo, inB, outB);
 
492
          dofft(nfo, inC, outC);
 
493
          raadd(tmp, outB, outC, N);
 
494
          e = dmax(e, racmp(tmp, outA, N, "impulse", tol));
 
495
     }
 
496
     return e;
 
497
}
 
498
 
 
499
static double t_shift(int n, int vecn, info *nfo, 
 
500
                      R *inA, R *inB, R *outA, R *outB, R *tmp,
 
501
                      int rounds, double tol,
 
502
                      dim_stuff *d)
 
503
{
 
504
     double e = 0.0;
 
505
     int nb, na, dim, N = n * vecn;
 
506
     int i, j;
 
507
     bench_tensor *sz = nfo->probsz;
 
508
 
 
509
     /* test 3: check the time-shift property */
 
510
     /* the paper performs more tests, but this code should be fine too */
 
511
 
 
512
     nb = 1;
 
513
     na = n;
 
514
 
 
515
     /* check shifts across all SZ dimensions */
 
516
     for (dim = 0; dim < sz->rnk; ++dim) {
 
517
          int ncur = sz->dims[dim].n;
 
518
 
 
519
          na /= ncur;
 
520
 
 
521
          for (j = 0; j < rounds; ++j) {
 
522
               rarand(inA, N);
 
523
 
 
524
               for (i = 0; i < vecn; ++i) {
 
525
                    rarolr(inB + i * n, inA + i*n, ncur, nb,na, 
 
526
                          nfo->p->k[dim]);
 
527
               }
 
528
               dofft(nfo, inA, outA);
 
529
               dofft(nfo, inB, outB);
 
530
               for (i = 0; i < vecn; ++i) 
 
531
                    raphase_shift(tmp + i * n, outA + i * n, ncur, 
 
532
                                 nb, na, d[dim].n0, d[dim].k0, d[dim].ts);
 
533
               e = dmax(e, racmp(tmp, outB, N, "time shift", tol));
 
534
          }
 
535
 
 
536
          nb *= ncur;
 
537
     }
 
538
     return e;
 
539
}
 
540
 
 
541
/***********************************************************************/
 
542
 
 
543
void verify_r2r(bench_problem *p, int rounds, double tol, errors *e)
 
544
{
 
545
     R *inA, *inB, *inC, *outA, *outB, *outC, *tmp;
 
546
     info nfo;
 
547
     int n, vecn, N;
 
548
     double impulse_amp = 1.0;
 
549
     dim_stuff *d;
 
550
     int i;
 
551
 
 
552
     if (rounds == 0)
 
553
          rounds = 20;  /* default value */
 
554
 
 
555
     n = tensor_sz(p->sz);
 
556
     vecn = tensor_sz(p->vecsz);
 
557
     N = n * vecn;
 
558
 
 
559
     d = (dim_stuff *) bench_malloc(sizeof(dim_stuff) * p->sz->rnk);
 
560
     for (i = 0; i < p->sz->rnk; ++i) {
 
561
          int n0, i0, k0;
 
562
          trigfun ti, ts;
 
563
 
 
564
          d[i].n = n0 = p->sz->dims[i].n;
 
565
          if (p->k[i] > R2R_DHT)
 
566
               n0 = 2 * (n0 + (p->k[i] == R2R_REDFT00 ? -1 : 
 
567
                               (p->k[i] == R2R_RODFT00 ? 1 : 0)));
 
568
          
 
569
          switch (p->k[i]) {
 
570
              case R2R_R2HC:
 
571
                   i0 = k0 = 0;
 
572
                   ti = realhalf;
 
573
                   ts = coshalf;
 
574
                   break;
 
575
              case R2R_DHT:
 
576
                   i0 = k0 = 0;
 
577
                   ti = unity;
 
578
                   ts = cos00;
 
579
                   break;
 
580
              case R2R_HC2R:
 
581
                   i0 = k0 = 0;
 
582
                   ti = unity;
 
583
                   ts = cos00;
 
584
                   break;
 
585
              case R2R_REDFT00:
 
586
                   i0 = k0 = 0;
 
587
                   ti = ts = cos00;
 
588
                   break;
 
589
              case R2R_REDFT01:
 
590
                   i0 = k0 = 0;
 
591
                   ti = ts = cos01;
 
592
                   break;
 
593
              case R2R_REDFT10:
 
594
                   i0 = k0 = 0;
 
595
                   ti = cos10; impulse_amp *= 2.0;
 
596
                   ts = cos00;
 
597
                   break;
 
598
              case R2R_REDFT11:
 
599
                   i0 = k0 = 0;
 
600
                   ti = cos11; impulse_amp *= 2.0;
 
601
                   ts = cos01;
 
602
                   break;
 
603
              case R2R_RODFT00:
 
604
                   i0 = k0 = 1;
 
605
                   ti = sin00; impulse_amp *= 2.0;
 
606
                   ts = cos00;
 
607
                   break;
 
608
              case R2R_RODFT01:
 
609
                   i0 = 1; k0 = 0;
 
610
                   ti = sin01; impulse_amp *= n == 1 ? 1.0 : 2.0;
 
611
                   ts = cos01;
 
612
                   break;
 
613
              case R2R_RODFT10:
 
614
                   i0 = 0; k0 = 1;
 
615
                   ti = sin10; impulse_amp *= 2.0;
 
616
                   ts = cos00;
 
617
                   break;
 
618
              case R2R_RODFT11:
 
619
                   i0 = k0 = 0;
 
620
                   ti = sin11; impulse_amp *= 2.0;
 
621
                   ts = cos01;
 
622
                   break;
 
623
              default:
 
624
                   BENCH_ASSERT(0);
 
625
                   return;
 
626
          }
 
627
 
 
628
          d[i].n0 = n0;
 
629
          d[i].i0 = i0;
 
630
          d[i].k0 = k0;
 
631
          d[i].ti = ti;
 
632
          d[i].ts = ts;
 
633
     }
 
634
 
 
635
 
 
636
     inA = (R *) bench_malloc(N * sizeof(R));
 
637
     inB = (R *) bench_malloc(N * sizeof(R));
 
638
     inC = (R *) bench_malloc(N * sizeof(R));
 
639
     outA = (R *) bench_malloc(N * sizeof(R));
 
640
     outB = (R *) bench_malloc(N * sizeof(R));
 
641
     outC = (R *) bench_malloc(N * sizeof(R));
 
642
     tmp = (R *) bench_malloc(N * sizeof(R));
 
643
 
 
644
     nfo.p = p;
 
645
     nfo.probsz = p->sz;
 
646
     nfo.totalsz = tensor_append(p->vecsz, nfo.probsz);
 
647
     nfo.pckdsz = verify_pack(nfo.totalsz, 1);
 
648
     nfo.pckdvecsz = verify_pack(p->vecsz, tensor_sz(nfo.probsz));
 
649
 
 
650
     e->i = rimpulse(d, impulse_amp, n, vecn, &nfo,
 
651
                     inA, inB, inC, outA, outB, outC, tmp, rounds, tol);
 
652
     e->l = rlinear(N, &nfo, inA, inB, inC, outA, outB, outC, tmp, rounds,tol);
 
653
     e->s = t_shift(n, vecn, &nfo, inA, inB, outA, outB, tmp, 
 
654
                    rounds, tol, d);
 
655
 
 
656
     /* grr, verify-lib.c:preserves_input() only works for complex */
 
657
     if (!p->in_place && !p->destroy_input) {
 
658
          bench_tensor *totalsz_swap, *pckdsz_swap;
 
659
          totalsz_swap = tensor_copy_swapio(nfo.totalsz);
 
660
          pckdsz_swap = tensor_copy_swapio(nfo.pckdsz);
 
661
 
 
662
          for (i = 0; i < rounds; ++i) {
 
663
               rarand(inA, N);
 
664
               dofft(&nfo, inA, outB);
 
665
               cpyr((R *) nfo.p->in, totalsz_swap, inB, pckdsz_swap);
 
666
               racmp(inB, inA, N, "preserves_input", 0.0);
 
667
          }
 
668
 
 
669
          tensor_destroy(totalsz_swap);
 
670
          tensor_destroy(pckdsz_swap);
 
671
     }
 
672
 
 
673
     tensor_destroy(nfo.totalsz);
 
674
     tensor_destroy(nfo.pckdsz);
 
675
     tensor_destroy(nfo.pckdvecsz);
 
676
     bench_free(tmp);
 
677
     bench_free(outC);
 
678
     bench_free(outB);
 
679
     bench_free(outA);
 
680
     bench_free(inC);
 
681
     bench_free(inB);
 
682
     bench_free(inA);
 
683
     bench_free(d);
 
684
}
 
685
 
 
686
 
 
687
typedef struct {
 
688
     dofft_closure k;
 
689
     bench_problem *p;
 
690
     int n0;
 
691
} dofft_r2r_closure;
 
692
 
 
693
static void cpyr1(int n, R *in, int is, R *out, int os, R scale)
 
694
{
 
695
     int i;
 
696
     for (i = 0; i < n; ++i)
 
697
          out[i * os] = in[i * is] * scale;
 
698
}
 
699
 
 
700
static void mke00(C *a, int n, int c)
 
701
{
 
702
     int i;
 
703
     for (i = 1; i + i < n; ++i)
 
704
          a[n - i][c] = a[i][c];
 
705
}
 
706
 
 
707
static void mkre00(C *a, int n)
 
708
{
 
709
     mkreal(a, n);
 
710
     mke00(a, n, 0);
 
711
}
 
712
 
 
713
static void mkimag(C *a, int n)
 
714
{
 
715
     int i;
 
716
     for (i = 0; i < n; ++i)
 
717
          c_re(a[i]) = 0.0;
 
718
}
 
719
 
 
720
static void mko00(C *a, int n, int c)
 
721
{
 
722
     int i;
 
723
     a[0][c] = 0.0;
 
724
     for (i = 1; i + i < n; ++i)
 
725
          a[n - i][c] = -a[i][c];
 
726
     if (i + i == n)
 
727
          a[i][c] = 0.0;
 
728
}
 
729
 
 
730
static void mkro00(C *a, int n)
 
731
{
 
732
     mkreal(a, n);
 
733
     mko00(a, n, 0);
 
734
}
 
735
 
 
736
static void mkio00(C *a, int n)
 
737
{
 
738
     mkimag(a, n);
 
739
     mko00(a, n, 1);
 
740
}
 
741
 
 
742
static void mkre01(C *a, int n) /* n should be be multiple of 4 */
 
743
{
 
744
     R a0;
 
745
     a0 = c_re(a[0]);
 
746
     mko00(a, n/2, 0);
 
747
     c_re(a[n/2]) = -(c_re(a[0]) = a0);
 
748
     mkre00(a, n);
 
749
}
 
750
 
 
751
static void mkro01(C *a, int n) /* n should be be multiple of 4 */
 
752
{
 
753
     c_re(a[0]) = c_im(a[0]) = 0.0;
 
754
     mkre00(a, n/2);
 
755
     mkro00(a, n);
 
756
}
 
757
 
 
758
static void mkoddonly(C *a, int n)
 
759
{
 
760
     int i;
 
761
     for (i = 0; i < n; i += 2)
 
762
          c_re(a[i]) = c_im(a[i]) = 0.0;
 
763
}
 
764
 
 
765
static void mkre10(C *a, int n)
 
766
{
 
767
     mkoddonly(a, n);
 
768
     mkre00(a, n);
 
769
}
 
770
 
 
771
static void mkio10(C *a, int n)
 
772
{
 
773
     mkoddonly(a, n);
 
774
     mkio00(a, n);
 
775
}
 
776
 
 
777
static void mkre11(C *a, int n)
 
778
{
 
779
     mkoddonly(a, n);
 
780
     mko00(a, n/2, 0);
 
781
     mkre00(a, n);
 
782
}
 
783
 
 
784
static void mkro11(C *a, int n)
 
785
{
 
786
     mkoddonly(a, n);
 
787
     mkre00(a, n/2);
 
788
     mkro00(a, n);
 
789
}
 
790
 
 
791
static void mkio11(C *a, int n)
 
792
{
 
793
     mkoddonly(a, n);
 
794
     mke00(a, n/2, 1);
 
795
     mkio00(a, n);
 
796
}
 
797
 
 
798
static void r2r_apply(dofft_closure *k_, bench_complex *in, bench_complex *out)
 
799
{
 
800
     dofft_r2r_closure *k = (dofft_r2r_closure *)k_;
 
801
     bench_problem *p = k->p;
 
802
     bench_real *ri, *ro;
 
803
     int n, is, os;
 
804
 
 
805
     n = p->sz->dims[0].n;
 
806
     is = p->sz->dims[0].is;
 
807
     os = p->sz->dims[0].os;
 
808
 
 
809
     ri = (bench_real *) p->in;
 
810
     ro = (bench_real *) p->out;
 
811
 
 
812
     switch (p->k[0]) {
 
813
         case R2R_R2HC:
 
814
              cpyr1(n, &c_re(in[0]), 2, ri, is, 1.0);
 
815
              break;
 
816
         case R2R_HC2R:
 
817
              cpyr1(n/2 + 1, &c_re(in[0]), 2, ri, is, 1.0);
 
818
              cpyr1((n+1)/2 - 1, &c_im(in[n-1]), -2, ri + is*(n-1), -is, 1.0);
 
819
              break;
 
820
         case R2R_REDFT00:
 
821
              cpyr1(n, &c_re(in[0]), 2, ri, is, 1.0);
 
822
              break;
 
823
         case R2R_RODFT00:
 
824
              cpyr1(n, &c_re(in[1]), 2, ri, is, 1.0);
 
825
              break;
 
826
         case R2R_REDFT01:
 
827
              cpyr1(n, &c_re(in[0]), 2, ri, is, 1.0);
 
828
              break;
 
829
         case R2R_REDFT10:
 
830
              cpyr1(n, &c_re(in[1]), 4, ri, is, 1.0);
 
831
              break;
 
832
         case R2R_RODFT01:
 
833
              cpyr1(n, &c_re(in[1]), 2, ri, is, 1.0);
 
834
              break;
 
835
         case R2R_RODFT10:
 
836
              cpyr1(n, &c_im(in[1]), 4, ri, is, 1.0);
 
837
              break;
 
838
         case R2R_REDFT11:
 
839
              cpyr1(n, &c_re(in[1]), 4, ri, is, 1.0);
 
840
              break;
 
841
         case R2R_RODFT11:
 
842
              cpyr1(n, &c_re(in[1]), 4, ri, is, 1.0);
 
843
              break;
 
844
         default:
 
845
              BENCH_ASSERT(0); /* not yet implemented */
 
846
     }
 
847
 
 
848
     doit(1, p);
 
849
 
 
850
     switch (p->k[0]) {
 
851
         case R2R_R2HC:
 
852
              if (k->k.recopy_input)
 
853
                   cpyr1(n, ri, is, &c_re(in[0]), 2, 1.0);
 
854
              cpyr1(n/2 + 1, ro, os, &c_re(out[0]), 2, 1.0);
 
855
              cpyr1((n+1)/2 - 1, ro + os*(n-1), -os, &c_im(out[1]), 2, 1.0);
 
856
              c_im(out[0]) = 0.0;
 
857
              if (n % 2 == 0)
 
858
                   c_im(out[n/2]) = 0.0;
 
859
              mkhermitian1(out, n);
 
860
              break;
 
861
         case R2R_HC2R:
 
862
              if (k->k.recopy_input) {
 
863
                   cpyr1(n/2 + 1, ri, is, &c_re(in[0]), 2, 1.0);
 
864
                   cpyr1((n+1)/2 - 1, ri + is*(n-1), -is, &c_im(in[1]), 2,1.0);
 
865
              }
 
866
              cpyr1(n, ro, os, &c_re(out[0]), 2, 1.0);
 
867
              mkreal(out, n);
 
868
              break;
 
869
         case R2R_REDFT00:
 
870
              if (k->k.recopy_input)
 
871
                   cpyr1(n, ri, is, &c_re(in[0]), 2, 1.0);
 
872
              cpyr1(n, ro, os, &c_re(out[0]), 2, 1.0);
 
873
              mkre00(out, k->n0);
 
874
              break;
 
875
         case R2R_RODFT00:
 
876
              if (k->k.recopy_input)
 
877
                   cpyr1(n, ri, is, &c_im(in[1]), 2, -1.0);
 
878
              cpyr1(n, ro, os, &c_im(out[1]), 2, -1.0);
 
879
              mkio00(out, k->n0);
 
880
              break;
 
881
         case R2R_REDFT01:
 
882
              if (k->k.recopy_input)
 
883
                   cpyr1(n, ri, is, &c_re(in[0]), 2, 1.0);
 
884
              cpyr1(n, ro, os, &c_re(out[1]), 4, 2.0);
 
885
              mkre10(out, k->n0);
 
886
              break;
 
887
         case R2R_REDFT10:
 
888
              if (k->k.recopy_input)
 
889
                   cpyr1(n, ri, is, &c_re(in[1]), 4, 2.0);
 
890
              cpyr1(n, ro, os, &c_re(out[0]), 2, 1.0);
 
891
              mkre01(out, k->n0);
 
892
              break;
 
893
         case R2R_RODFT01:
 
894
              if (k->k.recopy_input)
 
895
                   cpyr1(n, ri, is, &c_re(in[1]), 2, 1.0);
 
896
              cpyr1(n, ro, os, &c_im(out[1]), 4, -2.0);
 
897
              mkio10(out, k->n0);
 
898
              break;
 
899
         case R2R_RODFT10:
 
900
              if (k->k.recopy_input)
 
901
                   cpyr1(n, ri, is, &c_im(in[1]), 4, -2.0);
 
902
              cpyr1(n, ro, os, &c_re(out[1]), 2, 1.0);
 
903
              mkro01(out, k->n0);
 
904
              break;
 
905
         case R2R_REDFT11:
 
906
              if (k->k.recopy_input)
 
907
                   cpyr1(n, ri, is, &c_re(in[1]), 4, 2.0);
 
908
              cpyr1(n, ro, os, &c_re(out[1]), 4, 2.0);
 
909
              mkre11(out, k->n0);
 
910
              break;
 
911
         case R2R_RODFT11:
 
912
              if (k->k.recopy_input)
 
913
                   cpyr1(n, ri, is, &c_im(in[1]), 4, -2.0);
 
914
              cpyr1(n, ro, os, &c_im(out[1]), 4, -2.0);
 
915
              mkio11(out, k->n0);
 
916
              break;
 
917
         default:
 
918
              BENCH_ASSERT(0); /* not yet implemented */
 
919
     }
 
920
}
 
921
 
 
922
void accuracy_r2r(bench_problem *p, int rounds, int impulse_rounds,
 
923
                  double t[6])
 
924
{
 
925
     dofft_r2r_closure k;
 
926
     int n, n0 = 1;
 
927
     C *a, *b;
 
928
     aconstrain constrain = 0;
 
929
 
 
930
     BENCH_ASSERT(p->kind == PROBLEM_R2R);
 
931
     BENCH_ASSERT(p->sz->rnk == 1);
 
932
     BENCH_ASSERT(p->vecsz->rnk == 0);
 
933
 
 
934
     k.k.apply = r2r_apply;
 
935
     k.k.recopy_input = 0;
 
936
     k.p = p;
 
937
     n = tensor_sz(p->sz);
 
938
     
 
939
     switch (p->k[0]) {
 
940
         case R2R_R2HC: constrain = mkreal; n0 = n; break;
 
941
         case R2R_HC2R: constrain = mkhermitian1; n0 = n; break;
 
942
         case R2R_REDFT00: constrain = mkre00; n0 = 2*(n-1); break;
 
943
         case R2R_RODFT00: constrain = mkro00; n0 = 2*(n+1); break;
 
944
         case R2R_REDFT01: constrain = mkre01; n0 = 4*n; break;
 
945
         case R2R_REDFT10: constrain = mkre10; n0 = 4*n; break;
 
946
         case R2R_RODFT01: constrain = mkro01; n0 = 4*n; break;
 
947
         case R2R_RODFT10: constrain = mkio10; n0 = 4*n; break;
 
948
         case R2R_REDFT11: constrain = mkre11; n0 = 8*n; break;
 
949
         case R2R_RODFT11: constrain = mkro11; n0 = 8*n; break;
 
950
         default: BENCH_ASSERT(0); /* not yet implemented */
 
951
     }
 
952
     k.n0 = n0;
 
953
 
 
954
     a = (C *) bench_malloc(n0 * sizeof(C));
 
955
     b = (C *) bench_malloc(n0 * sizeof(C));
 
956
     accuracy_test(&k.k, constrain, -1, n0, a, b, rounds, impulse_rounds, t);
 
957
     bench_free(b);
 
958
     bench_free(a);
 
959
}