FourierTransform.h

Go to the documentation of this file.

/* Copyright 2016 Kristofer Björnson
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
 
#ifndef COM_DAFER45_TBTK_FOURIER_TRANSFORM
#define COM_DAFER45_TBTK_FOURIER_TRANSFORM
 
#include "TBTK/CArray.h"
#include "TBTK/Index.h"
 
#include <fftw3.h>
 
#include <complex>
#include <vector>
 
namespace TBTK{
 
class FourierTransform{
public:
    template<typename DataType>
    class Plan{
    public:
        Plan(
            const CArray<DataType> &in,
            CArray<DataType> &out,
            const std::vector<unsigned int> &ranges,
            int sign
        );
 
        Plan(const Plan &plan) = delete;
 
        Plan(Plan &&plan);
 
        ~Plan();
 
        Plan& operator=(const Plan &plan) = delete;
 
        Plan& operator=(Plan &&plan);
 
        void setNormalizationFactor(double normalizationFactor);
 
        double getNormalizationFactor() const;
    private:
        fftw_plan *plan;
 
        double normalizationFactor;
 
        unsigned int size;
 
        const CArray<DataType> &input;
 
        CArray<DataType> &output;
 
        fftw_plan& getFFTWPlan();
 
        unsigned int getSize() const;
 
        CArray<DataType>& getInput();
 
        CArray<DataType>& getOutput();
 
        friend class FourierTransform;
    };
 
    template<typename DataType>
    class ForwardPlan : public Plan<DataType>{
    public:
        ForwardPlan(
            const CArray<DataType> &in,
            CArray<DataType> &out,
            const std::vector<unsigned int> &ranges
        ) : Plan<DataType>(in, out, ranges, -1){}
    };
 
    template<typename DataType>
    class InversePlan : public Plan<DataType>{
    public:
        InversePlan(
            const CArray<DataType> &in,
            CArray<DataType> &out,
            const std::vector<unsigned int> &ranges
        ) : Plan<DataType>(in, out, ranges, 1
        ){}
    };
 
    static void transform(
        const CArray<std::complex<double>> &in,
        CArray<std::complex<double>> &out,
        const std::vector<unsigned int> &ranges,
        int sign
    );
 
    template<typename DataType>
    static void transform(Plan<DataType> &plan);
 
    static void forward(
        const CArray<std::complex<double>> &in,
        CArray<std::complex<double>> &out,
        const std::vector<unsigned int> &ranges
    );
 
    static void inverse(
        const CArray<std::complex<double>> &in,
        CArray<std::complex<double>> &out,
        const std::vector<unsigned int> &ranges
    );
 
private:
};
 
template<typename DataType>
inline void FourierTransform::transform(Plan<DataType> &plan){
    fftw_execute(plan.getFFTWPlan());
 
    double normalizationFactor = plan.getNormalizationFactor();
    if(normalizationFactor != 1.){
        CArray<DataType> &output = plan.getOutput();
        for(unsigned int n = 0; n < plan.getSize(); n++)
            output[n] /= normalizationFactor;
    }
}
 
inline void FourierTransform::forward(
    const CArray<std::complex<double>> &in,
    CArray<std::complex<double>> &out,
    const std::vector<unsigned int> &ranges
){
    transform(in, out, ranges, -1);
}
 
inline void FourierTransform::inverse(
    const CArray<std::complex<double>> &in,
    CArray<std::complex<double>> &out,
    const std::vector<unsigned int> &ranges
){
    transform(in, out, ranges, 1);
}
 
template<typename DataType>
inline FourierTransform::Plan<DataType>::Plan(Plan &&plan){
    this->plan = plan.plan;
    plan.plan = nullptr;
 
    normalizationFactor = plan.normalizationFactor;
    size = plan.size;
    input = plan.input;
    output = plan.output;
}
 
template<typename DataType>
inline FourierTransform::Plan<DataType>::~Plan(){
    if(plan != nullptr){
        #pragma omp critical (TBTK_FOURIER_TRANSFORM)
        fftw_destroy_plan(*plan);
 
        delete plan;
    }
 
}
 
template<typename DataType>
inline FourierTransform::Plan<DataType>& FourierTransform::Plan<
    DataType
>::operator=(Plan &&rhs){
    if(this != &rhs){
        if(this->plan != nullptr){
            #pragma omp critical (TBTK_FOURIER_TRANSFORM)
            fftw_destroy_plan(*this->plan);
 
            delete this->plan;
 
            this->plan = rhs.plan;
 
            normalizationFactor = rhs.normalizationFactor;
            size = rhs.size;
            input = rhs.input;
            output = rhs.output;
        }
    }
 
    return *this;
}
 
template<typename DataType>
inline void FourierTransform::Plan<DataType>::setNormalizationFactor(
    double normalizationFactor
){
    this->normalizationFactor = normalizationFactor;
}
 
template<typename DataType>
inline double FourierTransform::Plan<DataType>::getNormalizationFactor() const{
    return normalizationFactor;
}
 
template<typename DataType>
inline fftw_plan& FourierTransform::Plan<DataType>::getFFTWPlan(){
    return *plan;
}
 
template<typename DataType>
inline unsigned int FourierTransform::Plan<DataType>::getSize() const{
    return size;
}
 
template<typename DataType>
inline CArray<DataType>& FourierTransform::Plan<DataType>::getInput(){
    return input;
}
 
template<typename DataType>
inline CArray<DataType>& FourierTransform::Plan<DataType>::getOutput(){
    return output;
}
 
};  //End of namespace TBTK
 
#endif