testFastDiscreteConvolution.cpp 5.77 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
// ===================================================================================
// 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
// ===================================================================================

// ==== CMAKE =====
18
// @FUSE_FFT
19
// @FUSE_BLAS
20
// ================
21 22 23
// Keep in private GIT
// @SCALFMM_PRIVATE

24 25 26 27 28 29 30


#include <iostream>

#include <stdlib.h>
#include <time.h>

31 32 33
#include "Utils/FTic.hpp"
#include "Utils/FMath.hpp"
#include "Utils/FBlas.hpp" 
34 35

// 
36
#include "Utils/stdComplex.hpp"
37
#include "Utils/FDft.hpp"
38

39
#include "Utils/FParameterNames.hpp"
40 41


42 43
int main(int argc, char ** argv){
    FHelpDescribeAndExit(argc, argv, "Test the FFT (only the code is interesting).");
44

45
    typedef double FReal;
46
    const FReal FRandMax = FReal(RAND_MAX);
47

48
    FTic time;
49

50 51
    const unsigned int ORDER = 3;
    const unsigned int N=(2*ORDER-1);
52 53


54 55 56
    //////////////////////////////////////////////////////////////////////////////
    // Application of Convolution Theorem
    // Let there be a Circulant Toeplitz matrix K of size N
57

58 59 60
    FReal K[N*N];
    FReal Y[N];
    FReal X[N];
61

62 63
    for(unsigned int i=0; i<N; ++i) X[i]=0.0;
    for(unsigned int i=0; i<N; ++i) Y[i]=FReal(rand())/FRandMax;
64

65 66 67 68
    std::cout<< "Y: "<<std::endl;
    for(unsigned int i=0; i<N; ++i)
        std::cout<< Y[i] << ", ";
    std::cout<<std::endl;
69

70 71 72 73 74 75 76 77
    // Assemble Matrix K
    std::cout<< "Assemble Full Matrix K: ";
    time.tic();
    for(unsigned int i=0; i<N; ++i){
        for(unsigned int j=0; j<N; ++j){
            if(i>j) K[i*N+j]=i-j-1;
            else   K[i*N+j]=N+i-j-1;
        }
78
    }
79
    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;
80

81 82 83 84 85
    std::cout<< "Circulant Toeplitz matrix K: "<<std::endl;
    for(unsigned int i=0; i<N; ++i){
        for(unsigned int j=0; j<N; ++j){
            std::cout<< K[i*N+j] << ", ";
        } std::cout<<std::endl;
86
    } std::cout<<std::endl;
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103

    // Direct application of Circulant Matrix
    std::cout<< "Direct application of K: ";
    time.tic();
    for(unsigned int i=0; i<N; ++i){
        for(unsigned int j=0; j<N; ++j){
            X[i]+=K[i*N+j]*Y[j];
        }
    }
    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    std::cout<< "X=KY: "<<std::endl;
    for(unsigned int i=0; i<N; ++i)
        std::cout<< X[i] << ", ";
    std::cout<<std::endl;

    // now compute via DFT and use convolution theorem
104 105 106
    stdComplex<FReal> FK[N];
    stdComplex<FReal> FY[N];
    stdComplex<FReal> FX[N];
107 108 109 110 111 112 113
    FReal iFX[N];

    // Init DFTor
    std::cout<< "Set DFT: ";
    time.tic();
    const int dim = 1;
    //FDft<FReal> Dft(N);// direct version (Beware! Ordering of output differs from REAL valued-FFT)
114
    FFftw<FReal,stdComplex<FReal>,dim> Dft(N);// fast version
115 116 117 118
    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    // Initialize manually
    //  for(unsigned int s=0; s<N; ++s){
119
    //    FX[s] = FY[s] = FK[s] =stdComplex<FReal>(0.0,0.0); // init
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
    //    iFX[s]=0.0;
    //  }
    // ... or using Blas routines
    FBlas::c_setzero(N,reinterpret_cast<FReal*>(FX));
    FBlas::c_setzero(N,reinterpret_cast<FReal*>(FY));
    FBlas::c_setzero(N,reinterpret_cast<FReal*>(FK));
    FBlas::setzero(N,iFX);


    // Transform Y in Fourier space
    std::cout<< "Transform Y->FY: ";
    time.tic();
    Dft.applyDFT(Y,FY);
    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    std::cout<< "FY: "<<std::endl;
    for(unsigned int i=0; i<N; ++i){
        std::cout<< FY[i] << ", ";
    }std::cout<<std::endl;

    // Transform first column of K
    FReal tK[N];
    for(unsigned int i=0; i<N; ++i) tK[i]=K[i*N]; // first column
    //  for(unsigned int i=0; i<N; ++i) tK[i]=K[i]; // first row
    std::cout<< "Transform tK->FK: ";
    time.tic();
    Dft.applyDFT(tK,FK);
    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    //  std::cout<< "Transformed Matrix FK: "<<std::endl;
    //  for(unsigned int i=0; i<N; ++i){
    //    std::cout<< FK[i] << ", ";
    //  }std::cout<<std::endl;

    // Compute X=KY in Fourier space
    std::cout<< "Apply K in Fourier space (entrywise prod.): ";
    time.tic();
    for(unsigned int s=0; s<N; ++s){// TODO loop only over non zero entries
        FX[s]=FK[s];
        FX[s]*=FY[s];
160
    }
161 162 163 164 165 166 167 168 169 170
    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    std::cout<< "Transformed Local Exp FX=FK:FY: "<<std::endl;
    for(unsigned int i=0; i<N; ++i){
        std::cout<< FX[i] << ", ";
    }std::cout<<std::endl;

    // Transform FX back in Physical space
    std::cout<< "Transform FX->iF(FX): ";
    time.tic();
171
    Dft.applyIDFTNorm(FX,iFX);
172 173 174 175 176 177 178
    std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;

    std::cout<< "iF(FK:FY): "<<std::endl;
    for(unsigned int i=0; i<N; ++i){
        std::cout<< iFX[i] << ", ";
    }std::cout<<std::endl;

179
    std::cout << "Relative error   = " << FMath::FAccurater<FReal>( X, iFX, N) << std::endl;
180 181 182


    return 0;
183
}