Commit 2f1e1948 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 24/3/16 transfert code from test_suite to laguerre2bessel. Cleanup code....

(JEC) 24/3/16 transfert code from test_suite to laguerre2bessel. Cleanup code. remove using namespace in global space
parent fb5907a8
......@@ -140,6 +140,7 @@ CXXHDR = $(SRC)lagsht_exceptions.h \
$(SRC)walltimer.h \
$(SRC)lagsht_bessel.h \
$(SRC)laguerre2bessel.h \
$(SRC)lagsht_func.h \
$(SRC)quadinteg.h
CPPFLAGS += $(BLASYES) $(SHARPINC) $(BLASINC) $(BOOSTINC) $(FFTWINC)
......
......@@ -15,9 +15,9 @@
namespace LagSHT {
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
vector<r_8>& fijk,
vector< complex<r_8> >& flmk) {
void LaguerreSphericalTransform::Synthesis(const std::vector< std::complex<r_8> >& flmn,
std::vector<r_8>& fijk,
std::vector< std::complex<r_8> >& flmk) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
if(flmn.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
......@@ -25,7 +25,6 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
fijk.resize(Ntot3DPix);
//Perform the Inverse Laguerre Transform first
// vector< complex<r_8> > flmk(NtotCoeff);
flmk.resize(NtotCoeff);
......@@ -52,7 +51,7 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
off7k = n*Npix_;
#if DEBUG >= 2
cout << "Synthesis... at shell num: " << n << endl;
std::cout << "Synthesis... at shell num: " << n << std::endl;
#endif
sphere_->alm2map(reinterpret_cast<const r_8*>(&(flmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation
......@@ -61,15 +60,15 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
tstack_pop("SHT Synthesis");
#if DEBUG >= 1
cout << "SphLagTransform::Synthesis.... Done" << endl;
std::cout << "SphLagTransform::Synthesis.... Done" << std::endl;
#endif
}//Synthesis
void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
vector< complex<r_8> >& flmn,
vector< complex<r_8> >& almk
void LaguerreSphericalTransform::Analysis(const std::vector<r_8>& fijk,
std::vector< std::complex<r_8> >& flmn,
std::vector< std::complex<r_8> >& almk
) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
......@@ -79,7 +78,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
//done in Laguerre routine flmn.resize(NtotCoeff);
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Go..." << endl;
std::cout << "SphLagTransform::Analysis.... Go..." << std::endl;
#endif
//Perform the Spherical Hamonic analysis first
......@@ -93,7 +92,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
off7k = n*Npix_;
#if DEBUG >= 1
cout << "Start Harmonic Spheric... at n = "<< n << endl;
std::cout << "Start Harmonic Spheric... at n = "<< n << std::endl;
#endif
sphere_->map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(almk.data()[off7n])),false); //the false is to forbidd accumulation
......@@ -101,8 +100,8 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
}//loop on shell
#if DEBUG >= 1
cout << "End Harmonic Spheric..." << endl;
cout << "Start Laguerre..." << endl;
std::cout << "End Harmonic Spheric..." << std::endl;
std::cout << "Start Laguerre..." << std::endl;
#endif
tstack_pop("SHT Analysis");
......@@ -112,7 +111,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
lagTrans_->MultiAnalysis(almk, flmn, Nalm_);
#if DEBUG >= 1
cout << "End Laguerre..." << endl;
std::cout << "End Laguerre..." << std::endl;
#endif
tstack_pop("Laguerre MultiAnalysis");
......@@ -120,20 +119,22 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
#if DEBUG >= 2
for (int i=0;i<NtotCoeff; i++) {
cout << "flmn("<<i<<"): " <<flmn[i] << endl;
std::cout << "flmn("<<i<<"): " <<flmn[i] << std::endl;
}
#endif
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Done" << endl;
std::cout << "SphLagTransform::Analysis.... Done" << std::endl;
#endif
}//Analysis
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, const vector<complex<r_8> >& Blmn,
void LaguerreSphericalTransform::Synthesis(const std::vector< std::complex<r_8> >& Elmn,
const std::vector< std::complex<r_8> >& Blmn,
int spin,
vector<r_8>& fijk_re, vector<r_8>& fijk_im,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
std::vector<r_8>& fijk_re, std::vector<r_8>& fijk_im,
std::vector< std::complex<r_8> >& Elmk,
std::vector< std::complex<r_8> >& Blmk
) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
......@@ -148,8 +149,8 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, c
Blmk.resize(NtotCoeff);
#if DEBUG >= 1
cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << endl;
cout << "Start Inverse Laguerre...." << endl;
std::cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << std::endl;
std::cout << "Start Inverse Laguerre...." << std::endl;
#endif
tstack_push("Laguerre MultiSynthesis");
......@@ -157,17 +158,9 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, c
lagTrans_->MultiSynthesis(Blmn, Blmk, Nalm_);
tstack_pop("Laguerre MultiSynthesis");
// cout << "Elmk, Blmk" <<endl;
// for(int i=0; i<100; i++){
// cout << "Elmk["<<i<<"]: " << Elmk[i] << "; "
// << "Blmk["<<i<<"]: " << Blmk[i] << endl;
// }
#if DEBUG >= 1
cout << "End Inverse Laguerre...." << endl;
cout << "Start Inverse H.Spherical Transform " << endl;
std::cout << "End Inverse Laguerre...." << std::endl;
std::cout << "Start Inverse H.Spherical Transform " << std::endl;
#endif
tstack_push("SHT Synthesis");
......@@ -178,7 +171,7 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, c
off7k = n*Npix_;
#if DEBUG >= 2
cout << "Synthesis... at shell num: " << n << endl;
std::cout << "Synthesis... at shell num: " << n << std::endl;
#endif
sphere_->alm2map_spin(reinterpret_cast<const r_8*>(&(Elmk.data()[off7n])),
reinterpret_cast<const r_8*>(&(Blmk.data()[off7n])),
......@@ -192,15 +185,18 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, c
tstack_pop("SHT Synthesis");
#if DEBUG >= 1
cout << "SphLagTransform::Synthesis.... Done" << endl;
std::cout << "SphLagTransform::Synthesis.... Done" << std::endl;
#endif
}//Synthesis
void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vector<r_8>& fijk_im,
void LaguerreSphericalTransform::Analysis(const std::vector<r_8>& fijk_re,
const std::vector<r_8>& fijk_im,
int spin,
vector< complex<r_8> >& Elmn, vector<complex<r_8> >& Blmn,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
std::vector< std::complex<r_8> >& Elmn,
std::vector< std::complex<r_8> >& Blmn,
std::vector< std::complex<r_8> >& Elmk,
std::vector< std::complex<r_8> >& Blmk
) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
......@@ -215,7 +211,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vect
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Go..." << endl;
std::cout << "SphLagTransform::Analysis.... Go..." << std::endl;
#endif
//Perform the Spherical Hamonic analysis first
......@@ -229,7 +225,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vect
off7k = n*Npix_;
#if DEBUG >= 1
cout << "Start Harmonic Spheric... at n = "<< n << endl;
std::cout << "Start Harmonic Spheric... at n = "<< n << std::endl;
#endif
sphere_->map2alm_spin(&(fijk_re.data()[off7k]),
......@@ -243,8 +239,8 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vect
}//loop on shell
#if DEBUG >= 1
cout << "End Harmonic Spheric..." << endl;
cout << "Start Laguerre..." << endl;
std::cout << "End Harmonic Spheric..." << std::endl;
std::cout << "Start Laguerre..." << std::endl;
#endif
tstack_pop("SHT Analysis");
......@@ -255,18 +251,18 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vect
lagTrans_->MultiAnalysis(Blmk, Blmn, Nalm_);
#if DEBUG >= 1
cout << "End Laguerre..." << endl;
std::cout << "End Laguerre..." << std::endl;
#endif
tstack_pop("Laguerre MultiAnalysis");
#if DEBUG >= 2
for (int i=0;i<NtotCoeff; i++) {
cout << "flmn("<<i<<"): " << Elmn[i] << ", " << Blmn[i] << endl;
std::cout << "flmn("<<i<<"): " << Elmn[i] << ", " << Blmn[i] << std::endl;
}
#endif
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Done" << endl;
std::cout << "SphLagTransform::Analysis.... Done" << std::endl;
#endif
}//Analysis
......
......@@ -54,7 +54,7 @@ class LaguerreSphericalTransform {
// LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax,
// r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2, bool R2tau=true):
LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax,
LaguerreSphericalTransform(std::string spheregeom, int Lmax, int Mmax, int Nmax,
r_8 R, int Nrings = -1, int Nphi = -1, r_8 alpha=2, bool R2tau=true):
Lmax_(Lmax), Mmax_(Mmax), R_(R) {
......@@ -66,7 +66,7 @@ class LaguerreSphericalTransform {
} else if (spheregeom == "Healpix") {
sphere_ = new HealpixGeometry(Lmax, Mmax, Nrings); //here Nsides = Nrings
} else {
cout << "FAILURE: LaguerreSphericalTransform::Ctor wrong geometry type" << endl;
std::cout << "FAILURE: LaguerreSphericalTransform::Ctor wrong geometry type" << std::endl;
sphere_ = 0;
throw LagSHTError("LaguerreSphericalTransform::Ctor wrong geometry type");
}
......@@ -119,12 +119,12 @@ class LaguerreSphericalTransform {
//! Get the f(l,m,n) coefficient in the vector collection
//make a Flmn Class
inline complex<r_8> GetFlmn(const vector< complex<r_8> >& flmn, int l , int m, int n) const {
inline std::complex<r_8> GetFlmn(const std::vector< std::complex<r_8> >& flmn, int l , int m, int n) const {
int id= n*Nalm_ + l +m*Lmax_ -m*(m+1)/2;
return flmn[id];
}
//! Set the f(l,m,n) coefficient in the vector collection
inline void SetFlmn(vector< complex<r_8> >& flmn, int l , int m, int n, complex<r_8> value) const {
inline void SetFlmn(std::vector< std::complex<r_8> >& flmn, int l , int m, int n, std::complex<r_8> value) const {
int id= n*Nalm_ + l +m*Lmax_ -m*(m+1)/2;
flmn[id] = value;
}
......@@ -137,18 +137,18 @@ class LaguerreSphericalTransform {
\output fijk: 3D real function values
\output almk: 2D alm complex spherical coefficients for each subshell k
*/
void Synthesis(const vector< complex<r_8> >& flmn,
vector<r_8>& fijk,
vector< complex<r_8> >& almk);
void Synthesis(const std::vector< std::complex<r_8> >& flmn,
std::vector<r_8>& fijk,
std::vector< std::complex<r_8> >& almk);
//!Synthesis (spin = 0)
/*! \brief Coeffs -> Pixels with Input/Output using T floating representation
\input flmn: 3D complex spherical-laguerre coefficients
\output fijk: 3D real function values
*/
void Synthesis(const vector< complex<r_8> >& flmn, vector<r_8>& fijk) {
void Synthesis(const std::vector< std::complex<r_8> >& flmn, std::vector<r_8>& fijk) {
int NtotCoeff = Nalm_*N_;
vector< complex<r_8> > almk(NtotCoeff);
std::vector< std::complex<r_8> > almk(NtotCoeff);
Synthesis(flmn,fijk,almk);
}
......@@ -158,9 +158,9 @@ class LaguerreSphericalTransform {
\input fijk: 3D real function values
\output flmn: 3D complex spherical-laguerre coefficients
*/
void Analysis(const vector<r_8>& fijk,
vector< complex<r_8> >& flmn,
vector< complex<r_8> >& almk
void Analysis(const std::vector<r_8>& fijk,
std::vector< std::complex<r_8> >& flmn,
std::vector< std::complex<r_8> >& almk
);
//! Analysis (spin = 0)
......@@ -168,9 +168,9 @@ class LaguerreSphericalTransform {
\input fijk: 3D real function values
\output flmn: 3D complex spherical-laguerre coefficients
*/
void Analysis(const vector<r_8>& fijk, vector< complex<r_8> >& flmn){
void Analysis(const std::vector<r_8>& fijk, std::vector< std::complex<r_8> >& flmn){
int NtotCoeff = Nalm_*N_;
vector< complex<r_8> > almk(NtotCoeff);
std::vector< std::complex<r_8> > almk(NtotCoeff);
Analysis(fijk,flmn,almk);
}
......@@ -185,20 +185,22 @@ class LaguerreSphericalTransform {
\output Elmk: 2D alm complex spherical coefficients for each subshell k of the gradient E
\output Blmk: 2D alm complex spherical coefficients for each subshell k of the curl B
*/
void Synthesis(const vector< complex<r_8> >& Elmn, const vector<complex<r_8> >& Blmn,
void Synthesis(const std::vector< std::complex<r_8> >& Elmn,
const std::vector< std::complex<r_8> >& Blmn,
int spin,
vector<r_8>& fijk_re, vector<r_8>& fijk_im,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
std::vector<r_8>& fijk_re, std::vector<r_8>& fijk_im,
std::vector< std::complex<r_8> >& Elmk, std::vector< std::complex<r_8> >& Blmk
);
void Synthesis(const vector< complex<r_8> >& Elmn, const vector<complex<r_8> >& Blmn,
void Synthesis(const std::vector< std::complex<r_8> >& Elmn,
const std::vector< std::complex<r_8> >& Blmn,
int spin,
vector<r_8>& fijk_re, vector<r_8>& fijk_im
std::vector<r_8>& fijk_re, std::vector<r_8>& fijk_im
){
int NtotCoeff = Nalm_*N_;
vector< complex<r_8> > Elmk(NtotCoeff);
vector< complex<r_8> > Blmk(NtotCoeff);
std::vector< std::complex<r_8> > Elmk(NtotCoeff);
std::vector< std::complex<r_8> > Blmk(NtotCoeff);
Synthesis(Elmn,Blmn,spin,fijk_re,fijk_im,Elmk,Blmk);
}
......@@ -211,19 +213,19 @@ class LaguerreSphericalTransform {
\output Elmk: 2D alm complex spherical coefficients for each subshell k of the gradient E
\output Blmk: 2D alm complex spherical coefficients for each subshell k of the curl B
*/
void Analysis(const vector<r_8>& fijk_re, const vector<r_8>& fijk_im,
void Analysis(const std::vector<r_8>& fijk_re, const std::vector<r_8>& fijk_im,
int spin,
vector< complex<r_8> >& Elmn, vector<complex<r_8> >& Blmn,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
std::vector< std::complex<r_8> >& Elmn, std::vector< std::complex<r_8> >& Blmn,
std::vector< std::complex<r_8> >& Elmk, std::vector< std::complex<r_8> >& Blmk
);
void Analysis(const vector<r_8>& fijk_re, const vector<r_8>& fijk_im,
void Analysis(const std::vector<r_8>& fijk_re, const std::vector<r_8>& fijk_im,
int spin,
vector< complex<r_8> >& Elmn, vector<complex<r_8> >& Blmn
std::vector< std::complex<r_8> >& Elmn, std::vector< std::complex<r_8> >& Blmn
){
int NtotCoeff = Nalm_*N_;
vector< complex<r_8> > Elmk(NtotCoeff);
vector< complex<r_8> > Blmk(NtotCoeff);
std::vector< std::complex<r_8> > Elmk(NtotCoeff);
std::vector< std::complex<r_8> > Blmk(NtotCoeff);
Analysis(fijk_re, fijk_im, spin, Elmn, Blmn, Elmk, Blmk);
}
......
......@@ -5,7 +5,7 @@
namespace LagSHT {
void Bessel::BesselRoots() {
// using namespace std;
using namespace std;
string msg;
......
......@@ -33,7 +33,7 @@
#include "boost/math/special_functions/bessel.hpp"
#include "lagsht_numbers.h"
using namespace std;
namespace LagSHT {
......@@ -79,9 +79,10 @@ class Bessel {
/* if(l>lmax_ || n > nmax_) cout <<" ERROR Besel::Op(l,n) called with ("<<l<<","<<n<<") while " */
/* << " lmax = " << lmax_ << " nmax = " << nmax_ <<endl; */
return qln_[l*nmax_+n];} //qln(l,n)
void GetVecRoots(std::vector<r_8>& roots, int l) const {
roots.resize(nmax_);
copy(qln_.begin()+l*nmax_,qln_.begin()+(l+1)*nmax_,roots.begin());
std::copy(qln_.begin()+l*nmax_,qln_.begin()+(l+1)*nmax_,roots.begin());
}
//Find Bessel root (qln) values such that j_l(qln) = 0 for l=O,...,Lmax and n=1,..,Nmax
......@@ -93,7 +94,7 @@ class Bessel {
int lmax_; //!< l=0,...,lmax-1 order of Spherical Bessel indice
int nmax_; //!< n=0,...,nmax-1 index of the roots of Spherical Bessel func
vector<r_8> qln_; //!< roots: Care qln(l,n) l=0,...,lmax-1; n=0,...,nmax-1 but stands for n=1,..,nmax
std::vector<r_8> qln_; //!< roots: Care qln(l,n) l=0,...,lmax-1; n=0,...,nmax-1 but stands for n=1,..,nmax
private:
......
......@@ -27,6 +27,8 @@ void BaseGeometry::SetCommonMapFeatures() {
}//SetThetaPhiMap
void BaseGeometry::SetPixelPositions() {
using namespace std;
const sharp_geom_info *ginfo;
ginfo = sharpjob_.get_geom_info();
for (int i=0; i<ginfo->npairs; ++i) {
......@@ -87,6 +89,7 @@ void ECPGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
}//PixThetaPhi
void ECPGeometry::SetMap(int nrings, int nphi) {
using namespace std;
geomType_ = "ECP";
......@@ -146,6 +149,8 @@ void GaussGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
void GaussGeometry::SetMap(int nrings, int nphi) {
using namespace std;
geomType_ = "Gauss";
if(geomDone_){
......@@ -193,6 +198,7 @@ void HealpixGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
}//PixThetaPhi
void HealpixGeometry::SetMap(int nsides, int DUMMY) { //second argument not used
using namespace std;
geomType_ = "Healpix";
......
......@@ -29,7 +29,7 @@
#include <iostream>
#include <algorithm>
#include <utility> // std::pair, std::make_pair
using namespace std;
//using namespace std;
#include "lagsht_geom.h"
//Sharp
......@@ -79,7 +79,7 @@ class BaseGeometry {
//! Mmax
int Mmax() const { return M_; }
//! Type of geomery
string Geometry() const { return geomType_; }
std::string Geometry() const { return geomType_; }
//! Sorted List of pixel center (theta, phi) positions
geom_t GetPixelPositions() const { return vecPix_; }
......@@ -162,8 +162,8 @@ protected:
int Nalm_; //!< Total number of Alm coefficients
geom_t vecPix_; //!< sorted vector of pixel positions
vector<r_8> thetaPix_; //!< sorted vector of theta rings
string geomType_; //!<< Type of pixelization
std::vector<r_8> thetaPix_; //!< sorted vector of theta rings
std::string geomType_; //!<< Type of pixelization
bool geomDone_; //!< Set to true once the geometry is properly set (See SetThetaPhiMap(string) )
sharp_cxxjob<r_8> sharpjob_; //!< libsharp job // voir si r_8 -> r_16 extend the validity of the computation
......
......@@ -29,10 +29,6 @@ using namespace std;
using namespace LagSHT;
#define DEBUG 10
#define DEBVEC 0
#define NUM_OMP_THREADS 4
#define CHECK_BY_CCQUAD 0
//JEC 12/1/16: alpha becomes double instead of int
//-------- Parameters set in the main and used in the different test functions
......@@ -499,688 +495,14 @@ void TestJlnpComputation() {
clenshawDir,jlnpDir,recomputeJlnp);
tstack_pop("Ctor lag2bess");
cout << " ______________ TestJlnpComputation End ___________ " << endl;
}
//---------------------------------------------------------------------------------
class JBess : public ClassFunc1D {
public:
JBess(int l, r_8 kt):l_(l),kt_(kt) {norm_ = sqrt(2./M_PI);}
virtual double operator()(double x) const {
return norm_*Bessel::jn(l_,x*kt_);
}
virtual ~JBess() {}
private:
int l_;
r_8 kt_;
r_8 norm_;
};
//------------------
class KLag : public ClassFunc1D {
public:
KLag(LaguerreFuncQuad* Ln, r_8 norm):Ln_(Ln),norm_(norm){}
virtual double operator()(double x) const {
return norm_*x*x*Ln_->Value(x);
}
virtual ~KLag() {}
private:
LaguerreFuncQuad* Ln_; //no ownership
r_8 norm_;
};
//------------------
class IntFunc1D : public ClassFunc1D {
public:
IntFunc1D(int l, r_8 klp,LaguerreFuncQuad* Ln, r_8 norm=1.0):
l_(l),klp_(klp),Ln_(Ln),norm_(norm) {}
virtual double operator()(double x) const {
return norm_*(x*x)*(Bessel::jn(l_,x*klp_))*Ln_->Value(x);
}
virtual ~IntFunc1D() {}
private:
int l_;
r_8 klp_;
LaguerreFuncQuad* Ln_; //no ownership
r_8 norm_;
};//class IntFunc1D
//------------------
void ChebyshevSampling(int n, vector<r_8>& val, const ClassFunc1D& f, r_8 a, r_8 b){
r_8 bma = 0.5*(b-a); r_8 bpa = 0.5*(b+a);
r_8 cte = M_PI/((r_8)n);
int k;
r_8 x;
#pragma omp parallel for private(x) num_threads(NUM_OMP_THREADS)
for(k=0;k<=n;k++){
x = cos(k*cte)*bma+bpa;
val[k] = f(x);
}
}
//------------------
class FFTPlanning {
public:
FFTPlanning(int n, vector<r_8>& vec) {
createPlan(n,vec);
}
~FFTPlanning() { DestroyPlan(); }
void DestroyPlan() {
if(plan_) {
cout << "Destroy Plan ["<< this << "]" << endl;
fftw_destroy_plan(plan_);
plan_=NULL;
}
}
void Execute() const { fftw_execute(plan_); }
private:
void createPlan(int n, vector<r_8>& vec) {
cout << "Create Plan ["<< this << "]" << endl;
plan_ = fftw_plan_r2r_1d(n + 1,&vec.data()[0], &vec.data()[0], FFTW_REDFT00, FFTW_ESTIMATE);
}
fftw_plan plan_;
};
//------------------
//size of val = n+1
//void ChebyshevCoeffFFT(int n, const fftw_plan plan, vector<r_8>& val){
void ChebyshevCoeffFFT(int n, const FFTPlanning& plan, vector<r_8>& val){
tstack_push("FFTW plan exec...");
// fftw_execute ( plan );
plan.Execute();
tstack_pop("FFTW plan exec...");
/*
Chebyshev Coeff. the noramization of FFTW is propto 1/(val.size()-1) = 1/n
*/
tstack_push("FFTW norm...");
r_8 norm = 1./((r_8)n);
transform(val.begin(), val.end(), val.begin(),
std::bind1st(multiplies<r_8>(),norm));
val[0] *= 0.5;
val[n] *= 0.5;
tstack_pop("FFTW norm...");
}
//------------------
void InverseChebyshevCoeffFFT(int n, const FFTPlanning& plan, vector<r_8>& val){
val[0] *= 2.0;
val[n] *= 2.0;
// fftw_execute ( plan );
plan.Execute();
transform(val.begin(), val.end(), val.begin(),
std::bind1st(multiplies<r_8>(),0.5));
}
//------------------
void VectorDebug(const string& msg, const vector<r_8>&vec, int beg, int end){
cout << "(JEC) debug " << msg << endl;
for(int i=beg;i<end;i++) {
cout << vec[i] << ", ";
}