Commit de1d3704 authored by COULAUD Olivier's avatar COULAUD Olivier

Move complex in interpolation kernels to std::complex. Add imag() and real()...

Move complex in interpolation kernels to std::complex. Add imag() and real() methods for FComplex (compatibility with std::complex
parent 7890ce20
......@@ -4,12 +4,12 @@
#ifndef FUNIFCELL_HPP
#define FUNIFCELL_HPP
#include "Utils/stdComplex.hpp"
#include "./FUnifTensor.hpp"
#include "../../Components/FBasicCell.hpp"
#include "../../Extensions/FExtendCellType.hpp"
#include "../../Utils/FComplex.hpp"
/**
* @author Pierre Blanchard (pierre.blanchard@inria.fr)
......@@ -35,15 +35,15 @@ public:
struct exp_impl {
FReal exp[N * NVALS * VectorSize]; //< Multipole expansion
// Multipole expansion in Fourier space
FComplex<FReal> transformed_exp[N * NVALS * TransformedVectorSize];
stdComplex<FReal> transformed_exp[N * NVALS * TransformedVectorSize];
const FReal* get(const int inRhs) const
{ return this->exp + inRhs*VectorSize; }
FReal* get(const int inRhs)
{ return this->exp + inRhs*VectorSize; }
const FComplex<FReal>* getTransformed(const int inRhs) const
const stdComplex<FReal>* getTransformed(const int inRhs) const
{ return this->transformed_exp + inRhs*TransformedVectorSize; }
FComplex<FReal>* getTransformed(const int inRhs)
stdComplex<FReal>* getTransformed(const int inRhs)
{ return this->transformed_exp + inRhs*TransformedVectorSize; }
constexpr int getVectorSize() const {
......@@ -80,7 +80,7 @@ public:
FSize getSavedSize() const {
return N * NVALS * VectorSize * (FSize) sizeof(FReal)
+ N * NVALS * TransformedVectorSize * (FSize) sizeof(FComplex<FReal>);
+ N * NVALS * TransformedVectorSize * (FSize) sizeof(stdComplex<FReal>);
}
};
......
......@@ -131,7 +131,7 @@ public:
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FComplex<FReal> *const TransformedLocalExpansion = TargetExpansion->getTransformed(idxRhs);
stdComplex<FReal> *const TransformedLocalExpansion = TargetExpansion->getTransformed(idxRhs);
for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
const int idxNeigh = neighborPositions[idxExistingNeigh];
......
......@@ -16,7 +16,7 @@
#include "Utils/FDft.hpp"
#include "Utils/FSmartPointer.hpp"
#include "Utils/FComplex.hpp"
#include "Utils/stdComplex.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "FUnifTensor.hpp"
......@@ -24,7 +24,7 @@
/*! Precomputation of the 316 interactions by evaluation of the matrix kernel on the uniform grid and transformation into Fourier space.
PB: Compute() does not belong to the M2LHandler like it does in the Chebyshev kernel. This allows much nicer specialization of the M2LHandler class with respect to the homogeneity of the kernel of interaction like in the ChebyshevSym kernel.*/
template < class FReal,int ORDER, typename MatrixKernelClass>
static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal CellWidth, FComplex<FReal>* &FC, const int SeparationCriterion = 1)
static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal CellWidth, stdComplex<FReal>* &FC, const int SeparationCriterion = 1)
{
// allocate memory and store compressed M2L operators
if (FC) throw std::runtime_error("M2L operators are already set");
......@@ -41,19 +41,19 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
// allocate memory and compute 316 m2l operators (342 if separation equals 0, 343 if separation equals -1)
FReal *_C;
FComplex<FReal> *_FC;
stdComplex<FReal> *_FC;
// reduce storage from nnodes^2=order^6 to (2order-1)^3
const unsigned int rc = (2*order-1)*(2*order-1)*(2*order-1);
_C = new FReal [rc];
_FC = new FComplex<FReal> [rc * ninteractions];
_FC = new stdComplex<FReal> [rc * ninteractions];
// initialize root node ids pairs
unsigned int node_ids_pairs[rc][2];
TensorType::setNodeIdsPairs(node_ids_pairs);
// init Discrete Fourier Transformator
const int dimfft = 1; // unidim FFT since fully circulant embedding
FFftw<FReal,FComplex<FReal>,dimfft> Dft(rc);
FFftw<FReal,stdComplex<FReal>,dimfft> Dft(rc);
// get first column of K via permutation
unsigned int perm[rc];
TensorType::setStoragePermutation(perm);
......@@ -104,7 +104,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
// reduce storage if real valued kernel
const unsigned int opt_rc = rc/2+1;
// allocate M2L
FC = new FComplex<FReal>[343 * opt_rc];
FC = new stdComplex<FReal>[343 * opt_rc];
for (int i=-3; i<=3; ++i)
for (int j=-3; j<=3; ++j)
......@@ -159,10 +159,10 @@ public:
// DFT specific
static const int dimfft = 1; // unidim FFT since fully circulant embedding
typedef FFftw<FReal,FComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
typedef FFftw<FReal,stdComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
private:
/// M2L Operators (stored in Fourier space)
FSmartPointer< FComplex<FReal>,FSmartArrayMemory> FC;
FSmartPointer< stdComplex<FReal>,FSmartArrayMemory> FC;
/// Utils
unsigned int node_diff[nnodes*nnodes];
......@@ -223,12 +223,12 @@ public:
if (FC) throw std::runtime_error("M2L operator already set");
// Compute matrix of interactions
const FReal ReferenceCellWidth = FReal(2.);
FComplex<FReal>* pFC = NULL;
stdComplex<FReal>* pFC = NULL;
Compute<FReal,order>(MatrixKernel,ReferenceCellWidth,pFC,LeafLevelSeparationCriterion);
FC.assign(pFC);
// Compute memory usage
unsigned long sizeM2L = 343*opt_rc*sizeof(FComplex<FReal>);
unsigned long sizeM2L = 343*opt_rc*sizeof(stdComplex<FReal>);
// write info
......@@ -237,7 +237,7 @@ public:
}
unsigned long long getMemory() const {
return 343*opt_rc*sizeof(FComplex<FReal>);
return 343*opt_rc*sizeof(stdComplex<FReal>);
}
/**
......@@ -247,7 +247,7 @@ public:
* @param[in] X transformed local expansion of size \f$r\f$
* @param[out] x local expansion of size \f$\ell^3\f$
*/
void unapplyZeroPaddingAndDFT(const FComplex<FReal> *const FX, FReal *const x) const
void unapplyZeroPaddingAndDFT(const stdComplex<FReal> *const FX, FReal *const x) const
{
FReal Px[rc];
FBlas::setzero(rc,Px);
......@@ -273,13 +273,16 @@ public:
* @param[in] CellWidth needed for the scaling of the compressed M2L operators which are based on a homogeneous matrix kernel computed for the reference cell width \f$w=2\f$, ie in \f$[-1,1]^3\f$.
*/
void applyFC(const unsigned int idx, const unsigned int, const FReal scale,
const FComplex<FReal> *const FY, FComplex<FReal> *const FX) const
const stdComplex<FReal> *const FY, stdComplex<FReal> *const FX) const
{
// Perform entrywise product manually
for (unsigned int j=0; j<opt_rc; ++j){
FX[j].addMul(FComplex<FReal>(scale*FC[idx*opt_rc + j].getReal(),
scale*FC[idx*opt_rc + j].getImag()),
FY[j]);
// FX[j].addMul(newFComplex<FReal>(scale*FC[idx*opt_rc + j].real(),
// scale*FC[idx*opt_rc + j].imag()),
// FY[j]);
FX[j] += stdComplex<FReal>(scale*FC[idx*opt_rc + j].real(),
scale*FC[idx*opt_rc + j].imag())
*FY[j];
}
}
......@@ -291,7 +294,7 @@ public:
* @param[in] y multipole expansion of size \f$\ell^3\f$
* @param[out] Y transformed multipole expansion of size \f$r\f$
*/
void applyZeroPaddingAndDFT(FReal *const y, FComplex<FReal> *const FY) const
void applyZeroPaddingAndDFT(FReal *const y, stdComplex<FReal> *const FY) const
{
FReal Py[rc];
FBlas::setzero(rc,Py);
......@@ -303,7 +306,7 @@ public:
}
const FComplex<FReal>& getFc(const int i, const int j) const{
const stdComplex<FReal>& getFc(const int i, const int j) const{
return FC[i*opt_rc + j];
}
};
......@@ -319,7 +322,7 @@ class FUnifM2LHandler<FReal,ORDER,NON_HOMOGENEOUS>
rc = (2*ORDER-1)*(2*ORDER-1)*(2*ORDER-1)};
/// M2L Operators (stored in Fourier space for each level)
FSmartPointer< FComplex<FReal>*,FSmartArrayMemory> FC;
FSmartPointer< stdComplex<FReal>*,FSmartArrayMemory> FC;
/// Homogeneity specific variables
const unsigned int TreeHeight;
const FReal RootCellWidth;
......@@ -328,7 +331,7 @@ class FUnifM2LHandler<FReal,ORDER,NON_HOMOGENEOUS>
unsigned int node_diff[nnodes*nnodes];
/// DFT specific
static const int dimfft = 1; // unidim FFT since fully circulant embedding
typedef FFftw<FReal,FComplex<FReal>,dimfft> DftClass; // Fast real-valued Discrete Fourier Transformator
typedef FFftw<FReal,stdComplex<FReal>,dimfft> DftClass; // Fast real-valued Discrete Fourier Transformator
DftClass Dft;
const unsigned int opt_rc; // specific to real valued kernel
......@@ -357,7 +360,7 @@ public:
TensorType::setNodeIdsDiff(node_diff);
// init M2L operators
FC = new FComplex<FReal>*[TreeHeight];
FC = new stdComplex<FReal>*[TreeHeight];
FC[0] = NULL; FC[1] = NULL;
for (unsigned int l=2; l<TreeHeight; ++l) FC[l] = NULL;
......@@ -412,7 +415,7 @@ public:
}
// Compute memory usage
unsigned long sizeM2L = (TreeHeight-2)*343*opt_rc*sizeof(FComplex<FReal>);
unsigned long sizeM2L = (TreeHeight-2)*343*opt_rc*sizeof(stdComplex<FReal>);
// write info
std::cout << "Compute and set M2L operators ("<< long(sizeM2L/**1e-6*/) <<" B) in "
......@@ -420,7 +423,7 @@ public:
}
unsigned long long getMemory() const {
return (TreeHeight-2)*343*opt_rc*sizeof(FComplex<FReal>);
return (TreeHeight-2)*343*opt_rc*sizeof(stdComplex<FReal>);
}
/**
......@@ -430,7 +433,7 @@ public:
* @param[in] X transformed local expansion of size \f$r\f$
* @param[out] x local expansion of size \f$\ell^3\f$
*/
void unapplyZeroPaddingAndDFT(const FComplex<FReal> *const FX, FReal *const x) const
void unapplyZeroPaddingAndDFT(const stdComplex<FReal> *const FX, FReal *const x) const
{
FReal Px[rc];
FBlas::setzero(rc,Px);
......@@ -456,11 +459,12 @@ public:
* @param[in] CellWidth needed for the scaling of the compressed M2L operators which are based on a homogeneous matrix kernel computed for the reference cell width \f$w=2\f$, ie in \f$[-1,1]^3\f$.
*/
void applyFC(const unsigned int idx, const unsigned int TreeLevel, const FReal,
const FComplex<FReal> *const FY, FComplex<FReal> *const FX) const
const stdComplex<FReal> *const FY, stdComplex<FReal> *const FX) const
{
// Perform entrywise product manually
for (unsigned int j=0; j<opt_rc; ++j){
FX[j].addMul(FC[TreeLevel][idx*opt_rc + j],FY[j]);
// FX[j].addMul(FC[TreeLevel][idx*opt_rc + j],FY[j]);
FX[j] += (FC[TreeLevel][idx*opt_rc + j]*FY[j]);
}
}
......@@ -472,7 +476,7 @@ public:
* @param[in] y multipole expansion of size \f$\ell^3\f$
* @param[out] Y transformed multipole expansion of size \f$r\f$
*/
void applyZeroPaddingAndDFT(FReal *const y, FComplex<FReal> *const FY) const
void applyZeroPaddingAndDFT(FReal *const y, stdComplex<FReal> *const FY) const
{
FReal Py[rc];
FBlas::setzero(rc,Py);
......@@ -483,7 +487,7 @@ public:
Dft.applyDFT(Py,FY);
}
const FComplex<FReal>& getFc(const int i, const int j) const{
const stdComplex<FReal>& getFc(const int i, const int j) const{
return FC[i*opt_rc + j];
}
......
......@@ -172,7 +172,7 @@ public:
const int idxLoc = idxV*nLhs + idxLhs;
// load transformed local expansion
FComplex<FReal> *const TransformedLocalExpansion = TargetExpansion->getTransformed(idxLoc);
stdComplex<FReal> *const TransformedLocalExpansion = TargetExpansion->getTransformed(idxLoc);
// update idxRhs
const int idxRhs = idxLhs % nPV;
......
......@@ -15,7 +15,7 @@
#include "Utils/FTic.hpp"
#include "Utils/FDft.hpp"
#include "Utils/FComplex.hpp"
#include "Utils/stdComplex.hpp"
#include "FUnifTensor.hpp"
......@@ -27,7 +27,7 @@ template < class FReal, int ORDER, class MatrixKernelClass>
static void Compute(const MatrixKernelClass *const MatrixKernel,
const FReal CellWidth,
const FReal CellWidthExtension,
FComplex<FReal>** &FC,
stdComplex<FReal>** &FC,
const int SeparationCriterion = 1)
{
// dimensions of operators
......@@ -51,16 +51,16 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
// allocate memory and compute 316 m2l operators
FReal** _C;
FComplex<FReal>** _FC;
stdComplex<FReal>** _FC;
// reduce storage from nnodes^2=order^6 to (2order-1)^3
const unsigned int rc = (2*order-1)*(2*order-1)*(2*order-1);
_C = new FReal* [ncmp];
_FC = new FComplex<FReal>* [ncmp];
_FC = new stdComplex<FReal>* [ncmp];
for (unsigned int d=0; d<ncmp; ++d)
_C[d] = new FReal[rc];
for (unsigned int d=0; d<ncmp; ++d)
_FC[d] = new FComplex<FReal>[rc * ninteractions];
_FC[d] = new stdComplex<FReal>[rc * ninteractions];
// initialize root node ids pairs
unsigned int node_ids_pairs[rc][2];
......@@ -68,7 +68,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
// init Discrete Fourier Transformator
const int dimfft = 1; // unidim FFT since fully circulant embedding
FFftw<FReal,FComplex<FReal>,dimfft> Dft(rc);
FFftw<FReal,stdComplex<FReal>,dimfft> Dft(rc);
// get first column of K via permutation
unsigned int perm[rc];
......@@ -110,7 +110,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
Dft.applyDFT(_C[d],_FC[d]+counter*rc);
// increment interaction counter
counter++;
++counter;
}
}
}
......@@ -131,7 +131,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
const unsigned int opt_rc = rc/2+1;
// allocate M2L
for (unsigned int d=0; d<ncmp; ++d)
FC[d] = new FComplex<FReal>[343 * opt_rc];
FC[d] = new stdComplex<FReal>[343 * opt_rc];
for (int i=-3; i<=3; ++i)
for (int j=-3; j<=3; ++j)
......@@ -141,7 +141,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
for (unsigned int d=0; d<ncmp; ++d)
FBlas::c_copy(opt_rc, reinterpret_cast<FReal*>(_FC[d] + counter*rc),
reinterpret_cast<FReal*>(FC[d] + idx*opt_rc));
counter++;
++counter;
} else {
for (unsigned int d=0; d<ncmp; ++d)
FBlas::c_setzero(opt_rc, reinterpret_cast<FReal*>(FC[d] + idx*opt_rc));
......@@ -181,7 +181,7 @@ class FUnifTensorialM2LHandler<FReal, ORDER,MatrixKernelClass,HOMOGENEOUS>
ncmp = MatrixKernelClass::NCMP};
/// M2L Operators (stored in Fourier space for each component of tensor)
FSmartPointer< FComplex<FReal>*,FSmartArrayMemory> FC;
FSmartPointer< stdComplex<FReal>*,FSmartArrayMemory> FC;
const FReal CellWidthExtension; //<! extension of cells width
......@@ -191,7 +191,7 @@ class FUnifTensorialM2LHandler<FReal, ORDER,MatrixKernelClass,HOMOGENEOUS>
/// DFT specific
static const int dimfft = 1; // unidim FFT since fully circulant embedding
typedef FFftw<FReal,FComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
typedef FFftw<FReal,stdComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
DftClass Dft;
const unsigned int opt_rc; // specific to real valued kernel
......@@ -215,7 +215,7 @@ public:
{
// allocate FC
FC = new FComplex<FReal>*[ncmp];
FC = new stdComplex<FReal>*[ncmp];
for (unsigned int d=0; d<ncmp; ++d)
FC[d]=nullptr;
......@@ -264,7 +264,7 @@ public:
Compute<FReal,order>(MatrixKernel,ReferenceCellWidth, 0., FC, LeafLevelSeparationCriterion);
// Compute memory usage
unsigned long sizeM2L = 343*ncmp*opt_rc*sizeof(FComplex<FReal>);
unsigned long sizeM2L = 343*ncmp*opt_rc*sizeof(stdComplex<FReal>);
// write info
std::cout << "Compute and set M2L operators ("<< long(sizeM2L/**1e-6*/) <<" B) in "
......@@ -272,7 +272,7 @@ public:
}
unsigned long long getMemory() const {
return 343*ncmp*opt_rc*sizeof(FComplex<FReal>);
return 343*ncmp*opt_rc*sizeof(stdComplex<FReal>);
}
/**
......@@ -282,7 +282,7 @@ public:
* @param[in] X transformed local expansion of size \f$r\f$
* @param[out] x local expansion of size \f$\ell^3\f$
*/
void unapplyZeroPaddingAndDFT(const FComplex<FReal> *const FX, FReal *const x) const
void unapplyZeroPaddingAndDFT(const stdComplex<FReal> *const FX, FReal *const x) const
{
FReal Px[rc];
FBlas::setzero(rc,Px);
......@@ -309,11 +309,11 @@ public:
*/
void applyFC(const unsigned int idx, const unsigned int ,
const FReal scale, const unsigned int d,
const FComplex<FReal> *const FY, FComplex<FReal> *const FX) const
const stdComplex<FReal> *const FY, stdComplex<FReal> *const FX) const
{
// Perform entrywise product manually
for (unsigned int j=0; j<opt_rc; ++j){
FX[j].addMul(FComplex<FReal>(scale*FC[d][idx*opt_rc + j].getReal(),
FX[j].addMul(stdComplex<FReal>(scale*FC[d][idx*opt_rc + j].getReal(),
scale*FC[d][idx*opt_rc + j].getImag()),
FY[j]);
}
......@@ -328,7 +328,7 @@ public:
* @param[in] y multipole expansion of size \f$\ell^3\f$
* @param[out] Y transformed multipole expansion of size \f$r\f$
*/
void applyZeroPaddingAndDFT(FReal *const y, FComplex<FReal> *const FY) const
void applyZeroPaddingAndDFT(FReal *const y, stdComplex<FReal> *const FY) const
{
FReal Py[rc];
FBlas::setzero(rc,Py);
......@@ -355,7 +355,7 @@ class FUnifTensorialM2LHandler<FReal,ORDER,MatrixKernelClass,NON_HOMOGENEOUS>
/// M2L Operators (stored in Fourier space for each component and each level)
FSmartPointer< FComplex<FReal>**,FSmartArrayMemory> FC;
FSmartPointer< stdComplex<FReal>**,FSmartArrayMemory> FC;
/// Homogeneity specific variables
const unsigned int TreeHeight;
......@@ -368,7 +368,7 @@ class FUnifTensorialM2LHandler<FReal,ORDER,MatrixKernelClass,NON_HOMOGENEOUS>
/// DFT specific
static const int dimfft = 1; // unidim FFT since fully circulant embedding
typedef FFftw<FReal,FComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
typedef FFftw<FReal,stdComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
DftClass Dft;
const unsigned int opt_rc; // specific to real valued kernel
......@@ -394,9 +394,9 @@ public:
{
// allocate FC
FC = new FComplex<FReal>**[TreeHeight];
FC = new stdComplex<FReal>**[TreeHeight];
for (unsigned int l=0; l<TreeHeight; ++l){
FC[l] = new FComplex<FReal>*[ncmp];
FC[l] = new stdComplex<FReal>*[ncmp];
for (unsigned int d=0; d<ncmp; ++d)
FC[l][d]=nullptr;
}
......@@ -453,7 +453,7 @@ public:
}
// Compute memory usage
unsigned long sizeM2L = (TreeHeight-2)*343*ncmp*opt_rc*sizeof(FComplex<FReal>);
unsigned long sizeM2L = (TreeHeight-2)*343*ncmp*opt_rc*sizeof(stdComplex<FReal>);
// write info
std::cout << "Compute and set M2L operators ("<< long(sizeM2L/**1e-6*/) <<" B) in "
......@@ -461,7 +461,7 @@ public:
}
unsigned long long getMemory() const {
return (TreeHeight-2)*343*ncmp*opt_rc*sizeof(FComplex<FReal>);
return (TreeHeight-2)*343*ncmp*opt_rc*sizeof(stdComplex<FReal>);
}
/**
......@@ -471,7 +471,7 @@ public:
* @param[in] X transformed local expansion of size \f$r\f$
* @param[out] x local expansion of size \f$\ell^3\f$
*/
void unapplyZeroPaddingAndDFT(const FComplex<FReal> *const FX, FReal *const x) const
void unapplyZeroPaddingAndDFT(const stdComplex<FReal> *const FX, FReal *const x) const
{
FReal Px[rc];
FBlas::setzero(rc,Px);
......@@ -498,11 +498,12 @@ public:
*/
void applyFC(const unsigned int idx, const unsigned int TreeLevel,
const FReal, const unsigned int d,
const FComplex<FReal> *const FY, FComplex<FReal> *const FX) const
const stdComplex<FReal> *const FY, stdComplex<FReal> *const FX) const
{
// Perform entrywise product manually
for (unsigned int j=0; j<opt_rc; ++j){
FX[j].addMul(FC[TreeLevel][d][idx*opt_rc + j], FY[j]);
// FX[j].addMul(FC[TreeLevel][d][idx*opt_rc + j], FY[j]);
FX[j] += (FC[TreeLevel][d][idx*opt_rc + j]*FY[j]);
}
}
......@@ -515,7 +516,7 @@ public:
* @param[in] y multipole expansion of size \f$\ell^3\f$
* @param[out] Y transformed multipole expansion of size \f$r\f$
*/
void applyZeroPaddingAndDFT(FReal *const y, FComplex<FReal> *const FY) const
void applyZeroPaddingAndDFT(FReal *const y, stdComplex<FReal> *const FY) const
{
FReal Py[rc];
FBlas::setzero(rc,Py);
......
......@@ -74,11 +74,19 @@ public:
return !(*this == other);
}
/** Get imaginary */
FReal imag() const{
return this->complex[1];
}
/** Get imaginary */
FReal getImag() const{
return this->complex[1];
}
/** Get complex[0] */
FReal real() const{
return this->complex[0];
}
/** Get complex[0] */
FReal getReal() const{
return this->complex[0];
......
This diff is collapsed.
// ===================================================================================
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program 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 and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef STDCOMPLEXE_HPP
#define STDCOMPLEXE_HPP
#include <complex>
template<typename T>
using stdComplex = std::complex<T> ;
#endif //STDCOMPLEXE_HPP
......@@ -28,7 +28,8 @@
#endif
#include "Utils/FGlobal.hpp"
#include "Utils/FComplex.hpp"
#include "Utils/FMath.hpp"
#include "Utils/stdComplex.hpp"
#include "Utils/FTic.hpp"
......@@ -59,9 +60,9 @@ int main(int argc, char** argv)
// input/output types
typedef FReal FftR2CInputType;
typedef FComplex<FReal> FftR2COutputType;
typedef FComplex<FReal> FftC2CInputType;
typedef FComplex<FReal> FftC2COutputType;
typedef stdComplex<FReal> FftR2COutputType;
typedef stdComplex<FReal> FftC2CInputType;
typedef stdComplex<FReal> FftC2COutputType;
// fftw arrays
FftR2CInputType* FftR2CInput = new FftR2CInputType[total_size];
......@@ -74,8 +75,8 @@ int main(int argc, char** argv)
// fftw wrappers
typedef FFftw<FReal ,FComplex<FReal>,dim> FftwR2CClass;
typedef FFftw<FComplex<FReal>,FComplex<FReal>,dim> FftwC2CClass;
typedef FFftw<FReal ,stdComplex<FReal>,dim> FftwR2CClass;
typedef FFftw<stdComplex<FReal>,stdComplex<FReal>,dim> FftwC2CClass;
std::cout<< "Init FFTW wrappers: ";
time.tic();
......@@ -160,8 +161,8 @@ int main(int argc, char** argv)
FReal* FftC2CInvOutReal = new FReal[total_size]; FReal* FftC2CInputReal = new FReal[total_size];
FReal* FftC2CInvOutImag = new FReal[total_size]; FReal* FftC2CInputImag = new FReal[total_size];
for( int s=0; s<total_size; ++s){
FftC2CInvOutReal[s] = FftC2CInvOut[s].getReal(); FftC2CInputReal[s] = FftC2CInput[s].getReal();
FftC2CInvOutImag[s] = FftC2CInvOut[s].getImag(); FftC2CInputImag[s] = FftC2CInput[s].getImag();
FftC2CInvOutReal[s] = FftC2CInvOut[s].real(); FftC2CInputReal[s] = FftC2CInput[s].real();
FftC2CInvOutImag[s] = FftC2CInvOut[s].imag(); FftC2CInputImag[s] = FftC2CInput[s].imag();
}
if(FMath::FAccurater<FReal>(FftR2CInput, FftR2CInvOut, total_size).getRelativeInfNorm() > threshold) isWorking=false;
......
......@@ -33,7 +33,7 @@
#include "Utils/FBlas.hpp"
//
#include "Utils/FComplex.hpp"
#include "Utils/stdComplex.hpp"
#include "Utils/FDft.hpp"
#include "Utils/FParameterNames.hpp"
......@@ -101,9 +101,9 @@ int main(int argc, char ** argv){
std::cout<<std::endl;
// now compute via DFT and use convolution theorem
FComplex<FReal> FK[N];
FComplex<FReal> FY[N];
FComplex<FReal> FX[N];
stdComplex<FReal> FK[N];
stdComplex<FReal> FY[N];
stdComplex<FReal> FX[N];
FReal iFX[N];
// Init DFTor
......@@ -111,12 +111,12 @@ int main(int argc, char ** argv){
time.tic();
const int dim = 1;
//FDft<FReal> Dft(N);// direct version (Beware! Ordering of output differs from REAL valued-FFT)
FFftw<FReal,FComplex<FReal>,dim> Dft(N);// fast version
FFftw<FReal,stdComplex<FReal>,dim> Dft(N);// fast version
std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;
// Initialize manually
// for(unsigned int s=0; s<N; ++s){
// FX[s] = FY[s] = FK[s] =FComplex<FReal>(0.0,0.0); // init
// FX[s] = FY[s] = FK[s] =stdComplex<FReal>(0.0,0.0); // init
// iFX[s]=0.0;
// }
// ... or using Blas routines
......
......@@ -42,7 +42,7 @@
#include "Kernels/Uniform/FUnifTensor.hpp"
// Check DFT
#include "Utils/FComplex.hpp"
#include "Utils/stdComplex.hpp"
#include "Utils/FDft.hpp"
......@@ -333,11 +333,11 @@ int main(int argc, char ** argv){
/////////////////////////////////////////////////////////////////////////////////////
// Efficient application of the Toeplitz system in FOURIER SPACE
FComplex<FReal> FPMultExp[rc];
FComplex<FReal> FPLocalExp[rc];
stdComplex<FReal> FPMultExp[rc];
stdComplex<FReal> FPLocalExp[rc];
FReal PLocalExp[rc];
//for (unsigned int n=0; n<rc; ++n) FPLocalExp[n]=FComplex<FReal>(0.0,0.0);
//for (unsigned int n=0; n<rc; ++n) FPLocalExp[n]=stdComplex<FReal>(0.0,0.0);
FBlas::c_setzero(rc,reinterpret_cast<FReal*>(FPLocalExp));
FBlas::setzero(rc,PLocalExp);
......@@ -345,7 +345,7 @@ int main(int argc, char ** argv){
// Init DFT
const int dimfft = 1;
//FDft Dft(rc); // direct version
FFftw<FReal,FComplex<FReal>,dimfft> Dft(rc); // fast version
FFftw<FReal,stdComplex<FReal>,dimfft> Dft(rc); // fast version
// Get first COLUMN of K and Store in T
FReal T[rc];
......@@ -372,8 +372,8 @@ int main(int argc, char ** argv){
// std::cout<<std::endl;
// Apply DFT to T
FComplex<FReal> FT[rc];
// for (unsigned int n=0; n<rc; ++n) FT[n]=FComplex<FReal>(0.0,0.0);
stdComplex<FReal> FT[rc];
// for (unsigned int n=0; n<rc; ++n) FT[n]=stdComplex<FReal>(0.0,0.0);
FBlas::c_setzero(rc,reinterpret_cast<FReal*>(FT));
// if first COLUMN (T) of C is used
......@@ -397,7 +397,7 @@ int main(int argc, char ** argv){
// Set transformed MultExp to 0
// for (unsigned int n=0; n<rc; ++n) FPMultExp[n]=FComplex<FReal>(0.0,0.0);
// for (unsigned int n=0; n<rc; ++n) FPMultExp[n]=stdComplex<FReal>(0.0,0.0);
FBlas::c_setzero(rc,reinterpret_cast<FReal*>(FPMultExp));
// Transform PaddedMultExp
......
......@@ -42,7 +42,7 @@
#include "Kernels/Uniform/FUnifTensor.hpp"
// Check DFT
#include "Utils/FComplex.hpp"
#include "Utils/stdComplex.hpp"
#include "Utils/FDft.hpp"
......@@ -362,7 +362,7 @@ int main(int argc, char ** argv){
// Init DFT
const int dimfft = 1;
//FDft Dft(rc); // direct version
FFftw<FReal,FComplex<FReal>,dimfft> Dft(rc); // fast version
FFftw<FReal,stdComplex<FReal>,dimfft> Dft(rc); // fast version
// Get first COLUMN of K and Store in T
FReal T[ncmp*rc];
......@@ -390,8 +390,8 @@ int main(int argc, char ** argv){
// std::cout<<std::endl;
// Apply DFT to T
FComplex<FReal> FT[ncmp*rc];
// for (unsigned int n=0; n<rc; ++n) FT[n]=FComplex<FReal>(0.0,0.0);
stdComplex<FReal> FT[ncmp*rc];
// for (unsigned int n=0; n<rc; ++n) FT[n]=stdComplex<FReal>(0.0,0.0);
FBlas::c_setzero(ncmp*rc,reinterpret_cast<FReal*>(FT));
// if first COLUMN (T) of C is used
......@@ -400,11 +400,11 @@ int main(int argc, char ** argv){
// // if first ROW of C is used
// Dft.applyDFT(C,FT);