Homec4science
Diffusion µSpectre prototype implementation (add_eigenstrain) Browse

muSpectre_proto/src/system/fftw_engine_c2c.cc
54cf92443a69add_eigenstrain

Last Change
Raw File

fftw_engine_c2c.cc
View Options

/**
* file fftw_engine_c2c.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 11 May 2017
*
* @brief fft_engine usinge FFTW
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <complex>
#include <array>
#include "system/fftw_engine_c2c.hh"
#include <unsupported/Eigen/CXX11/Tensor>
#include <memory>
namespace muSpectre {
template<Dim_t DimS, Dim_t DimM>
FFTW_EngineC2C<DimS, DimM>::FFTW_EngineC2C(Index nb_pixels)
:FFT_Engine<DimS, DimM>(nb_pixels)
{
this->i_work_space =
fftw_alloc_complex(this->tot_nb_pixels* DimM*DimM);
this->o_work_space =
fftw_alloc_complex(this->tot_nb_pixels* DimM*DimM);
}
//----------------------------------------------------------------------------//
template<Dim_t DimS, Dim_t DimM>
FFTW_EngineC2C<DimS, DimM>::~FFTW_EngineC2C()
{
fftw_destroy_plan(this->plan_fft);
fftw_destroy_plan(this->plan_ifft);
fftw_free(this->i_work_space);
fftw_free(this->o_work_space);
fftw_cleanup();
}
//----------------------------------------------------------------------------//
template<Dim_t DimS, Dim_t DimM>
void FFTW_EngineC2C<DimS, DimM>::init(FFT_PlanFlags inflags) {
const int & rank = DimS;
std::array<int, DimS> narr;
const int * const n = &narr[0];
std::copy(this->nb_pixels.begin(), this->nb_pixels.end(), narr.begin());
int howmany = DimM*DimM;
fftw_complex * in = reinterpret_cast<fftw_complex*>(this->i_work_space);
const int * const inembed = n;
int istride = howmany;
int idist = 1;
fftw_complex * out = in;
const int * const onembed = n;
int ostride = howmany;
int odist = idist;
int sign = FFTW_FORWARD;
unsigned int flags;
if (inflags == FFT_PlanFlags::estimate) {
flags = FFTW_ESTIMATE;
}
if (inflags == FFT_PlanFlags::measure) {
flags = FFTW_MEASURE;
}
if (inflags == FFT_PlanFlags::patient) {
flags = FFTW_PATIENT;
}
this->plan_fft = fftw_plan_many_dft(rank, n, howmany, in, inembed, istride,
idist, out, onembed, ostride, odist, sign,
flags);
sign = FFTW_BACKWARD;
in = reinterpret_cast<fftw_complex*>(this->o_work_space);
out = reinterpret_cast<fftw_complex*>(this->o_work_space);
this->plan_ifft = fftw_plan_many_dft(rank, n, howmany, in, inembed, istride,
idist, out, onembed, ostride, odist, sign,
flags);
}
//----------------------------------------------------------------------------//
template<Dim_t DimS, Dim_t DimM>
Real * FFTW_EngineC2C<DimS, DimM>::convolve(const Real * const Ghat,
const Real * const arr) {
std::copy(arr, arr+this->tot_nb_pixels*DimM*DimM,
reinterpret_cast<Complex*>(this->i_work_space));
fftw_execute(this->plan_fft);
//----------------------------------------------------------------------------//
const std::array<Eigen::IndexPair<int>,2>
prod_dims {Eigen::IndexPair<int>{2, 1},
Eigen::IndexPair<int>{3, 0}};
//set output to zero
Eigen::Map<Eigen::Matrix<Complex,Eigen::Dynamic, 1>>
o_space_handle(reinterpret_cast<Complex*>(this->o_work_space),
this->tot_nb_pixels* DimM*DimM) ;
o_space_handle.setZero();
for (size_t pix_id = 0; pix_id < this->tot_nb_pixels; ++pix_id) {
using T4shape = Eigen::Sizes<DimM, DimM,
DimM,DimM>;
using T2shape = Eigen::Sizes<DimM, DimM>;
using T4type = Eigen::TensorFixedSize<Real, T4shape>;
using T2type = Eigen::TensorFixedSize<Complex, T2shape>;
using T4map = Eigen::TensorMap<T4type>;
using T2map = Eigen::TensorMap<T2type>;
const Real* const G_ptr = &Ghat[pix_id*DimM*DimM*DimM*DimM];
const T4map G(const_cast<Real*>(G_ptr),DimM, DimM, DimM,DimM);
Complex * P_ptr =
reinterpret_cast<Complex*>(&this->i_work_space[pix_id*DimM*DimM]);
const T2map P(P_ptr, DimM, DimM);
Complex* Pout_ptr =
reinterpret_cast<Complex*>(&this->o_work_space[pix_id*DimM*DimM]);
T2map Pout(Pout_ptr, DimM, DimM);
////TODO manual contraction to work around complex-
//for (Dim_t i = 0; i < DimM; ++i) {
// for (Dim_t j = 0; j < DimM; ++j) {
// for (Dim_t k = 0; k < DimM; ++k) {
// for (Dim_t l = 0; l < DimM; ++l) {
// Pout(i,j) = G(i,j,k,l)*P(l,k);
// }
// }
// }
//}
Eigen::TensorFixedSize<Complex, T4shape> Gc = G.template cast<Complex>();
Pout = Gc.contract(P, prod_dims);
}
//----------------------------------------------------------------------------//
fftw_execute(this->plan_ifft);
return reinterpret_cast<Real*>(this->o_work_space);
}
template class FFTW_EngineC2C<1, 1>;
template class FFTW_EngineC2C<2, 2>;
template class FFTW_EngineC2C<3, 3>;
} // muSpectre
c4science · Help