test_fftw.C 4.49 KB
Newer Older
Francois Gygi committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
////////////////////////////////////////////////////////////////////////////////
//
// test_fftw.C
//
////////////////////////////////////////////////////////////////////////////////

#include "Timer.h"

#include <iostream>
#include <complex>
#include <valarray>
using namespace std;
#include <cassert>

#include "fftw.h"

Francois Gygi committed
17 18 19 20 21 22 23 24 25
#ifdef IA32
#include "readTSC.h"
long long clk, clk_bwd, clk_fwd;
#endif

#if USE_APC
#include "apc.h"
#endif

Francois Gygi committed
26 27
int main(int argc, char**argv)
{
Francois Gygi committed
28 29 30 31
#if USE_APC
  ApcInit();
#endif

Francois Gygi committed
32
  const int niter = 10;
33 34
  const int np = atoi(argv[1]);
  const int nvec = atoi(argv[2]);
Francois Gygi committed
35
  const int ldz = np + 1;
Francois Gygi committed
36

37
  fftw_plan fwplan, bwplan;
Francois Gygi committed
38 39

  // resize array zvec holding columns
40
  valarray<complex<double> > zvec(nvec * ldz);
Francois Gygi committed
41
  //cout << "zvec ptr: " << &zvec[0] << endl;
Francois Gygi committed
42 43 44
  
  // initialization of FFT libs

45
// #define FFTWMEASURE 1
Francois Gygi committed
46 47
#if FFTWMEASURE
  // FFTWMEASURE
48 49
  fwplan = fftw_create_plan(np,FFTW_FORWARD,FFTW_MEASURE|FFTW_IN_PLACE);
  bwplan = fftw_create_plan(np,FFTW_BACKWARD,FFTW_MEASURE|FFTW_IN_PLACE);
Francois Gygi committed
50 51
#else
  // FFTW_ESTIMATE
52 53
  fwplan = fftw_create_plan(np,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_IN_PLACE);
  bwplan = fftw_create_plan(np,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_IN_PLACE);
Francois Gygi committed
54 55 56 57
#endif

  Timer t_fwd,t_bwd;

Francois Gygi committed
58
#ifdef IA32
Francois Gygi committed
59 60
  clk_bwd = 0;
  clk_fwd = 0;
Francois Gygi committed
61
#endif
Francois Gygi committed
62 63 64
  for ( int iter = 0; iter < niter; iter++ )
  {
  t_bwd.start();
Francois Gygi committed
65

Francois Gygi committed
66 67 68 69 70
   /* 
    * void fftw(fftw_plan plan, int howmany,
    *    FFTW_COMPLEX *in, int istride, int idist,
    *    FFTW_COMPLEX *out, int ostride, int odist);
    */
71
  int ntrans = nvec;
Francois Gygi committed
72 73
  int inc1 = 1;
  int inc2 = ldz;
Francois Gygi committed
74
#ifdef IA32
Francois Gygi committed
75
  clk = readTSC();
Francois Gygi committed
76 77 78 79
#endif
#if USE_APC
  ApcStart(1);
#endif
80
  fftw(bwplan,ntrans,(FFTW_COMPLEX*)&zvec[0],inc1,inc2,
Francois Gygi committed
81
                      (FFTW_COMPLEX*)0,0,0);
Francois Gygi committed
82 83 84 85
#if USE_APC
  ApcStop(1);
#endif
#ifdef IA32
Francois Gygi committed
86
  clk_bwd += readTSC() - clk;
Francois Gygi committed
87
#endif
Francois Gygi committed
88 89
  t_bwd.stop();
  t_fwd.start();
Francois Gygi committed
90
#ifdef IA32
Francois Gygi committed
91
  clk = readTSC();
Francois Gygi committed
92 93 94 95
#endif
#if USE_APC
  ApcStart(2);
#endif
96
  fftw(fwplan,ntrans,(FFTW_COMPLEX*)&zvec[0],inc1,inc2,
Francois Gygi committed
97
                      (FFTW_COMPLEX*)0,0,0);
Francois Gygi committed
98 99 100 101
#if USE_APC
  ApcStop(2);
#endif
#ifdef IA32
Francois Gygi committed
102
  clk_fwd += readTSC() - clk;
Francois Gygi committed
103
#endif
Francois Gygi committed
104 105 106
  t_fwd.stop();
  }

107 108
  fftw_destroy_plan(fwplan);
  fftw_destroy_plan(bwplan);
Francois Gygi committed
109

110 111 112 113
  cout << " fwd: time per transform (in-place,generic)" 
#if FFTWMEASURE
       << "(fftw-measure)"
#endif
Francois Gygi committed
114 115 116 117 118
       << ": " << 1.e6*t_fwd.real()/(niter*nvec) << " microseconds"
#ifdef IA32
       << "  " << clk_fwd/(niter*nvec) << " cycles"
#endif
       << endl;
119 120 121 122 123

  cout << " bwd: time per transform (in-place,generic)" 
#if FFTWMEASURE
       << "(fftw-measure)"
#endif
Francois Gygi committed
124 125 126 127 128
       << ": " << 1.e6*t_bwd.real()/(niter*nvec) << " microseconds"
#ifdef IA32
       << "  " << clk_bwd/(niter*nvec) << " cycles"
#endif
       << endl;
129

130 131 132 133 134 135 136 137 138 139 140 141 142 143
#if 1
  // Use out-of-place, specific plan
  
  valarray<complex<double> > zvec_out(zvec.size());
  t_bwd.reset();
  t_fwd.reset();

  fwplan = fftw_create_plan_specific(np,
    FFTW_FORWARD,FFTW_ESTIMATE|FFTW_OUT_OF_PLACE,
    (FFTW_COMPLEX*)&zvec[0],1,(FFTW_COMPLEX*)&zvec_out[0],1);
  bwplan = fftw_create_plan_specific(np,
    FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_OUT_OF_PLACE,
    (FFTW_COMPLEX*)&zvec[0],1,(FFTW_COMPLEX*)&zvec_out[0],1);
    
Francois Gygi committed
144
#ifdef IA32
Francois Gygi committed
145 146
  clk_bwd = 0;
  clk_fwd = 0;
Francois Gygi committed
147
#endif
148 149 150 151 152 153 154
  for ( int iter = 0; iter < niter; iter++ )
  {

    int ntrans = nvec;
    int inc1 = 1;
    int inc2 = ldz;
    t_bwd.start();
Francois Gygi committed
155 156 157 158 159 160
#ifdef IA32
  clk = readTSC();
#endif
#if USE_APC
  ApcStart(3);
#endif
161 162
    fftw(bwplan,ntrans,(FFTW_COMPLEX*)&zvec[0],inc1,inc2,
                       (FFTW_COMPLEX*)&zvec_out[0],inc1,inc2);
Francois Gygi committed
163 164 165 166 167 168
#if USE_APC
  ApcStop(3);
#endif
#ifdef IA32
  clk_bwd += readTSC() - clk;
#endif
169 170 171
    t_bwd.stop();
  
    t_fwd.start();
Francois Gygi committed
172 173 174 175 176 177
#ifdef IA32
  clk = readTSC();
#endif
#if USE_APC
  ApcStart(4);
#endif
178 179
    fftw(fwplan,ntrans,(FFTW_COMPLEX*)&zvec[0],inc1,inc2,
                       (FFTW_COMPLEX*)&zvec_out[0],inc1,inc2);
Francois Gygi committed
180 181 182 183 184 185
#if USE_APC
  ApcStop(4);
#endif
#ifdef IA32
  clk_fwd += readTSC() - clk;
#endif
186 187 188
    t_fwd.stop();

  }
Francois Gygi committed
189

190 191 192
  fftw_destroy_plan(fwplan);
  fftw_destroy_plan(bwplan);
  
193 194 195 196
  cout << " fwd: time per transform (out-of-place,specific)" 
#if FFTWMEASURE
       << "(fftw-measure)"
#endif
Francois Gygi committed
197 198 199 200 201
       << ": " << 1.e6*t_fwd.real()/(niter*nvec) << " microseconds"
#ifdef IA32
       << "  " << clk_fwd/(niter*nvec) << " cycles"
#endif
       << endl;
202 203 204 205 206

  cout << " bwd: time per transform (out-of-place,specific)" 
#if FFTWMEASURE
       << "(fftw-measure)"
#endif
Francois Gygi committed
207 208 209 210 211 212 213 214
       << ": " << 1.e6*t_bwd.real()/(niter*nvec) << " microseconds"
#ifdef IA32
       << "  " << clk_bwd/(niter*nvec) << " cycles"
#endif
       << endl;
#endif
#if USE_APC
  ApcFinalize();
215
#endif
Francois Gygi committed
216 217
  return 0;
}