← Back to branch summary

~ubuntu-branches/ubuntu/raring/fftw3/raring-proposed

« back to all changes in this revision

Viewing changes to rdft/generic.c

  • Committer: Bazaar Package Importer
  • Author(s): Paul Brossier
  • Date: 2006-05-31 13:44:05 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20060531134405-ol9hrbg6bh81sg0c
Tags: 3.1.1-1
http://bugs.debian.org/350327
http://bugs.debian.org/338487
http://bugs.debian.org/338501
* New upstream release (closes: #350327, #338487, #338501)
* Add --enable-portable-binary to use -mtune instead of -march
* Use --with-gcc-arch=G5 / pentium4 on powerpc / i386
* Updated Standards-Version

Show diffs side-by-side

added added

removed removed

Lines of Context:
rdft/generic.c
1
1
/*
2
 
 * Copyright (c) 2003 Matteo Frigo
3
 
 * Copyright (c) 2003 Massachusetts Institute of Technology
 
2
 * Copyright (c) 2003, 2006 Matteo Frigo
 
3
 * Copyright (c) 2003, 2006 Massachusetts Institute of Technology
4
4
 *
5
5
 * This program is free software; you can redistribute it and/or modify
6
6
 * it under the terms of the GNU General Public License as published by
27
27
 
28
28
typedef struct {
29
29
     plan_rdft super;
30
 
     plan *cld;
31
30
     twid *td;
32
 
     int os;
33
 
     int r, m;
 
31
     INT n, is, os;
34
32
     rdft_kind kind;
35
33
} P;
36
34
 
37
35
/***************************************************************************/
38
36
 
39
 
static void apply_dit(const plan *ego_, R *I, R *O)
40
 
{
41
 
     const P *ego = (const P *) ego_;
42
 
     int n, m, r;
43
 
     int i, j, k;
44
 
     int os, osm;
45
 
     E *buf;
46
 
     const R *W;
47
 
     R *X, *YO, *YI;
48
 
     E rsum, isum;
49
 
     int wp, wincr;
50
 
 
51
 
     {
52
 
          plan_rdft *cld = (plan_rdft *) ego->cld;
53
 
          cld->apply((plan *) cld, I, O);
54
 
     }
55
 
 
56
 
     r = ego->r;
57
 
 
58
 
     STACK_MALLOC(E *, buf, r * 2 * sizeof(E));
59
 
     
60
 
     osm = (m = ego->m) * (os = ego->os);
61
 
     n = m * r;
62
 
     W = ego->td->W;
63
 
 
64
 
     X = O;
65
 
     YO = O + r * osm;
66
 
     YI = O + osm;
67
 
 
68
 
     /* compute the transform of the r 0th elements (which are real) */
69
 
     for (i = 0; i + i < r; ++i) {
70
 
          rsum = K(0.0);
71
 
          isum = K(0.0);
72
 
          wincr = m * i;
73
 
          for (j = 0, wp = 0; j < r; ++j) {
74
 
               E tw_r = W[2*wp];
75
 
               E tw_i = W[2*wp+1] ;
76
 
               E re = X[j * osm];
77
 
               rsum += re * tw_r;
78
 
               isum += re * tw_i;
79
 
               wp += wincr;
80
 
               if (wp >= n)
81
 
                    wp -= n;
82
 
          }
83
 
          buf[2*i] = rsum;
84
 
          buf[2*i+1] = isum;
85
 
     }
86
 
 
87
 
     /* store the transform back onto the A array */
88
 
     X[0] = buf[0];
89
 
     for (i = 1; i + i < r; ++i) {
90
 
          X[i * osm] = buf[2*i];
91
 
          YO[-i * osm] = buf[2*i+1];
92
 
     }
93
 
 
94
 
     X += os;
95
 
     YI -= os;
96
 
     YO -= os;
97
 
 
98
 
     /* compute the transform of the middle elements (which are complex) */
99
 
     for (k = 1; k + k < m; ++k, X += os, YI -= os, YO -= os) {
100
 
          for (i = 0; i < r; ++i) {
101
 
               rsum = K(0.0);
102
 
               isum = K(0.0);
103
 
               wincr = k + m * i;
104
 
               for (j = 0, wp = 0; j < r; ++j) {
105
 
                    E tw_r = W[2*wp];
106
 
                    E tw_i = W[2*wp+1] ;
107
 
                    E re = X[j * osm];
108
 
                    E im = YI[j * osm];
109
 
                    rsum += re * tw_r - im * tw_i;
110
 
                    isum += re * tw_i + im * tw_r;
111
 
                    wp += wincr;
112
 
                    if (wp >= n)
113
 
                         wp -= n;
114
 
               }
115
 
               buf[2*i] = rsum;
116
 
               buf[2*i+1] = isum;
117
 
          }
118
 
 
119
 
          /* store the transform back onto the A array */
120
 
          for (i = 0; i + i < r; ++i) {
121
 
               X[i * osm] = buf[2*i];
122
 
               YO[-i * osm] = buf[2*i+1];
123
 
          }
124
 
          for (; i < r; ++i) {
125
 
               X[i * osm] = -buf[2*i+1];
126
 
               YO[-i * osm] = buf[2*i];
127
 
          }
128
 
     }
129
 
 
130
 
     /* no final element, since m is odd */
131
 
 
132
 
     STACK_FREE(buf);
133
 
}
134
 
 
135
 
static void apply_dif(const plan *ego_, R *I, R *O)
136
 
{
137
 
     const P *ego = (const P *) ego_;
138
 
     int n, m, r;
139
 
     int i, j, k;
140
 
     int is, ism;
141
 
     E *buf;
142
 
     const R *W;
143
 
     R *X, *YO, *YI;
144
 
     E rsum, isum;
145
 
     int wp, wincr;
146
 
 
147
 
     r = ego->r;
148
 
 
149
 
     STACK_MALLOC(E *, buf, r * 2 * sizeof(E));
150
 
     
151
 
     ism = (m = ego->m) * (is = ego->os);
152
 
     n = m * r;
153
 
     W = ego->td->W;
154
 
 
155
 
     X = I;
156
 
     YI = I + r * ism;
157
 
     YO = I + ism;
158
 
 
159
 
     /* 
160
 
      * compute the transform of the r 0th elements (which are halfcomplex)
161
 
      * yielding real numbers
162
 
      */
163
 
     /* copy the input into the temporary array */
164
 
     buf[0] = X[0];
165
 
     for (i = 1; i + i < r; ++i) {
166
 
          buf[2*i] = X[i * ism];
167
 
          buf[2*i+1] = YI[-i * ism];
168
 
     }
169
 
 
170
 
     for (i = 0; i < r; ++i) {
171
 
          rsum = K(0.0);
172
 
          wincr = m * i;
173
 
          for (j = 1, wp = wincr; j + j < r; ++j) {
174
 
               E tw_r = W[2*wp];
175
 
               E tw_i = W[2*wp+1];
176
 
               E re = buf[2*j];
177
 
               E im = buf[2*j+1];
178
 
               rsum += re * tw_r + im * tw_i;
179
 
               wp += wincr;
180
 
               if (wp >= n)
181
 
                    wp -= n;
182
 
          }
183
 
          X[i * ism] = K(2.0) * rsum + buf[0];
184
 
     }
185
 
 
186
 
     X += is;
187
 
     YI -= is;
188
 
     YO -= is;
189
 
 
190
 
     /* compute the transform of the middle elements (which are complex) */
191
 
     for (k = 1; k + k < m; ++k, X += is, YI -= is, YO -= is) {
192
 
          /* copy the input into the temporary array */
193
 
          for (i = 0; i + i < r; ++i) {
194
 
               buf[2*i] = X[i * ism];
195
 
               buf[2*i+1] = YI[-i * ism];
196
 
          }
197
 
          for (; i < r; ++i) {
198
 
               buf[2*i+1] = -X[i * ism];
199
 
               buf[2*i] = YI[-i * ism];
200
 
          }
201
 
 
202
 
          for (i = 0; i < r; ++i) {
203
 
               rsum = K(0.0);
204
 
               isum = K(0.0);
205
 
               wincr = m * i;
206
 
               for (j = 0, wp = k * i; j < r; ++j) {
207
 
                    E tw_r = W[2*wp];
208
 
                    E tw_i = W[2*wp+1];
209
 
                    E re = buf[2*j];
210
 
                    E im = buf[2*j+1];
211
 
                    rsum += re * tw_r + im * tw_i;
212
 
                    isum += im * tw_r - re * tw_i;
213
 
                    wp += wincr;
214
 
                    if (wp >= n)
215
 
                         wp -= n;
216
 
               }
217
 
               X[i * ism] = rsum;
218
 
               YO[i * ism] = isum;
219
 
          }
220
 
     }
221
 
 
222
 
     /* no final element, since m is odd */
223
 
 
224
 
     STACK_FREE(buf);
225
 
 
226
 
     {
227
 
          plan_rdft *cld = (plan_rdft *) ego->cld;
228
 
          cld->apply((plan *) cld, I, O);
229
 
     }
230
 
 
231
 
}
 
37
static void cdot_r2hc(INT n, const E *x, const R *w, R *or0, R *oi1)
 
38
{
 
39
     INT i;
 
40
 
 
41
     E rr = x[0], ri = 0;
 
42
     x += 1;
 
43
     for (i = 1; i + i < n; ++i) {
 
44
          rr += x[0] * w[0];
 
45
          ri += x[1] * w[1];
 
46
          x += 2; w += 2;
 
47
     }
 
48
     *or0 = rr;
 
49
     *oi1 = ri;
 
50
}
 
51
 
 
52
static void hartley_r2hc(INT n, const R *xr, INT xs, E *o, R *pr)
 
53
{
 
54
     INT i;
 
55
     E sr;
 
56
     o[0] = sr = xr[0]; o += 1;
 
57
     for (i = 1; i + i < n; ++i) {
 
58
          R a, b;
 
59
          a = xr[i * xs];
 
60
          b =  xr[(n - i) * xs];
 
61
          sr += (o[0] = a + b);
 
62
#if FFT_SIGN == -1
 
63
          o[1] = b - a;
 
64
#else
 
65
          o[1] = a - b;
 
66
#endif
 
67
          o += 2;
 
68
     }
 
69
     *pr = sr;
 
70
}
 
71
                    
 
72
static void apply_r2hc(const plan *ego_, R *I, R *O)
 
73
{
 
74
     const P *ego = (const P *) ego_;
 
75
     INT i;
 
76
     INT n = ego->n, is = ego->is, os = ego->os;
 
77
     const R *W = ego->td->W;
 
78
     E *buf;
 
79
 
 
80
     STACK_MALLOC(E *, buf, n * sizeof(E));
 
81
     hartley_r2hc(n, I, is, buf, O);
 
82
 
 
83
     for (i = 1; i + i < n; ++i) {
 
84
          cdot_r2hc(n, buf, W, O + i * os, O + (n - i) * os);
 
85
          W += n - 1;
 
86
     }
 
87
 
 
88
     STACK_FREE(buf);
 
89
}
 
90
 
 
91
 
 
92
static void cdot_hc2r(INT n, const E *x, const R *w, R *or0, R *or1)
 
93
{
 
94
     INT i;
 
95
 
 
96
     E rr = x[0], ii = 0; 
 
97
     x += 1;
 
98
     for (i = 1; i + i < n; ++i) {
 
99
          rr += x[0] * w[0];
 
100
          ii += x[1] * w[1];
 
101
          x += 2; w += 2;
 
102
     }
 
103
#if FFT_SIGN == -1
 
104
     *or0 = rr - ii;
 
105
     *or1 = rr + ii;
 
106
#else
 
107
     *or0 = rr + ii;
 
108
     *or1 = rr - ii;
 
109
#endif
 
110
}
 
111
 
 
112
static void hartley_hc2r(INT n, const R *x, INT xs, E *o, R *pr)
 
113
{
 
114
     INT i;
 
115
     E sr;
 
116
 
 
117
     o[0] = sr = x[0]; o += 1;
 
118
     for (i = 1; i + i < n; ++i) {
 
119
          sr += (o[0] = x[i * xs] + x[i * xs]);
 
120
          o[1] = x[(n - i) * xs] + x[(n - i) * xs];
 
121
          o += 2;
 
122
     }
 
123
     *pr = sr;
 
124
}
 
125
 
 
126
static void apply_hc2r(const plan *ego_, R *I, R *O)                
 
127
{
 
128
     const P *ego = (const P *) ego_;
 
129
     INT i;
 
130
     INT n = ego->n, is = ego->is, os = ego->os;
 
131
     const R *W = ego->td->W;
 
132
     E *buf;
 
133
 
 
134
     STACK_MALLOC(E *, buf, n * sizeof(E));
 
135
     hartley_hc2r(n, I, is, buf, O);
 
136
 
 
137
     for (i = 1; i + i < n; ++i) {
 
138
          cdot_hc2r(n, buf, W, O + i * os, O + (n - i) * os);
 
139
          W += n - 1;
 
140
     }
 
141
 
 
142
     STACK_FREE(buf);
 
143
}
 
144
 
232
145
 
233
146
/***************************************************************************/
234
147
 
235
 
static void awake(plan *ego_, int flg)
 
148
static void awake(plan *ego_, enum wakefulness wakefulness)
236
149
{
237
150
     P *ego = (P *) ego_;
238
 
     static const tw_instr generic_tw[] = {
239
 
          { TW_GENERIC, 0, 0 },
 
151
     static const tw_instr half_tw[] = {
 
152
          { TW_HALF, 1, 0 },
240
153
          { TW_NEXT, 1, 0 }
241
154
     };
242
155
 
243
 
     AWAKE(ego->cld, flg);
244
 
     /* FIXME: can we get away with fewer twiddles? */
245
 
     X(twiddle_awake)(flg, &ego->td, generic_tw,
246
 
                      ego->r * ego->m, ego->r, ego->m);
247
 
}
248
 
 
249
 
static void destroy(plan *ego_)
250
 
{
251
 
     P *ego = (P *) ego_;
252
 
     X(plan_destroy_internal)(ego->cld);
 
156
     X(twiddle_awake)(wakefulness, &ego->td, half_tw, ego->n, ego->n,
 
157
                      (ego->n - 1) / 2);
253
158
}
254
159
 
255
160
static void print(const plan *ego_, printer *p)
256
161
{
257
162
     const P *ego = (const P *) ego_;
258
163
 
259
 
     p->print(p, "(rdft-generic-%s-%d%(%p%))", 
260
 
              ego->kind == R2HC ? "r2hc-dit" : "hc2r-dif",
261
 
              ego->r, ego->cld);
 
164
     p->print(p, "(rdft-generic-%s-%D)", 
 
165
              ego->kind == R2HC ? "r2hc" : "hc2r", 
 
166
              ego->n);
262
167
}
263
168
 
264
 
static int applicable0(const solver *ego_, const problem *p_)
 
169
static int applicable0(const S *ego, const problem *p_)
265
170
{
266
 
     if (RDFTP(p_)) {
267
 
          const S *ego = (const S *) ego_;
268
 
          const problem_rdft *p = (const problem_rdft *) p_;
269
 
          return (1
270
 
                  && p->sz->rnk == 1
271
 
                  && p->vecsz->rnk == 0
272
 
                  && p->sz->dims[0].n > 1
273
 
                  && p->sz->dims[0].n % 2 /* ensure r and n/r odd */
274
 
                  && p->kind[0] == ego->kind
275
 
               );
276
 
     }
277
 
 
278
 
     return 0;
 
171
     const problem_rdft *p = (const problem_rdft *) p_;
 
172
     return (1
 
173
             && p->sz->rnk == 1
 
174
             && p->vecsz->rnk == 0
 
175
             && (p->sz->dims[0].n % 2) == 1 
 
176
             && X(is_prime)(p->sz->dims[0].n)
 
177
             && p->kind[0] == ego->kind
 
178
          );
279
179
}
280
180
 
281
 
static int applicable(const solver *ego_, const problem *p_, 
 
181
static int applicable(const S *ego, const problem *p_, 
282
182
                      const planner *plnr)
283
183
{
284
 
     if (NO_UGLYP(plnr)) return 0; /* always ugly */
285
 
     if (!applicable0(ego_, p_)) return 0;
 
184
     if (NO_SLOWP(plnr)) return 0;
 
185
     if (!applicable0(ego, p_)) return 0;
286
186
 
287
187
     if (NO_LARGE_GENERICP(plnr)) {
288
188
          const problem_rdft *p = (const problem_rdft *) p_;
289
 
          if (X(first_divisor)(p->sz->dims[0].n) >= GENERIC_MIN_BAD) return 0; 
 
189
          if (p->sz->dims[0].n >= GENERIC_MIN_BAD) return 0; 
290
190
     }
291
191
     return 1;
292
192
}
293
193
 
294
 
static plan *mkplan(const solver *ego, const problem *p_, planner *plnr)
 
194
static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
295
195
{
296
 
     const problem_rdft *p = (const problem_rdft *) p_;
297
 
     P *pln = 0;
298
 
     int n, r, m;
299
 
     int is, os;
300
 
     plan *cld = (plan *) 0;
301
 
     problem *cldp;
 
196
     const S *ego = (const S *)ego_;
 
197
     const problem_rdft *p;
 
198
     P *pln;
 
199
     INT n;
302
200
 
303
201
     static const plan_adt padt = {
304
 
          X(rdft_solve), awake, print, destroy
 
202
          X(rdft_solve), awake, print, X(plan_null_destroy)
305
203
     };
306
204
 
307
205
     if (!applicable(ego, p_, plnr))
308
 
          goto nada;
309
 
 
310
 
     n = p->sz->dims[0].n;
311
 
     is = p->sz->dims[0].is;
312
 
     os = p->sz->dims[0].os;
313
 
 
314
 
     r = X(first_divisor)(n);
315
 
     m = n / r;
316
 
 
317
 
     if (R2HC_KINDP(p->kind[0])) {
318
 
          cldp = X(mkproblem_rdft_d)(X(mktensor_1d)(m, r * is, os),
319
 
                                     X(mktensor_1d)(r, is, m * os),
320
 
                                     p->I, p->O, p->kind);
321
 
     }
322
 
     else {
323
 
          cldp = X(mkproblem_rdft_d)(X(mktensor_1d)(m, is, r * os),
324
 
                                     X(mktensor_1d)(r, m * is, os),
325
 
                                     p->I, p->O, p->kind);
326
 
     }
327
 
     if (!(cld = X(mkplan_d)(plnr, cldp))) goto nada;
328
 
 
329
 
     pln = MKPLAN_RDFT(P, &padt, R2HC_KINDP(p->kind[0]) ? apply_dit:apply_dif);
330
 
 
331
 
     pln->os = R2HC_KINDP(p->kind[0]) ? os : is;
332
 
     pln->r = r;
333
 
     pln->m = m;
334
 
     pln->cld = cld;
 
206
          return (plan *)0;
 
207
 
 
208
     p = (const problem_rdft *) p_;
 
209
     pln = MKPLAN_RDFT(P, &padt, 
 
210
                       R2HC_KINDP(p->kind[0]) ? apply_r2hc : apply_hc2r);
 
211
 
 
212
     pln->n = n = p->sz->dims[0].n;
 
213
     pln->is = p->sz->dims[0].is;
 
214
     pln->os = p->sz->dims[0].os;
335
215
     pln->td = 0;
336
 
     pln->kind = p->kind[0];
 
216
     pln->kind = ego->kind;
337
217
 
338
 
     X(ops_zero)(&pln->super.super.ops);
339
 
     pln->super.super.ops.add = 4 * r * r;
340
 
     pln->super.super.ops.mul = 4 * r * r;
341
 
     /* loads + stores, minus loads + stores for all DIT codelets */
342
 
     pln->super.super.ops.other = 4 * r + 4 * r * r - (6*r - 2);
343
 
     X(ops_madd)((m - 1)/2, &pln->super.super.ops, &cld->ops,
344
 
                 &pln->super.super.ops);
345
 
     pln->super.super.ops.add += 2 * r * r;
346
 
     pln->super.super.ops.mul += 2 * r * r;
347
 
     pln->super.super.ops.other += 3 * r + 3 * r * r - 2*r;
 
218
     pln->super.super.ops.add = (n-1) * 2.5;
 
219
     pln->super.super.ops.mul = 0;
 
220
     pln->super.super.ops.fma = 0.5 * (n-1) * (n-1) ;
 
221
#if 0 /* these are nice pipelined sequential loads and should cost nothing */
 
222
     pln->super.super.ops.other = (n-1)*(2 + 1 + (n-1));  /* approximate */
 
223
#endif
348
224
 
349
225
     return &(pln->super.super);
350
 
 
351
 
 nada:
352
 
     X(plan_destroy_internal)(cld);
353
 
     X(ifree0)(pln);
354
 
     return (plan *) 0;
355
226
}
356
227
 
357
 
/* constructors */
358
 
 
359
228
static solver *mksolver(rdft_kind kind)
360
229
{
361
 
     static const solver_adt sadt = { mkplan };
 
230
     static const solver_adt sadt = { PROBLEM_RDFT, mkplan };
362
231
     S *slv = MKSOLVER(S, &sadt);
363
232
     slv->kind = kind;
364
233
     return &(slv->super);