testFFTW.cpp 4.98 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
// ===================================================================================
// Ce LOGICIEL "ScalFmm" est couvert par le copyright Inria 20xx-2012.
// Inria détient tous les droits de propriété sur le LOGICIEL, et souhaite que
// la communauté scientifique l'utilise afin de le tester et de l'évaluer.
// Inria donne gracieusement le droit d'utiliser ce LOGICIEL. Toute utilisation
// dans un but lucratif ou à des fins commerciales est interdite sauf autorisation
// expresse et préalable d'Inria.
// Toute utilisation hors des limites précisées ci-dessus et réalisée sans l'accord
// expresse préalable d'Inria constituerait donc le délit de contrefaçon.
// Le LOGICIEL étant un produit en cours de développement, Inria ne saurait assurer
// aucune responsabilité et notamment en aucune manière et en aucun cas, être tenu
// de répondre d'éventuels dommages directs ou indirects subits par l'utilisateur.
// Tout utilisateur du LOGICIEL s'engage à communiquer à Inria ses remarques
// relatives à l'usage du LOGICIEL
// ===================================================================================

17 18 19 20
// ==== CMAKE =====
// @FUSE_FFT
// ================

21 22 23
#include <iostream>
#include <stdlib.h>

24
#include <fftw3.h>
25 26

#include "../../Src/Utils/FGlobal.hpp"
27
#include "../../Src/Utils/FComplex.hpp"
28 29 30

#include "../../Src/Utils/FTic.hpp"

31
#include "../../Src/Utils/FParameterNames.hpp"
32 33 34 35

//#include "../../Src/Utils/FDft.hpp"


36
int main(int argc, char** argv)
37
{
38 39
    FHelpDescribeAndExit(argc, argv, "Test the FFTw (only the code is interesting).");

40
    typedef double FReal;
41 42 43 44 45 46 47 48 49 50 51 52
    const FReal FRandMax = FReal(RAND_MAX);

    FTic time;

    //////////////////////////////////////////////////////////////////////////////
    // INITIALIZATION

    // size (pick a power of 2 for better performance of the FFT algorithm)
    unsigned int nsteps_ = 500;

    // fftw arrays
    FReal* fftR_;
53
    FComplex<FReal>* fftC_;
54 55

    fftR_ = (FReal*) fftw_malloc(sizeof(FReal) * nsteps_);
56
    fftC_ = (FComplex<FReal>*) fftw_malloc(sizeof(FComplex<FReal>) * nsteps_);
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111

    // fftw plans
    // use routine defined in file:
    // /PATH/TO/mkl/interfaces/fftw3xf/wrappers/fftw_plan_dft_c2r_1d.c
    fftw_plan plan_c2r_; // backward FFT plan
    fftw_plan plan_r2c_; // forward FFT plan

    std::cout<< "Init FFTW plans: ";
    time.tic();

    plan_c2r_ =
            fftw_plan_dft_c2r_1d(nsteps_,
                                 reinterpret_cast<fftw_complex*>(fftC_),
                                 fftR_,
                                 FFTW_MEASURE);
    plan_r2c_ =
            fftw_plan_dft_r2c_1d(nsteps_,
                                 fftR_,
                                 reinterpret_cast<fftw_complex*>(fftC_),
                                 FFTW_MEASURE);

    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    // alternative choice of plan type
    //  plan_c2r_ =
    //    fftw_plan_dft_1d(nsteps_, fftC_, fftR_, FFTW_BACKWARD, FFTW_MEASURE);
    //  plan_r2c_ =
    //    fftw_plan_dft_1d(nsteps_, fftR_, fftC_, FFTW_FORWARD, FFTW_MEASURE);

    //////////////////////////////////////////////////////////////////////////////
    // EXECUTION
    // generate random physical data
    for(unsigned int s=0; s<nsteps_; ++s)
        fftR_[s] = FReal(rand())/FRandMax;

    //  // display data in  physical space
    //  std::cout<< "Physical data: "<<std::endl;
    //  for(unsigned int s=0; s<nsteps_; ++s)
    //    std::cout<< fftR_[s] << ", ";
    //  std::cout<<std::endl;

    // perform fft
    std::cout<< "Perform Forward FFT: ";
    time.tic();
    fftw_execute( plan_r2c_ );
    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    //  // display transform in Fourier space
    //  // beware the real data FFT stores only N/2+1 complex output values
    //  std::cout<< "Transformed data : "<<std::endl;
    //  for(unsigned int s=0; s<nsteps_; ++s)
    //    std::cout<< fftC_[s] << ", ";
    //  std::cout<<std::endl;

    //  for(unsigned int s=0; s<nsteps_/2+1; ++s){
112
    //    fftC_[nsteps_-s]=FComplex<FReal>(fftC_[s].getReal(),-fftC_[s].getImag());
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
    //  }
    //
    //  std::cout<< "Full Transformed data : "<<std::endl;
    //  for(unsigned int s=0; s<nsteps_; ++s)
    //    std::cout<< fftC_[s] << ", ";
    //  std::cout<<std::endl;

    // perform ifft of tranformed data (in order to get physical values back)
    std::cout<< "Perform Backward FFT: ";
    time.tic();
    fftw_execute( plan_c2r_ );
    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    //  // display data in physical space
    //  std::cout<< "Physical data (from 1/N*IFFT(FFT(Physical data))): "<<std::endl;
    //  for(unsigned int s=0; s<nsteps_; ++s)
    //    std::cout<< fftR_[s]/nsteps_ << ", ";
    //  std::cout<<std::endl;

    //////////////////////////////////////////////////////////////////////////////
    // VALIDATION
    // TODO

    //free memory
    fftw_destroy_plan(plan_r2c_);
    fftw_destroy_plan(plan_c2r_);
    fftw_free(fftR_);
    fftw_free(fftC_);
141 142 143


}// end test