← Back to branch summary

~diresu/blender/blender-command-port

« back to all changes in this revision

Viewing changes to extern/fftw/reodft/reodft010e-r2hc.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/reodft/reodft010e-r2hc.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
/* $Id: reodft010e-r2hc.c,v 1.37 2006-01-27 02:10:50 athena Exp $ */
 
22
 
 
23
/* Do an R{E,O}DFT{01,10} problem via an R2HC problem, with some
 
24
   pre/post-processing ala FFTPACK. */
 
25
 
 
26
#include "reodft.h"
 
27
 
 
28
typedef struct {
 
29
     solver super;
 
30
} S;
 
31
 
 
32
typedef struct {
 
33
     plan_rdft super;
 
34
     plan *cld;
 
35
     twid *td;
 
36
     INT is, os;
 
37
     INT n;
 
38
     INT vl;
 
39
     INT ivs, ovs;
 
40
     rdft_kind kind;
 
41
} P;
 
42
 
 
43
/* A real-even-01 DFT operates logically on a size-4N array:
 
44
                   I 0 -r(I*) -I 0 r(I*),
 
45
   where r denotes reversal and * denotes deletion of the 0th element.
 
46
   To compute the transform of this, we imagine performing a radix-4
 
47
   (real-input) DIF step, which turns the size-4N DFT into 4 size-N
 
48
   (contiguous) DFTs, two of which are zero and two of which are
 
49
   conjugates.  The non-redundant size-N DFT has halfcomplex input, so
 
50
   we can do it with a size-N hc2r transform.  (In order to share
 
51
   plans with the re10 (inverse) transform, however, we use the DHT
 
52
   trick to re-express the hc2r problem as r2hc.  This has little cost
 
53
   since we are already pre- and post-processing the data in {i,n-i}
 
54
   order.)  Finally, we have to write out the data in the correct
 
55
   order...the two size-N redundant (conjugate) hc2r DFTs correspond
 
56
   to the even and odd outputs in O (i.e. the usual interleaved output
 
57
   of DIF transforms); since this data has even symmetry, we only
 
58
   write the first half of it.
 
59
 
 
60
   The real-even-10 DFT is just the reverse of these steps, i.e. a
 
61
   radix-4 DIT transform.  There, however, we just use the r2hc
 
62
   transform naturally without resorting to the DHT trick.
 
63
 
 
64
   A real-odd-01 DFT is very similar, except that the input is
 
65
   0 I (rI)* 0 -I -(rI)*.  This format, however, can be transformed
 
66
   into precisely the real-even-01 format above by sending I -> rI
 
67
   and shifting the array by N.  The former swap is just another
 
68
   transformation on the input during preprocessing; the latter
 
69
   multiplies the even/odd outputs by i/-i, which combines with
 
70
   the factor of -i (to take the imaginary part) to simply flip
 
71
   the sign of the odd outputs.  Vice-versa for real-odd-10.
 
72
 
 
73
   The FFTPACK source code was very helpful in working this out.
 
74
   (They do unnecessary passes over the array, though.)  The same
 
75
   algorithm is also described in:
 
76
 
 
77
      John Makhoul, "A fast cosine transform in one and two dimensions,"
 
78
      IEEE Trans. on Acoust. Speech and Sig. Proc., ASSP-28 (1), 27--34 (1980).
 
79
 
 
80
   Note that Numerical Recipes suggests a different algorithm that
 
81
   requires more operations and uses trig. functions for both the pre-
 
82
   and post-processing passes.
 
83
*/
 
84
 
 
85
static void apply_re01(const plan *ego_, R *I, R *O)
 
86
{
 
87
     const P *ego = (const P *) ego_;
 
88
     INT is = ego->is, os = ego->os;
 
89
     INT i, n = ego->n;
 
90
     INT iv, vl = ego->vl;
 
91
     INT ivs = ego->ivs, ovs = ego->ovs;
 
92
     R *W = ego->td->W;
 
93
     R *buf;
 
94
 
 
95
     buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
96
 
 
97
     for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
 
98
          buf[0] = I[0];
 
99
          for (i = 1; i < n - i; ++i) {
 
100
               E a, b, apb, amb, wa, wb;
 
101
               a = I[is * i];
 
102
               b = I[is * (n - i)];
 
103
               apb = a + b;
 
104
               amb = a - b;
 
105
               wa = W[2*i];
 
106
               wb = W[2*i + 1];
 
107
               buf[i] = wa * amb + wb * apb; 
 
108
               buf[n - i] = wa * apb - wb * amb; 
 
109
          }
 
110
          if (i == n - i) {
 
111
               buf[i] = K(2.0) * I[is * i] * W[2*i];
 
112
          }
 
113
          
 
114
          {
 
115
               plan_rdft *cld = (plan_rdft *) ego->cld;
 
116
               cld->apply((plan *) cld, buf, buf);
 
117
          }
 
118
          
 
119
          O[0] = buf[0];
 
120
          for (i = 1; i < n - i; ++i) {
 
121
               E a, b;
 
122
               INT k;
 
123
               a = buf[i];
 
124
               b = buf[n - i];
 
125
               k = i + i;
 
126
               O[os * (k - 1)] = a - b;
 
127
               O[os * k] = a + b;
 
128
          }
 
129
          if (i == n - i) {
 
130
               O[os * (n - 1)] = buf[i];
 
131
          }
 
132
     }
 
133
 
 
134
     X(ifree)(buf);
 
135
}
 
136
 
 
137
/* ro01 is same as re01, but with i <-> n - 1 - i in the input and
 
138
   the sign of the odd output elements flipped. */
 
139
static void apply_ro01(const plan *ego_, R *I, R *O)
 
140
{
 
141
     const P *ego = (const P *) ego_;
 
142
     INT is = ego->is, os = ego->os;
 
143
     INT i, n = ego->n;
 
144
     INT iv, vl = ego->vl;
 
145
     INT ivs = ego->ivs, ovs = ego->ovs;
 
146
     R *W = ego->td->W;
 
147
     R *buf;
 
148
 
 
149
     buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
150
 
 
151
     for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
 
152
          buf[0] = I[is * (n - 1)];
 
153
          for (i = 1; i < n - i; ++i) {
 
154
               E a, b, apb, amb, wa, wb;
 
155
               a = I[is * (n - 1 - i)];
 
156
               b = I[is * (i - 1)];
 
157
               apb = a + b;
 
158
               amb = a - b;
 
159
               wa = W[2*i];
 
160
               wb = W[2*i+1];
 
161
               buf[i] = wa * amb + wb * apb; 
 
162
               buf[n - i] = wa * apb - wb * amb; 
 
163
          }
 
164
          if (i == n - i) {
 
165
               buf[i] = K(2.0) * I[is * (i - 1)] * W[2*i];
 
166
          }
 
167
          
 
168
          {
 
169
               plan_rdft *cld = (plan_rdft *) ego->cld;
 
170
               cld->apply((plan *) cld, buf, buf);
 
171
          }
 
172
          
 
173
          O[0] = buf[0];
 
174
          for (i = 1; i < n - i; ++i) {
 
175
               E a, b;
 
176
               INT k;
 
177
               a = buf[i];
 
178
               b = buf[n - i];
 
179
               k = i + i;
 
180
               O[os * (k - 1)] = b - a;
 
181
               O[os * k] = a + b;
 
182
          }
 
183
          if (i == n - i) {
 
184
               O[os * (n - 1)] = -buf[i];
 
185
          }
 
186
     }
 
187
 
 
188
     X(ifree)(buf);
 
189
}
 
190
 
 
191
static void apply_re10(const plan *ego_, R *I, R *O)
 
192
{
 
193
     const P *ego = (const P *) ego_;
 
194
     INT is = ego->is, os = ego->os;
 
195
     INT i, n = ego->n;
 
196
     INT iv, vl = ego->vl;
 
197
     INT ivs = ego->ivs, ovs = ego->ovs;
 
198
     R *W = ego->td->W;
 
199
     R *buf;
 
200
 
 
201
     buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
202
 
 
203
     for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
 
204
          buf[0] = I[0];
 
205
          for (i = 1; i < n - i; ++i) {
 
206
               E u, v;
 
207
               INT k = i + i;
 
208
               u = I[is * (k - 1)];
 
209
               v = I[is * k];
 
210
               buf[n - i] = u;
 
211
               buf[i] = v;
 
212
          }
 
213
          if (i == n - i) {
 
214
               buf[i] = I[is * (n - 1)];
 
215
          }
 
216
          
 
217
          {
 
218
               plan_rdft *cld = (plan_rdft *) ego->cld;
 
219
               cld->apply((plan *) cld, buf, buf);
 
220
          }
 
221
          
 
222
          O[0] = K(2.0) * buf[0];
 
223
          for (i = 1; i < n - i; ++i) {
 
224
               E a, b, wa, wb;
 
225
               a = K(2.0) * buf[i];
 
226
               b = K(2.0) * buf[n - i];
 
227
               wa = W[2*i];
 
228
               wb = W[2*i + 1];
 
229
               O[os * i] = wa * a + wb * b;
 
230
               O[os * (n - i)] = wb * a - wa * b;
 
231
          }
 
232
          if (i == n - i) {
 
233
               O[os * i] = K(2.0) * buf[i] * W[2*i];
 
234
          }
 
235
     }
 
236
 
 
237
     X(ifree)(buf);
 
238
}
 
239
 
 
240
/* ro10 is same as re10, but with i <-> n - 1 - i in the output and
 
241
   the sign of the odd input elements flipped. */
 
242
static void apply_ro10(const plan *ego_, R *I, R *O)
 
243
{
 
244
     const P *ego = (const P *) ego_;
 
245
     INT is = ego->is, os = ego->os;
 
246
     INT i, n = ego->n;
 
247
     INT iv, vl = ego->vl;
 
248
     INT ivs = ego->ivs, ovs = ego->ovs;
 
249
     R *W = ego->td->W;
 
250
     R *buf;
 
251
 
 
252
     buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
253
 
 
254
     for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
 
255
          buf[0] = I[0];
 
256
          for (i = 1; i < n - i; ++i) {
 
257
               E u, v;
 
258
               INT k = i + i;
 
259
               u = -I[is * (k - 1)];
 
260
               v = I[is * k];
 
261
               buf[n - i] = u;
 
262
               buf[i] = v;
 
263
          }
 
264
          if (i == n - i) {
 
265
               buf[i] = -I[is * (n - 1)];
 
266
          }
 
267
          
 
268
          {
 
269
               plan_rdft *cld = (plan_rdft *) ego->cld;
 
270
               cld->apply((plan *) cld, buf, buf);
 
271
          }
 
272
          
 
273
          O[os * (n - 1)] = K(2.0) * buf[0];
 
274
          for (i = 1; i < n - i; ++i) {
 
275
               E a, b, wa, wb;
 
276
               a = K(2.0) * buf[i];
 
277
               b = K(2.0) * buf[n - i];
 
278
               wa = W[2*i];
 
279
               wb = W[2*i + 1];
 
280
               O[os * (n - 1 - i)] = wa * a + wb * b;
 
281
               O[os * (i - 1)] = wb * a - wa * b;
 
282
          }
 
283
          if (i == n - i) {
 
284
               O[os * (i - 1)] = K(2.0) * buf[i] * W[2*i];
 
285
          }
 
286
     }
 
287
 
 
288
     X(ifree)(buf);
 
289
}
 
290
 
 
291
static void awake(plan *ego_, enum wakefulness wakefulness)
 
292
{
 
293
     P *ego = (P *) ego_;
 
294
     static const tw_instr reodft010e_tw[] = {
 
295
          { TW_COS, 0, 1 },
 
296
          { TW_SIN, 0, 1 },
 
297
          { TW_NEXT, 1, 0 }
 
298
     };
 
299
 
 
300
     X(plan_awake)(ego->cld, wakefulness);
 
301
 
 
302
     X(twiddle_awake)(wakefulness, &ego->td, reodft010e_tw, 
 
303
                      4*ego->n, 1, ego->n/2+1);
 
304
}
 
305
 
 
306
static void destroy(plan *ego_)
 
307
{
 
308
     P *ego = (P *) ego_;
 
309
     X(plan_destroy_internal)(ego->cld);
 
310
}
 
311
 
 
312
static void print(const plan *ego_, printer *p)
 
313
{
 
314
     const P *ego = (const P *) ego_;
 
315
     p->print(p, "(%se-r2hc-%D%v%(%p%))",
 
316
              X(rdft_kind_str)(ego->kind), ego->n, ego->vl, ego->cld);
 
317
}
 
318
 
 
319
static int applicable0(const solver *ego_, const problem *p_)
 
320
{
 
321
     const problem_rdft *p = (const problem_rdft *) p_;
 
322
     UNUSED(ego_);
 
323
 
 
324
     return (1
 
325
             && p->sz->rnk == 1
 
326
             && p->vecsz->rnk <= 1
 
327
             && (p->kind[0] == REDFT01 || p->kind[0] == REDFT10
 
328
                 || p->kind[0] == RODFT01 || p->kind[0] == RODFT10)
 
329
          );
 
330
}
 
331
 
 
332
static int applicable(const solver *ego, const problem *p, const planner *plnr)
 
333
{
 
334
     return (!NO_SLOWP(plnr) && applicable0(ego, p));
 
335
}
 
336
 
 
337
static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
 
338
{
 
339
     P *pln;
 
340
     const problem_rdft *p;
 
341
     plan *cld;
 
342
     R *buf;
 
343
     INT n;
 
344
     opcnt ops;
 
345
 
 
346
     static const plan_adt padt = {
 
347
          X(rdft_solve), awake, print, destroy
 
348
     };
 
349
 
 
350
     if (!applicable(ego_, p_, plnr))
 
351
          return (plan *)0;
 
352
 
 
353
     p = (const problem_rdft *) p_;
 
354
 
 
355
     n = p->sz->dims[0].n;
 
356
     buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
357
 
 
358
     cld = X(mkplan_d)(plnr, X(mkproblem_rdft_1_d)(X(mktensor_1d)(n, 1, 1),
 
359
                                                   X(mktensor_0d)(),
 
360
                                                   buf, buf, R2HC));
 
361
     X(ifree)(buf);
 
362
     if (!cld)
 
363
          return (plan *)0;
 
364
 
 
365
     switch (p->kind[0]) {
 
366
         case REDFT01: pln = MKPLAN_RDFT(P, &padt, apply_re01); break;
 
367
         case REDFT10: pln = MKPLAN_RDFT(P, &padt, apply_re10); break;
 
368
         case RODFT01: pln = MKPLAN_RDFT(P, &padt, apply_ro01); break;
 
369
         case RODFT10: pln = MKPLAN_RDFT(P, &padt, apply_ro10); break;
 
370
         default: A(0); return (plan*)0;
 
371
     }
 
372
 
 
373
     pln->n = n;
 
374
     pln->is = p->sz->dims[0].is;
 
375
     pln->os = p->sz->dims[0].os;
 
376
     pln->cld = cld;
 
377
     pln->td = 0;
 
378
     pln->kind = p->kind[0];
 
379
     
 
380
     X(tensor_tornk1)(p->vecsz, &pln->vl, &pln->ivs, &pln->ovs);
 
381
     
 
382
     X(ops_zero)(&ops);
 
383
     ops.other = 4 + (n-1)/2 * 10 + (1 - n % 2) * 5;
 
384
     if (p->kind[0] == REDFT01 || p->kind[0] == RODFT01) {
 
385
          ops.add = (n-1)/2 * 6;
 
386
          ops.mul = (n-1)/2 * 4 + (1 - n % 2) * 2;
 
387
     }
 
388
     else { /* 10 transforms */
 
389
          ops.add = (n-1)/2 * 2;
 
390
          ops.mul = 1 + (n-1)/2 * 6 + (1 - n % 2) * 2;
 
391
     }
 
392
     
 
393
     X(ops_zero)(&pln->super.super.ops);
 
394
     X(ops_madd2)(pln->vl, &ops, &pln->super.super.ops);
 
395
     X(ops_madd2)(pln->vl, &cld->ops, &pln->super.super.ops);
 
396
 
 
397
     return &(pln->super.super);
 
398
}
 
399
 
 
400
/* constructor */
 
401
static solver *mksolver(void)
 
402
{
 
403
     static const solver_adt sadt = { PROBLEM_RDFT, mkplan };
 
404
     S *slv = MKSOLVER(S, &sadt);
 
405
     return &(slv->super);
 
406
}
 
407
 
 
408
void X(reodft010e_r2hc_register)(planner *p)
 
409
{
 
410
     REGISTER_SOLVER(p, mksolver());
 
411
}