~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to scipy/fftpack/src/zfft_fftw3.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
 
 
2
// This cache uses FFTW_MEASURE for the plans, and do not copy the data.
 
3
GEN_CACHE(zfftw3,(int n,int d)
 
4
        ,int direction;
 
5
        fftw_plan plan;
 
6
    fftw_complex *wrk;
 
7
        ,((caches_zfftw3[i].n==n) &&
 
8
            (caches_zfftw3[i].direction==d))
 
9
        ,caches_zfftw3[id].direction = d;
 
10
        // This working buffer is only used to compute the plan: we need it
 
11
        // since FFTW_MEASURE destroys its input when computing a plan
 
12
            caches_zfftw3[id].wrk = fftw_malloc(n * sizeof(double) * 2); 
 
13
            caches_zfftw3[id].plan = fftw_plan_dft_1d(n, 
 
14
                caches_zfftw3[id].wrk,
 
15
                caches_zfftw3[id].wrk,
 
16
                (d>0?FFTW_FORWARD:FFTW_BACKWARD),
 
17
                FFTW_ESTIMATE | FFTW_UNALIGNED);
 
18
        ,//fftw_print_plan(caches_zfftw3[id].plan);
 
19
    fftw_destroy_plan(caches_zfftw3[id].plan);
 
20
        fftw_free(caches_zfftw3[id].wrk);
 
21
    //fflush(stdout);
 
22
    //fprintf(stderr, "aligned count %d\n", countaligned);
 
23
        ,10)
 
24
 
 
25
static void zfft_fftw3(complex_double * inout, int n, int dir, int
 
26
howmany, int normalize)
 
27
{
 
28
        fftw_complex    *ptr = (fftw_complex*)inout;
 
29
    fftw_complex    *ptrm;
 
30
        fftw_plan       plan = NULL;
 
31
    double          factor = 1./n;
 
32
 
 
33
        int i;
 
34
 
 
35
        plan = caches_zfftw3[get_cache_id_zfftw3(n, dir)].plan;
 
36
 
 
37
        switch (dir) {
 
38
        case 1:
 
39
                for (i = 0; i < howmany; ++i, ptr += n) {
 
40
                        fftw_execute_dft(plan, ptr, ptr);
 
41
                }
 
42
                break;
 
43
 
 
44
        case -1:
 
45
                for (i = 0; i < howmany; ++i, ptr += n) {
 
46
                        fftw_execute_dft(plan, ptr, ptr);
 
47
                }
 
48
                break;
 
49
 
 
50
        default:
 
51
                fprintf(stderr, "zfft: invalid dir=%d\n", dir);
 
52
        }
 
53
 
 
54
        if (normalize) {
 
55
        ptr =(fftw_complex*)inout;
 
56
                for (i = n * howmany - 1; i >= 0; --i) {
 
57
                        *((double *) (ptr)) *= factor;
 
58
                        *((double *) (ptr++) + 1) *= factor;
 
59
                        //*((double *) (ptr)) /= n;
 
60
                        //*((double *) (ptr++) + 1) /= n;
 
61
                }
 
62
        }
 
63
}
 
64
#if 0
 
65
GEN_CACHE(zfftw3,(int n,int d)
 
66
        ,int direction;
 
67
        fftw_plan plan;
 
68
        fftw_complex* ptr;
 
69
        ,((caches_zfftw3[i].n==n) &&
 
70
            (caches_zfftw3[i].direction==d))
 
71
        ,caches_zfftw3[id].direction = d;
 
72
        caches_zfftw3[id].ptr = fftw_malloc(sizeof(fftw_complex)*(n));
 
73
            caches_zfftw3[id].plan = fftw_plan_dft_1d(n, caches_zfftw3[id].ptr,
 
74
        caches_zfftw3[id].ptr,
 
75
                (d>0?FFTW_FORWARD:FFTW_BACKWARD),
 
76
                FFTW_ESTIMATE);
 
77
        ,fftw_destroy_plan(caches_zfftw3[id].plan);
 
78
        fftw_free(caches_zfftw3[id].ptr);
 
79
        ,10)
 
80
 
 
81
static void zfft_fftw3(complex_double * inout, int n, int dir, int
 
82
howmany, int normalize)
 
83
{
 
84
        complex_double *ptr = inout;
 
85
        fftw_complex *ptrm = NULL;
 
86
        fftw_plan plan = NULL;
 
87
 
 
88
        int i;
 
89
 
 
90
        plan = caches_zfftw3[get_cache_id_zfftw3(n, dir)].plan;
 
91
 
 
92
        switch (dir) {
 
93
        case 1:
 
94
                for (i = 0; i < howmany; ++i, ptr += n) {
 
95
                        ptrm =
 
96
                            caches_zfftw3[get_cache_id_zfftw3(n, dir)].ptr;
 
97
                        memcpy(ptrm, ptr, sizeof(double) * 2 * n);
 
98
                        fftw_execute(plan);
 
99
                        memcpy(ptr, ptrm, sizeof(double) * 2 * n);
 
100
                }
 
101
                break;
 
102
 
 
103
        case -1:
 
104
                for (i = 0; i < howmany; ++i, ptr += n) {
 
105
                        ptrm =
 
106
                            caches_zfftw3[get_cache_id_zfftw3(n, dir)].ptr;
 
107
                        memcpy(ptrm, ptr, sizeof(double) * 2 * n);
 
108
                        fftw_execute(plan);
 
109
                        memcpy(ptr, ptrm, sizeof(double) * 2 * n);
 
110
                }
 
111
                break;
 
112
 
 
113
        default:
 
114
                fprintf(stderr, "zfft: invalid dir=%d\n", dir);
 
115
        }
 
116
 
 
117
        if (normalize) {
 
118
                ptr = inout;
 
119
                for (i = n * howmany - 1; i >= 0; --i) {
 
120
                        *((double *) (ptr)) /= n;
 
121
                        *((double *) (ptr++) + 1) /= n;
 
122
                }
 
123
        }
 
124
}
 
125
#endif