← Back to branch summary

~ubuntu-branches/ubuntu/maverick/blender/maverick

« back to all changes in this revision

Viewing changes to extern/fftw/reodft/redft00e-r2hc.c

  • Committer: Bazaar Package Importer
  • Author(s): Khashayar Naderehvandi, Khashayar Naderehvandi, Alessio Treglia
  • Date: 2009-01-22 16:53:59 UTC
  • mfrom: (14.1.1 experimental)
  • Revision ID: james.westby@ubuntu.com-20090122165359-v0996tn7fbit64ni
Tags: 2.48a+dfsg-1ubuntu1
https://launchpad.net/bugs/320045
[ Khashayar Naderehvandi ]
* Merge from debian experimental (LP: #320045), Ubuntu remaining changes:
  - Add patch correcting header file locations.
  - Add libvorbis-dev and libgsm1-dev to Build-Depends.
  - Use avcodec_decode_audio2() in source/blender/src/hddaudio.c

[ Alessio Treglia ]
* Add missing previous changelog entries.

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/fftw/reodft/redft00e-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: redft00e-r2hc.c,v 1.32 2006-01-27 02:10:50 athena Exp $ */
 
22
 
 
23
/* Do a REDFT00 problem via an R2HC problem, with some pre/post-processing.
 
24
 
 
25
   This code uses the trick from FFTPACK, also documented in a similar
 
26
   form by Numerical Recipes.  Unfortunately, this algorithm seems to
 
27
   have intrinsic numerical problems (similar to those in
 
28
   reodft11e-r2hc.c), possibly due to the fact that it multiplies its
 
29
   input by a cosine, causing a loss of precision near the zero.  For
 
30
   transforms of 16k points, it has already lost three or four decimal
 
31
   places of accuracy, which we deem unacceptable.
 
32
 
 
33
   So, we have abandoned this algorithm in favor of the one in
 
34
   redft00-r2hc-pad.c, which unfortunately sacrifices 30-50% in speed.
 
35
   The only other alternative in the literature that does not have
 
36
   similar numerical difficulties seems to be the direct adaptation of
 
37
   the Cooley-Tukey decomposition for symmetric data, but this would
 
38
   require a whole new set of codelets and it's not clear that it's
 
39
   worth it at this point.  However, we did implement the latter
 
40
   algorithm for the specific case of odd n (logically adapting the
 
41
   split-radix algorithm); see reodft00e-splitradix.c. */
 
42
 
 
43
#include "reodft.h"
 
44
 
 
45
typedef struct {
 
46
     solver super;
 
47
} S;
 
48
 
 
49
typedef struct {
 
50
     plan_rdft super;
 
51
     plan *cld;
 
52
     twid *td;
 
53
     INT is, os;
 
54
     INT n;
 
55
     INT vl;
 
56
     INT ivs, ovs;
 
57
} P;
 
58
 
 
59
static void apply(const plan *ego_, R *I, R *O)
 
60
{
 
61
     const P *ego = (const P *) ego_;
 
62
     INT is = ego->is, os = ego->os;
 
63
     INT i, n = ego->n;
 
64
     INT iv, vl = ego->vl;
 
65
     INT ivs = ego->ivs, ovs = ego->ovs;
 
66
     R *W = ego->td->W;
 
67
     R *buf;
 
68
     E csum;
 
69
 
 
70
     buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
71
 
 
72
     for (iv = 0; iv < vl; ++iv, I += ivs, O += ovs) {
 
73
          buf[0] = I[0] + I[is * n];
 
74
          csum = I[0] - I[is * n];
 
75
          for (i = 1; i < n - i; ++i) {
 
76
               E a, b, apb, amb;
 
77
               a = I[is * i];
 
78
               b = I[is * (n - i)];
 
79
               csum += W[2*i] * (amb = K(2.0)*(a - b));
 
80
               amb = W[2*i+1] * amb;
 
81
               apb = (a + b);
 
82
               buf[i] = apb - amb;
 
83
               buf[n - i] = apb + amb;
 
84
          }
 
85
          if (i == n - i) {
 
86
               buf[i] = K(2.0) * I[is * i];
 
87
          }
 
88
          
 
89
          {
 
90
               plan_rdft *cld = (plan_rdft *) ego->cld;
 
91
               cld->apply((plan *) cld, buf, buf);
 
92
          }
 
93
          
 
94
          /* FIXME: use recursive/cascade summation for better stability? */
 
95
          O[0] = buf[0];
 
96
          O[os] = csum;
 
97
          for (i = 1; i + i < n; ++i) {
 
98
               INT k = i + i;
 
99
               O[os * k] = buf[i];
 
100
               O[os * (k + 1)] = O[os * (k - 1)] - buf[n - i];
 
101
          }
 
102
          if (i + i == n) {
 
103
               O[os * n] = buf[i];
 
104
          }
 
105
     }
 
106
 
 
107
     X(ifree)(buf);
 
108
}
 
109
 
 
110
static void awake(plan *ego_, enum wakefulness wakefulness)
 
111
{
 
112
     P *ego = (P *) ego_;
 
113
     static const tw_instr redft00e_tw[] = {
 
114
          { TW_COS, 0, 1 },
 
115
          { TW_SIN, 0, 1 },
 
116
          { TW_NEXT, 1, 0 }
 
117
     };
 
118
 
 
119
     X(plan_awake)(ego->cld, wakefulness);
 
120
     X(twiddle_awake)(wakefulness,
 
121
                      &ego->td, redft00e_tw, 2*ego->n, 1, (ego->n+1)/2);
 
122
}
 
123
 
 
124
static void destroy(plan *ego_)
 
125
{
 
126
     P *ego = (P *) ego_;
 
127
     X(plan_destroy_internal)(ego->cld);
 
128
}
 
129
 
 
130
static void print(const plan *ego_, printer *p)
 
131
{
 
132
     const P *ego = (const P *) ego_;
 
133
     p->print(p, "(redft00e-r2hc-%D%v%(%p%))", ego->n + 1, ego->vl, ego->cld);
 
134
}
 
135
 
 
136
static int applicable0(const solver *ego_, const problem *p_)
 
137
{
 
138
     const problem_rdft *p = (const problem_rdft *) p_;
 
139
     UNUSED(ego_);
 
140
 
 
141
     return (1
 
142
             && p->sz->rnk == 1
 
143
             && p->vecsz->rnk <= 1
 
144
             && p->kind[0] == REDFT00
 
145
             && p->sz->dims[0].n > 1  /* n == 1 is not well-defined */
 
146
          );
 
147
}
 
148
 
 
149
static int applicable(const solver *ego, const problem *p, const planner *plnr)
 
150
{
 
151
     return (!NO_SLOWP(plnr) && applicable0(ego, p));
 
152
}
 
153
 
 
154
static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
 
155
{
 
156
     P *pln;
 
157
     const problem_rdft *p;
 
158
     plan *cld;
 
159
     R *buf;
 
160
     INT n;
 
161
     opcnt ops;
 
162
 
 
163
     static const plan_adt padt = {
 
164
          X(rdft_solve), awake, print, destroy
 
165
     };
 
166
 
 
167
     if (!applicable(ego_, p_, plnr))
 
168
          return (plan *)0;
 
169
 
 
170
     p = (const problem_rdft *) p_;
 
171
 
 
172
     n = p->sz->dims[0].n - 1;
 
173
     A(n > 0);
 
174
     buf = (R *) MALLOC(sizeof(R) * n, BUFFERS);
 
175
 
 
176
     cld = X(mkplan_d)(plnr, X(mkproblem_rdft_1_d)(X(mktensor_1d)(n, 1, 1), 
 
177
                                                   X(mktensor_0d)(), 
 
178
                                                   buf, buf, R2HC));
 
179
     X(ifree)(buf);
 
180
     if (!cld)
 
181
          return (plan *)0;
 
182
 
 
183
     pln = MKPLAN_RDFT(P, &padt, apply);
 
184
 
 
185
     pln->n = n;
 
186
     pln->is = p->sz->dims[0].is;
 
187
     pln->os = p->sz->dims[0].os;
 
188
     pln->cld = cld;
 
189
     pln->td = 0;
 
190
 
 
191
     X(tensor_tornk1)(p->vecsz, &pln->vl, &pln->ivs, &pln->ovs);
 
192
     
 
193
     X(ops_zero)(&ops);
 
194
     ops.other = 8 + (n-1)/2 * 11 + (1 - n % 2) * 5;
 
195
     ops.add = 2 + (n-1)/2 * 5;
 
196
     ops.mul = (n-1)/2 * 3 + (1 - n % 2) * 1;
 
197
 
 
198
     X(ops_zero)(&pln->super.super.ops);
 
199
     X(ops_madd2)(pln->vl, &ops, &pln->super.super.ops);
 
200
     X(ops_madd2)(pln->vl, &cld->ops, &pln->super.super.ops);
 
201
 
 
202
     return &(pln->super.super);
 
203
}
 
204
 
 
205
/* constructor */
 
206
static solver *mksolver(void)
 
207
{
 
208
     static const solver_adt sadt = { PROBLEM_RDFT, mkplan };
 
209
     S *slv = MKSOLVER(S, &sadt);
 
210
     return &(slv->super);
 
211
}
 
212
 
 
213
void X(redft00e_r2hc_register)(planner *p)
 
214
{
 
215
     REGISTER_SOLVER(p, mksolver());
 
216
}