← Back to branch summary

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

« back to all changes in this revision

Viewing changes to extern/fftw/rdft/dft-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/rdft/dft-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: dft-r2hc.c,v 1.37 2006-01-27 02:10:50 athena Exp $ */
 
22
 
 
23
/* Compute the complex DFT by combining R2HC RDFTs on the real
 
24
   and imaginary parts.   This could be useful for people just wanting
 
25
   to link to the real codelets and not the complex ones.  It could
 
26
   also even be faster than the complex algorithms for split (as opposed
 
27
   to interleaved) real/imag complex data. */
 
28
 
 
29
#include "rdft.h"
 
30
#include "dft.h"
 
31
 
 
32
typedef struct {
 
33
     solver super;
 
34
} S;
 
35
 
 
36
typedef struct {
 
37
     plan_dft super;
 
38
     plan *cld;
 
39
     INT ishift, oshift;
 
40
     INT os;
 
41
     INT n;
 
42
} P;
 
43
 
 
44
static void apply(const plan *ego_, R *ri, R *ii, R *ro, R *io)
 
45
{
 
46
     const P *ego = (const P *) ego_;
 
47
     INT n;
 
48
 
 
49
     UNUSED(ii);
 
50
 
 
51
     { /* transform vector of real & imag parts: */
 
52
          plan_rdft *cld = (plan_rdft *) ego->cld;
 
53
          cld->apply((plan *) cld, ri + ego->ishift, ro + ego->oshift);
 
54
     }
 
55
 
 
56
     n = ego->n;
 
57
     if (n > 1) {
 
58
          INT i, os = ego->os;
 
59
          for (i = 1; i < (n + 1)/2; ++i) {
 
60
               E rop, iop, iom, rom;
 
61
               rop = ro[os * i];
 
62
               iop = io[os * i];
 
63
               rom = ro[os * (n - i)];
 
64
               iom = io[os * (n - i)];
 
65
               ro[os * i] = rop - iom;
 
66
               io[os * i] = iop + rom;
 
67
               ro[os * (n - i)] = rop + iom;
 
68
               io[os * (n - i)] = iop - rom;
 
69
          }
 
70
     }
 
71
}
 
72
 
 
73
static void awake(plan *ego_, enum wakefulness wakefulness)
 
74
{
 
75
     P *ego = (P *) ego_;
 
76
     X(plan_awake)(ego->cld, wakefulness);
 
77
}
 
78
 
 
79
static void destroy(plan *ego_)
 
80
{
 
81
     P *ego = (P *) ego_;
 
82
     X(plan_destroy_internal)(ego->cld);
 
83
}
 
84
 
 
85
static void print(const plan *ego_, printer *p)
 
86
{
 
87
     const P *ego = (const P *) ego_;
 
88
     p->print(p, "(dft-r2hc-%D%(%p%))", ego->n, ego->cld);
 
89
}
 
90
 
 
91
 
 
92
static int applicable0(const problem *p_)
 
93
{
 
94
     const problem_dft *p = (const problem_dft *) p_;
 
95
     return ((p->sz->rnk == 1 && p->vecsz->rnk == 0)
 
96
             || (p->sz->rnk == 0 && FINITE_RNK(p->vecsz->rnk))
 
97
          );
 
98
}
 
99
 
 
100
static int splitp(R *r, R *i, INT n, INT s)
 
101
{
 
102
     return ((r > i ? (r - i) : (i - r)) >= n * (s > 0 ? s : 0-s));
 
103
}
 
104
 
 
105
static int applicable(const problem *p_, const planner *plnr)
 
106
{
 
107
     if (!applicable0(p_)) return 0;
 
108
 
 
109
     {
 
110
          const problem_dft *p = (const problem_dft *) p_;
 
111
 
 
112
          /* rank-0 problems are always OK */
 
113
          if (p->sz->rnk == 0) return 1;
 
114
 
 
115
          /* this solver is ok for split arrays */
 
116
          if (p->sz->rnk == 1 &&
 
117
              splitp(p->ri, p->ii, p->sz->dims[0].n, p->sz->dims[0].is) &&
 
118
              splitp(p->ro, p->io, p->sz->dims[0].n, p->sz->dims[0].os))
 
119
               return 1;
 
120
 
 
121
          return !(NO_DFT_R2HCP(plnr));
 
122
     }
 
123
}
 
124
 
 
125
static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
 
126
{
 
127
     P *pln;
 
128
     const problem_dft *p;
 
129
     plan *cld;
 
130
     INT ishift = 0, oshift = 0;
 
131
 
 
132
     static const plan_adt padt = {
 
133
          X(dft_solve), awake, print, destroy
 
134
     };
 
135
 
 
136
     UNUSED(ego_);
 
137
     if (!applicable(p_, plnr))
 
138
          return (plan *)0;
 
139
 
 
140
     p = (const problem_dft *) p_;
 
141
 
 
142
     {
 
143
          tensor *ri_vec = X(mktensor_1d)(2, p->ii - p->ri, p->io - p->ro);
 
144
          tensor *cld_vec = X(tensor_append)(ri_vec, p->vecsz);
 
145
          int i;
 
146
          for (i = 0; i < cld_vec->rnk; ++i) { /* make all istrides > 0 */
 
147
               if (cld_vec->dims[i].is < 0) {
 
148
                    INT nm1 = cld_vec->dims[i].n - 1;
 
149
                    ishift -= nm1 * (cld_vec->dims[i].is *= -1);
 
150
                    oshift -= nm1 * (cld_vec->dims[i].os *= -1);
 
151
               }
 
152
          }
 
153
          cld = X(mkplan_d)(plnr, 
 
154
                            X(mkproblem_rdft_1)(p->sz, cld_vec, 
 
155
                                                p->ri + ishift, 
 
156
                                                p->ro + oshift, R2HC));
 
157
          X(tensor_destroy2)(ri_vec, cld_vec);
 
158
     }
 
159
     if (!cld) return (plan *)0;
 
160
 
 
161
     pln = MKPLAN_DFT(P, &padt, apply);
 
162
 
 
163
     if (p->sz->rnk == 0) {
 
164
          pln->n = 1;
 
165
          pln->os = 0;
 
166
     }
 
167
     else {
 
168
          pln->n = p->sz->dims[0].n;
 
169
          pln->os = p->sz->dims[0].os;
 
170
     }
 
171
     pln->ishift = ishift;
 
172
     pln->oshift = oshift;
 
173
 
 
174
     pln->cld = cld;
 
175
     
 
176
     pln->super.super.ops = cld->ops;
 
177
     pln->super.super.ops.other += 8 * ((pln->n - 1)/2);
 
178
     pln->super.super.ops.add += 4 * ((pln->n - 1)/2);
 
179
     pln->super.super.ops.other += 1; /* estimator hack for nop plans */
 
180
 
 
181
     return &(pln->super.super);
 
182
}
 
183
 
 
184
/* constructor */
 
185
static solver *mksolver(void)
 
186
{
 
187
     static const solver_adt sadt = { PROBLEM_DFT, mkplan };
 
188
     S *slv = MKSOLVER(S, &sadt);
 
189
     return &(slv->super);
 
190
}
 
191
 
 
192
void X(dft_r2hc_register)(planner *p)
 
193
{
 
194
     REGISTER_SOLVER(p, mksolver());
 
195
}