test_fftw.C 4.59 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17 18 19 20 21 22 23 24 25 26 27 28
// test_fftw.C
//
////////////////////////////////////////////////////////////////////////////////

#include "Timer.h"

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

#include "fftw.h"

Francois Gygi committed
29 30 31 32 33
#ifdef IA32
#include "readTSC.h"
long long clk, clk_bwd, clk_fwd;
#endif

Francois Gygi committed
34 35 36
int main(int argc, char**argv)
{
  const int niter = 10;
37 38
  const int np = atoi(argv[1]);
  const int nvec = atoi(argv[2]);
Francois Gygi committed
39
  const int ldz = np + 1;
Francois Gygi committed
40

41
  fftw_plan fwplan, bwplan;
Francois Gygi committed
42 43

  // resize array zvec holding columns
44
  valarray<complex<double> > zvec(nvec * ldz);
Francois Gygi committed
45
  //cout << "zvec ptr: " << &zvec[0] << endl;
46

Francois Gygi committed
47 48
  // initialization of FFT libs

49
// #define FFTWMEASURE 1
Francois Gygi committed
50 51
#if FFTWMEASURE
  // FFTWMEASURE
52 53
  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
54 55
#else
  // FFTW_ESTIMATE
56 57
  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
58 59 60 61
#endif

  Timer t_fwd,t_bwd;

Francois Gygi committed
62
#ifdef IA32
Francois Gygi committed
63 64
  clk_bwd = 0;
  clk_fwd = 0;
Francois Gygi committed
65
#endif
Francois Gygi committed
66 67 68
  for ( int iter = 0; iter < niter; iter++ )
  {
  t_bwd.start();
Francois Gygi committed
69

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

99 100
  fftw_destroy_plan(fwplan);
  fftw_destroy_plan(bwplan);
Francois Gygi committed
101

102
  cout << " fwd: time per transform (in-place,generic)"
103 104 105
#if FFTWMEASURE
       << "(fftw-measure)"
#endif
Francois Gygi committed
106 107 108 109 110
       << ": " << 1.e6*t_fwd.real()/(niter*nvec) << " microseconds"
#ifdef IA32
       << "  " << clk_fwd/(niter*nvec) << " cycles"
#endif
       << endl;
111

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

122 123
#if 1
  // Use out-of-place, specific plan
124

125 126 127 128 129 130 131 132 133 134
  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);
135

Francois Gygi committed
136
#ifdef IA32
Francois Gygi committed
137 138
  clk_bwd = 0;
  clk_fwd = 0;
Francois Gygi committed
139
#endif
140 141 142 143 144 145 146
  for ( int iter = 0; iter < niter; iter++ )
  {

    int ntrans = nvec;
    int inc1 = 1;
    int inc2 = ldz;
    t_bwd.start();
Francois Gygi committed
147 148 149
#ifdef IA32
  clk = readTSC();
#endif
150 151
    fftw(bwplan,ntrans,(FFTW_COMPLEX*)&zvec[0],inc1,inc2,
                       (FFTW_COMPLEX*)&zvec_out[0],inc1,inc2);
Francois Gygi committed
152 153 154
#ifdef IA32
  clk_bwd += readTSC() - clk;
#endif
155
    t_bwd.stop();
156

157
    t_fwd.start();
Francois Gygi committed
158 159 160
#ifdef IA32
  clk = readTSC();
#endif
161 162
    fftw(fwplan,ntrans,(FFTW_COMPLEX*)&zvec[0],inc1,inc2,
                       (FFTW_COMPLEX*)&zvec_out[0],inc1,inc2);
Francois Gygi committed
163 164 165
#ifdef IA32
  clk_fwd += readTSC() - clk;
#endif
166 167 168
    t_fwd.stop();

  }
Francois Gygi committed
169

170 171
  fftw_destroy_plan(fwplan);
  fftw_destroy_plan(bwplan);
172 173

  cout << " fwd: time per transform (out-of-place,specific)"
174 175 176
#if FFTWMEASURE
       << "(fftw-measure)"
#endif
Francois Gygi committed
177 178 179 180 181
       << ": " << 1.e6*t_fwd.real()/(niter*nvec) << " microseconds"
#ifdef IA32
       << "  " << clk_fwd/(niter*nvec) << " cycles"
#endif
       << endl;
182

183
  cout << " bwd: time per transform (out-of-place,specific)"
184 185 186
#if FFTWMEASURE
       << "(fftw-measure)"
#endif
Francois Gygi committed
187 188 189 190 191
       << ": " << 1.e6*t_bwd.real()/(niter*nvec) << " microseconds"
#ifdef IA32
       << "  " << clk_bwd/(niter*nvec) << " cycles"
#endif
       << endl;
192
#endif
Francois Gygi committed
193 194
  return 0;
}