Changed strategy about ffmpeg fft usage. Using system ffmpeg rather than bundled...
[gsoc2010-fftw-neon:gsoc2010-fftw-neon.git] / dft / simd / nonportable / arm / neon-example.c
1 #include "dft.h"
2
3 /*
4  * Do not allow these files to improve runtimes vs. codelets (my intrinsics 
5  * implementation is now synonymous with codelets, will fix that later)
6  */
7 #if defined(FFTW_SINGLE) && defined(HAVE_NEON) && !defined(HAVE_NEON_INTRINSICS)
8
9 #ifndef _MATH_H
10 #include <math.h>
11 #endif
12 #ifndef _COMPLEX_H
13 #include <complex.h>
14 #endif
15
16 typedef struct {
17      solver super;
18 } S;
19
20 typedef struct {
21      plan_dft super;
22      INT n;     /* problem size */
23      INT m;    /* m = log2(n) */
24      R *w;      /* input */
25      R *W;      /* output */
26      plan *p;
27 } P;
28
29 static int _log2(unsigned int n) {
30         unsigned int m = 0;
31         for( m = 0; m < 32; m++ )
32                 if ( ((1<<m)^n) == 0 )
33                         break;
34         return ( m == 32 ) ? -1 : m;
35 }
36
37 // simple, recursive fft radix 2
38 typedef float f32;
39 typedef _Complex float cf32;
40 typedef unsigned int u32;
41 static void ditfft2( cf32 *X, cf32 *Y, u32 N, u32 s ) {
42   f32 th;
43   cf32 t;
44   u32 k;
45   const f32 pi = M_PI;
46   const cf32 j = _Complex_I;
47   if ( N == 1 ) {
48     Y[0] = X[0];
49   } else {
50     ditfft2( X, Y, N/2, 2*s );
51     ditfft2( &X[s], &Y[N/2], N/2, 2*s );
52     for ( k=0; k<N/2; k++ ) {
53       t = Y[k];
54       th = -2*pi*k/((float)N);
55       Y[k] = t + ( cosf(th) + j*sinf(th) ) * Y[k+N/2];
56       Y[k+N/2] = t - ( cosf(th) + j*sinf(th) ) * Y[k+N/2];
57     }
58   }
59 }
60
61 static void apply(const plan *ego_, R *ri, R *ii, R *ro, R *io)
62 {
63      const P *ego = (const P *) ego_;
64      INT n = ego->n;
65      R *w = (R *)ego->w, *W = (R *) ego->W;
66
67      ditfft2((cf32 *)w,(cf32 *)W,n,1);
68 }
69
70 static void awake(plan *ego_, enum wakefulness wakefulness)
71 {
72         (void) ego_; /* UNUSED */
73         (void) wakefulness; /* UNUSED */
74 }
75
76 static int applicable0(const problem *p_)
77 {
78      const problem_dft *p = (const problem_dft *) p_;
79      return (1
80              && p->sz->rnk == 1 /* 1d data */
81              && p->vecsz->rnk == 0 /* a single vector */
82              && _log2(p->sz->dims[0].n) != -1 /* power of two */
83              && p->io == &p->ro[1] /* contiguous memory */
84              && p->ii == &p->ri[1]
85              && p->sz->dims[0].is == 2
86              && p->sz->dims[0].os == 2
87              && p->ri != p->ro /* out-of-place */
88              && p->ii != p->io
89      );
90 }
91
92 static int applicable(const solver *ego, const problem *p_,
93                       const planner *plnr)
94 {
95      UNUSED(ego);
96      if (!applicable0(p_)) return 0;
97      return 1;
98 }
99
100 static void destroy(plan *ego_)
101 {
102      P *ego = (P *) ego_;
103      X(plan_destroy_internal)(ego->p);
104 }
105
106 static void print(const plan *ego_, printer *p)
107 {
108      const P *ego = (const P *)ego_;
109      p->print(p, "(dft-example-%D)",
110               ego->n);
111 }
112
113 static plan *mkplan(const solver *ego, const problem *p_, planner *plnr)
114 {
115      const problem_dft *p = (const problem_dft *) p_;
116      P *pln;
117      INT n, m;
118      plan *plan = 0;
119
120      static const plan_adt padt = {
121           X(dft_solve), awake, print, destroy
122      };
123
124      if (!applicable(ego, p_, plnr))
125           return NULL;
126
127      n = p->sz->dims[0].n;
128
129      plan = X(mkplan_f_d)(plnr,
130                           X(mkproblem_dft_d)(X(mktensor_1d)(2, 2, 2),
131                                              X(mktensor_1d)(1, 0, 0),
132                                              p->ri, p->ii,
133                                              p->ro, p->io),
134                           0, 0, 0);
135      if (!plan) goto nada;
136
137      pln = MKPLAN_DFT(P, &padt, apply);
138
139      pln->n = n;
140      pln->m = _log2(n);
141      pln->w = 0;
142      pln->W = 0;
143      pln->p = plan;
144
145      X(ops_add)(&plan->ops, &plan->ops, &pln->super.super.ops);
146      pln->super.super.ops.add += 2 * n * m;
147      pln->super.super.ops.mul += 5 * n * n;
148      pln->super.super.ops.other += n * m;
149
150      return &(pln->super.super);
151
152  nada:
153      X(plan_destroy_internal)(plan);
154      return NULL;
155 }
156
157
158 static solver *mksolver(void)
159 {
160      static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 };
161      S *slv = MKSOLVER(S, &sadt);
162      return &(slv->super);
163 }
164
165 void X(dft_neon_example_register)(planner *p)
166 {
167      REGISTER_SOLVER(p, mksolver());
168 }
169
170 #endif