1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009 Mark Borgerding mark a borgerding net
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#include <iostream>
11
12#include <bench/BenchUtil.h>
13#include <complex>
14#include <vector>
15#include <Eigen/Core>
16
17#include <unsupported/Eigen/FFT>
18
19using namespace Eigen;
20using namespace std;
21
22
23template <typename T>
24string nameof();
25
26template <> string nameof<float>() {return "float";}
27template <> string nameof<double>() {return "double";}
28template <> string nameof<long double>() {return "long double";}
29
30#ifndef TYPE
31#define TYPE float
32#endif
33
34#ifndef NFFT
35#define NFFT 1024
36#endif
37#ifndef NDATA
38#define NDATA 1000000
39#endif
40
41using namespace Eigen;
42
43template <typename T>
44void bench(int nfft,bool fwd,bool unscaled=false, bool halfspec=false)
45{
46    typedef typename NumTraits<T>::Real Scalar;
47    typedef typename std::complex<Scalar> Complex;
48    int nits = NDATA/nfft;
49    vector<T> inbuf(nfft);
50    vector<Complex > outbuf(nfft);
51    FFT< Scalar > fft;
52
53    if (unscaled) {
54        fft.SetFlag(fft.Unscaled);
55        cout << "unscaled ";
56    }
57    if (halfspec) {
58        fft.SetFlag(fft.HalfSpectrum);
59        cout << "halfspec ";
60    }
61
62
63    std::fill(inbuf.begin(),inbuf.end(),0);
64    fft.fwd( outbuf , inbuf);
65
66    BenchTimer timer;
67    timer.reset();
68    for (int k=0;k<8;++k) {
69        timer.start();
70        if (fwd)
71            for(int i = 0; i < nits; i++)
72                fft.fwd( outbuf , inbuf);
73        else
74            for(int i = 0; i < nits; i++)
75                fft.inv(inbuf,outbuf);
76        timer.stop();
77    }
78
79    cout << nameof<Scalar>() << " ";
80    double mflops = 5.*nfft*log2((double)nfft) / (1e6 * timer.value() / (double)nits );
81    if ( NumTraits<T>::IsComplex ) {
82        cout << "complex";
83    }else{
84        cout << "real   ";
85        mflops /= 2;
86    }
87
88
89    if (fwd)
90        cout << " fwd";
91    else
92        cout << " inv";
93
94    cout << " NFFT=" << nfft << "  " << (double(1e-6*nfft*nits)/timer.value()) << " MS/s  " << mflops << "MFLOPS\n";
95}
96
97int main(int argc,char ** argv)
98{
99    bench<complex<float> >(NFFT,true);
100    bench<complex<float> >(NFFT,false);
101    bench<float>(NFFT,true);
102    bench<float>(NFFT,false);
103    bench<float>(NFFT,false,true);
104    bench<float>(NFFT,false,true,true);
105
106    bench<complex<double> >(NFFT,true);
107    bench<complex<double> >(NFFT,false);
108    bench<double>(NFFT,true);
109    bench<double>(NFFT,false);
110    bench<complex<long double> >(NFFT,true);
111    bench<complex<long double> >(NFFT,false);
112    bench<long double>(NFFT,true);
113    bench<long double>(NFFT,false);
114    return 0;
115}
116