diff --git a/IFFTW/fftwserver.cc b/IFFTW/fftwserver.cc index 027f42568d665175eeba0df1dc17ec57eae4f454..d59c8e754177538e9f56bce93a9dfd078f9cd303 100644 --- a/IFFTW/fftwserver.cc +++ b/IFFTW/fftwserver.cc @@ -2,16 +2,20 @@ #include "FFTW/fftw.h" #include "FFTW/rfftw.h" -class FFTWServerPlan{ + +#define MAXND_FFTW 5 + +class FFTWServerPlan { public: FFTWServerPlan(int n, fftw_direction dir, bool fgreal=false); - FFTWServerPlan(int nx, int ny, fftw_direction dir, bool fgreal=false); + FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal=false); ~FFTWServerPlan(); void Recreate(int n); - void Recreate(int nx, int ny); + void Recreate(int nd, int * nxyz); - int _n; - int _nx, _ny; + int _n; // Array dimension for 1-d arrays + int _nd; // Nb of dimensions for n-d arrays + int _nxyz[MAXND_FFTW]; // Array dimensions for n-d arrays fftw_direction _dir; fftw_plan p; @@ -29,28 +33,34 @@ FFTWServerPlan::FFTWServerPlan(int n, fftw_direction dir, bool fgreal) rp = NULL; pnd = NULL; rpnd = NULL; - _nx = _ny = -10; + _nd = -1; + for(int k=0; k<MAXND_FFTW; k++) _nxyz[k] = -10; _n = n; _dir = dir; if (fgreal) rp = rfftw_create_plan(n, dir, FFTW_ESTIMATE); else p = fftw_create_plan(n, dir, FFTW_ESTIMATE); } -FFTWServerPlan::FFTWServerPlan(int nx, int ny, fftw_direction dir, bool fgreal) + +FFTWServerPlan::FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal) { - if ( (nx < 1) || (ny <1) ) - throw ParmError("FFTWServerPlan: Array size Nx or Ny <= 0 !"); + int k; + if (nd > MAXND_FFTW) + throw ParmError("FFTWServerPlan: Array rank (nd) > MAXND_FFTW !"); p = NULL; rp = NULL; pnd = NULL; rpnd = NULL; _n = -10; - _nx = nx; - _ny = ny; + _nd = nd; + for(k=0; k<nd; k++) { + if (nxyz[k] < 1) + throw ParmError("FFTWServerPlan: One of the Array size <= 0 !"); + _nxyz[k] = nxyz[k]; + } + for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10; _dir = dir; - int sz[2]; - sz[0]= nx; sz[1] = ny; - if (fgreal) rpnd = rfftwnd_create_plan(2, sz, dir, FFTW_ESTIMATE); - else pnd = fftwnd_create_plan(2, sz, dir, FFTW_ESTIMATE); + if (fgreal) rpnd = rfftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE); + else pnd = fftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE); } FFTWServerPlan::~FFTWServerPlan() @@ -66,8 +76,8 @@ FFTWServerPlan::Recreate(int n) { if (n < 1) throw ParmError("FFTWServerPlan::Recreate(n) n < 0 !"); - if ((_nx > 0) || (_ny > 0)) - throw ParmError("FFTWServerPlan::Recreate(n) Nx or Ny > 0 !"); + if (_nd > 0) + throw ParmError("FFTWServerPlan::Recreate(n) Multi-dimensional plan ! > 0 !"); if (n == _n) return; _n = n; if (p) { @@ -81,24 +91,33 @@ FFTWServerPlan::Recreate(int n) } void -FFTWServerPlan::Recreate(int nx, int ny) +FFTWServerPlan::Recreate(int nd, int * nxyz) { - if ( (nx < 1) || (ny <1) ) - throw ParmError("FFTWServerPlan:Recreate(nx, ny) size Nx or Ny <= 0 !"); if (_n > 0) - throw ParmError("FFTWServerPlan::Recreate(nx, ny) N > 0 !"); - if ((nx == _nx) && (ny == _ny)) return; - _nx = nx; - _ny = ny; - int sz[2]; - sz[0]= nx; sz[1] = ny; + throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) 1-dimensional plan !"); + int k; + if (nd == _nd) { + bool samepl = true; + for (int k=0; k<nd; k++) + if (nxyz[k] != _nxyz[k]) samepl = false; + if (samepl) return; + } + if (nd > MAXND_FFTW) + throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) Array rank (nd) > MAXND_FFTW !"); + _nd = nd; + for(k=0; k<nd; k++) { + if (nxyz[k] < 1) + throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) One of the Array size <= 0 !"); + _nxyz[k] = nxyz[k]; + } + for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10; if (pnd) { fftwnd_destroy_plan(pnd); - pnd = fftwnd_create_plan(2, sz,_dir, FFTW_ESTIMATE); + pnd = fftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE); } else { rfftwnd_destroy_plan(rpnd); - rpnd = rfftwnd_create_plan(2, sz, _dir, FFTW_ESTIMATE); + rpnd = rfftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE); } } @@ -107,16 +126,18 @@ FFTWServerPlan::Recreate(int nx, int ny) /* --Methode-- */ FFTWServer::FFTWServer() : FFTServerInterface("FFTServer using FFTW package") + , ckR4(true, false) , ckR8(true, false) + { _p1df = NULL; _p1db = NULL; - _p2df = NULL; - _p2db = NULL; + _pndf = NULL; + _pndb = NULL; _p1drf = NULL; _p1drb = NULL; - _p2drf = NULL; - _p2drb = NULL; + _pndrf = NULL; + _pndrb = NULL; } @@ -125,13 +146,13 @@ FFTWServer::~FFTWServer() { if (_p1df) delete _p1df ; if (_p1db) delete _p1db ; - if (_p2df) delete _p2df ; - if (_p2db) delete _p2db ; + if (_pndf) delete _pndf ; + if (_pndb) delete _pndb ; if (_p1drf) delete _p1drf ; if (_p1drb) delete _p1drb ; - if (_p2drf) delete _p2drf ; - if (_p2drb) delete _p2drb ; + if (_pndrf) delete _pndrf ; + if (_pndrb) delete _pndrb ; } /* --Methode-- */ @@ -142,51 +163,101 @@ FFTServerInterface * FFTWServer::Clone() /* --Methode-- */ void -FFTWServer::FFTForward(TVector< complex<double> > const & in, TVector< complex<double> > & out) +FFTWServer::FFTForward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out) { - if (_p1df) _p1df->Recreate(in.NElts()); - else _p1df = new FFTWServerPlan(in.NElts(), FFTW_FORWARD, false); - out.ReSize(in.NElts()); - fftw_one(_p1df->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) ); - if(this->getNormalize()) out=out/complex<double>(pow(in.NElts(),0.5),0.); + int rank = ckR8.CheckResize(in, out); + if (rank == 1) { // One dimensional transform + if (_p1df) _p1df->Recreate(in.Size()); + else _p1df = new FFTWServerPlan(in.Size(), FFTW_FORWARD, false); + fftw_one(_p1df->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) ); + } + else { // Multi dimensional + if (in.NbDimensions() > MAXND_FFTW) + throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !"); + int sz[MAXND_FFTW]; + int k1 = 0; + int k2 = 0; + for(k1=in.NbDimensions()-1; k1>=0; k1--) { + sz[k2] = in.Size(k1); k2++; + } + if (_pndf) _pndf->Recreate(in.NbDimensions(), sz); + else _pndf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_FORWARD, false); + fftwnd_one(_pndf->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) ); + } + // $CHECK$ Reza 9/2/2001 , Verifier normalisation + if(this->getNormalize()) out=out/complex<r_8>(sqrt((double)in.Size()),0.); + return; } + /* --Methode-- */ -void FFTWServer::FFTBackward(TVector< complex<double> > const & in, TVector< complex<double> > & out) +void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out) { - if (_p1db) _p1db->Recreate(in.NElts()); - else _p1db = new FFTWServerPlan(in.NElts(), FFTW_BACKWARD, false); - out.ReSize(in.NElts()); - fftw_one(_p1db->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) ); - if(this->getNormalize()) out=out/complex<double>(pow(in.NElts(),0.5),0.); - + int rank = ckR8.CheckResize(in, out); + if (rank == 1) { // One dimensional transform + if (_p1db) _p1db->Recreate(in.Size()); + else _p1db = new FFTWServerPlan(in.Size(), FFTW_BACKWARD, false); + fftw_one(_p1db->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) ); + } + else { // Multi dimensional + if (in.NbDimensions() > MAXND_FFTW) + throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !"); + int sz[MAXND_FFTW]; + int k1 = 0; + int k2 = 0; + for(k1=in.NbDimensions()-1; k1>=0; k1--) { + sz[k2] = in.Size(k1); k2++; + } + if (_pndb) _pndb->Recreate(in.NbDimensions(), sz); + else _pndb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_BACKWARD, false); + } + // $CHECK$ Reza 9/2/2001 , Verifier normalisation + if(this->getNormalize()) out=out/complex<r_8>(sqrt((double)in.Size()),0.); + return; } -void FFTWServer::FFTForward(TVector< double > const & in, TVector< complex<double> > & out) +void FFTWServer::FFTForward(TArray< r_8 > const & in, TArray< complex<r_8> > & out) { - int size = in.NElts()/2; - - if(in.NElts()%2 != 0) { size = in.NElts()/2 +1;} - else { size = in.NElts()/2 +1 ;} + int rank = ckR8.CheckResize(in, out); + TArray<r_8> outtemp(in, false); + + if (rank == 1) { // One dimensional transform + if (_p1drf) _p1drf->Recreate(in.Size()); + else _p1drf = new FFTWServerPlan(in.Size(), FFTW_REAL_TO_COMPLEX, true); + rfftw_one(_p1drf->rp, (fftw_real *)(in.Data()) , (fftw_real *)(outtemp.Data())); + ReShapetoCompl(outtemp, out); + } + else { // Multi dimensional + if (in.NbDimensions() > MAXND_FFTW) + throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !"); + int sz[MAXND_FFTW]; + int k1 = 0; + int k2 = 0; + for(k1=in.NbDimensions()-1; k1>=0; k1--) { + sz[k2] = in.Size(k1); k2++; + } + if (_pndrf) _pndrf->Recreate(in.NbDimensions(), sz); + else _pndrf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_REAL_TO_COMPLEX, true); + rfftwnd_one_real_to_complex(_pndrf->rpnd, (fftw_real *)(in.Data()) , + (fftw_complex *)(out.Data()) ); + } + // $CHECK$ Reza 9/2/2001 , Verifier normalisation + if(this->getNormalize()) out=out/complex<r_8>(sqrt((double)in.Size()),0.); + return; - TVector< double > const outTemp(in.NElts()); - out.ReSize(size); - if (_p1drf) _p1drf->Recreate(in.NElts()); - else _p1drf = new FFTWServerPlan(in.NElts(), FFTW_REAL_TO_COMPLEX, true); - rfftw_one(_p1drf->rp, (fftw_real *)(in.Data()) , (fftw_real *)(outTemp.Data())); - ReShapetoCompl(outTemp, out); - if(this->getNormalize()) out=out/complex<double>(pow(in.NElts(),0.5),0.); } -void FFTWServer::FFTBackward(TVector< complex<double> > const & in, TVector< double > & out) +void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< r_8 > & out) { + throw ParmError("FFTWServer::FFTBackward(TArray< complex<r_8> > ... Not implemented ... !"); + /* int size; if(in(in.NElts()).imag() == 0) { size = 2*in.NElts()-2;} else { size = 2*in.NElts()-1;} - TVector< double > inTemp(size); + TArray< r_8 > inTemp(size); out.ReSize(size); if (_p1drb) _p1drb->Recreate(size); @@ -195,80 +266,81 @@ void FFTWServer::FFTBackward(TVector< complex<double> > const & in, TVector< dou ReShapetoReal(in, inTemp); rfftw_one(_p1drb->rp, (fftw_real *)(inTemp.Data()) , (fftw_real *)(out.Data())); if(this->getNormalize()) out=out/pow(size,0.5); + */ } -/* --Methode-- */ -void FFTWServer::FFTForward(TMatrix< complex<double> > const & in, TMatrix< complex<double> > & out) +/* +void FFTWServer::FFTForward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out) { out.ReSize(in.NRows(),in.NCols()); - if (_p2df) _p2df->Recreate( in.NRows(),in.NCols()); - else _p2df = new FFTWServerPlan( in.NCols(),in.NRows(), FFTW_FORWARD, false); + if (_pndf) _pndf->Recreate( in.NRows(),in.NCols()); + else _pndf = new FFTWServerPlan( in.NCols(),in.NRows(), FFTW_FORWARD, false); - fftwnd_one(_p2df->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) ); - if(this->getNormalize()) out=out/complex<double>(pow(in.NRows()*in.NCols(),0.5),0.); + fftwnd_one(_pndf->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) ); + if(this->getNormalize()) out=out/complex<r_8>(pow(in.NRows()*in.NCols(),0.5),0.); } -/* --Methode-- */ -void FFTWServer::FFTBackward(TMatrix< complex<double> > const & in, TMatrix< complex<double> > & out) +void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out) { - if (_p2db) _p2db->Recreate(in.NCols(), in.NRows()); - else _p2db = new FFTWServerPlan(in.NCols(), in.NRows(), FFTW_BACKWARD, false); + if (_pndb) _pndb->Recreate(in.NCols(), in.NRows()); + else _pndb = new FFTWServerPlan(in.NCols(), in.NRows(), FFTW_BACKWARD, false); out.ReSize(in.NRows(), in.NCols()); - fftwnd_one(_p2db->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) ); - if(this->getNormalize()) out=out/complex<double>(pow(in.NRows()*in.NCols(),0.5),0.); + fftwnd_one(_pndb->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) ); + if(this->getNormalize()) out=out/complex<r_8>(pow(in.NRows()*in.NCols(),0.5),0.); } -/* --Methode-- */ -void FFTWServer::FFTForward(TMatrix< double > const & in, TMatrix< complex<double> > & out) +void FFTWServer::FFTForward(TArray< r_8 > const & in, TArray< complex<r_8> > & out) { - TMatrix< double > inNew(in.NCols(),in.NRows()); + TArray< r_8 > inNew(in.NCols(),in.NRows()); for(int i=0; i<in.NRows(); i++) for(int j=0;j<in.NCols(); j++) inNew(j,i) = in(i,j); - if (_p2drf) _p2drf->Recreate(inNew.NRows(),inNew.NCols()); - else _p2drf = new FFTWServerPlan(inNew.NRows(), inNew.NCols(),FFTW_REAL_TO_COMPLEX, true); + if (_pndrf) _pndrf->Recreate(inNew.NRows(),inNew.NCols()); + else _pndrf = new FFTWServerPlan(inNew.NRows(), inNew.NCols(),FFTW_REAL_TO_COMPLEX, true); // rfftwnd_plan p; - TMatrix< complex<double> > outTemp; + TArray< complex<r_8> > outTemp; outTemp.ReSize(in.NRows(),in.NCols()); - rfftwnd_one_real_to_complex(_p2drf->rpnd, (fftw_real *)(in.Data()) , (fftw_complex *)(out.Data()) ); + rfftwnd_one_real_to_complex(_pndrf->rpnd, (fftw_real *)(in.Data()) , (fftw_complex *)(out.Data()) ); } -/* --Methode-- */ -void FFTWServer::FFTBackward(TMatrix< complex<double> > const & in, TMatrix< double > & out) +void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< r_8 > & out) { - TMatrix< complex<double> > inNew(in.NCols(),in.NRows()); + TArray< complex<r_8> > inNew(in.NCols(),in.NRows()); for(int i=0; i<in.NRows(); i++) for(int j=0;j<in.NCols(); j++) inNew(j,i) = in(i,j); - if (_p2drb) _p2drb->Recreate(inNew.NRows(),inNew.NCols()); - else _p2drb = new FFTWServerPlan(inNew.NRows(), inNew.NCols(),FFTW_COMPLEX_TO_REAL, true); + if (_pndrb) _pndrb->Recreate(inNew.NRows(),inNew.NCols()); + else _pndrb = new FFTWServerPlan(inNew.NRows(), inNew.NCols(),FFTW_COMPLEX_TO_REAL, true); // rfftwnd_plan p; out.ReSize(in.NRows(),in.NCols()); - rfftwnd_one_complex_to_real(_p2drb->rpnd, (fftw_complex *)(in.Data()) , (fftw_real *)(out.Data()) ); + rfftwnd_one_complex_to_real(_pndrb->rpnd, (fftw_complex *)(in.Data()) , (fftw_real *)(out.Data()) ); cout << " in the function !!!" << endl; if(this->getNormalize()) { - double norm = (double)(in.NRows()*in.NCols()); + r_8 norm = (r_8)(in.NRows()*in.NCols()); out=out/norm; } } +*/ + /* --Methode-- */ -void FFTWServer::ReShapetoReal( TVector< complex<double> > const & in, TVector< double > & out) +void FFTWServer::ReShapetoReal( TArray< complex<r_8> > const & in, TArray< r_8 > & out) { - int N = in.NElts(); + int N = in.Size(); + /* int i; - if (in(in.NElts()).imag() == 0) + if (in(in.Size()).imag() == 0) { out(0) = in(0).real(); for(i=1; i<in.NElts(); i++) @@ -292,31 +364,27 @@ void FFTWServer::ReShapetoReal( TVector< complex<double> > const & in, TVector< out(i+in.NElts()-1) = in(in.NElts()-i).imag(); } } + */ + out[0] = in[0].real(); + int k=0; + for(int i=1; i<N; i++) { + out[i] = in[i].real(); + out[N-i] = in[i].imag(); + } + // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl; } /* --Methode-- */ -void FFTWServer::ReShapetoCompl(TVector< double > const & in, TVector< complex<double> > & out) +void FFTWServer::ReShapetoCompl(TArray< r_8 > const & in, TArray< complex<r_8> > & out) { - int N = in.NElts(); + int N = in.Size(); // for(int k=0; k<in.NElts(); k++) cout << "ReShapetoCompl in " << k << " " << in(k) << endl; - out(0) = complex<double> (in(0),0.); - if(in.NElts()%2 !=0) - { - for(int k=1;k<=N/2+1;k++) - { - out(k) = complex<double> (in(k),in(N-k)); - } - } - else - { - for(int k=1;k<N/2;k++) - { - out(k) = complex<double> (in(k),in(N-k)); - } - out(N/2) = complex<double> (in(N/2),0.); - } + out[0] = complex<r_8> (in[0],0.); + for(int k=1; k<N+1/2; k++) + out[k] = complex<r_8>(in[k], in[N-k]); + if (N%2 == 0) out[N/2] = complex<r_8>(in[N/2], 0.); // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoCompl out " << k << " " << out(k) << endl; } diff --git a/IFFTW/fftwserver.h b/IFFTW/fftwserver.h index e775e8c0ec04bf9ab6ea0a7fd1af90be3ae6e252..fbbf935ae5a912ade6afb52bc9485a640712e96b 100644 --- a/IFFTW/fftwserver.h +++ b/IFFTW/fftwserver.h @@ -18,32 +18,29 @@ class FFTWServer : public FFTServerInterface { virtual FFTServerInterface * Clone(); - // Transformee unidimensionnelle - virtual void FFTForward(TVector< complex<double> > const & in, TVector< complex<double> > & out); - virtual void FFTBackward(TVector< complex<double> > const & in, TVector< complex<double> > & out); - virtual void FFTForward(TVector< double > const & in, TVector< complex<double> > & out); - virtual void FFTBackward(TVector< complex<double> > const & in, TVector< double > & out); - - // Transforme 2D - virtual void FFTForward(TMatrix< complex<double> > const & in, TMatrix< complex<double> > & out); - virtual void FFTBackward(TMatrix< complex<double> > const & in, TMatrix< complex<double> > & out); - virtual void FFTForward(TMatrix< double > const & in, TMatrix< complex<double> > & out); - virtual void FFTBackward(TMatrix< complex<double> > const & in, TMatrix< double > & out); + // Transforme unidimensionnelle , N-dimensionnel + virtual void FFTForward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out); + virtual void FFTBackward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out); + virtual void FFTForward(TArray< r_8 > const & in, TArray< complex<r_8> > & out); + virtual void FFTBackward(TArray< complex<r_8> > const & in, TArray< r_8 > & out); // Methodes statiques pour reordonner les donnees - virtual void ReShapetoReal( TVector< complex<r_8> > const & in, TVector< r_8 > & out); - virtual void ReShapetoCompl(TVector< r_8 > const & in, TVector< complex<r_8> > & out); + virtual void ReShapetoReal( TArray< complex<r_8> > const & in, TArray< r_8 > & out); + virtual void ReShapetoCompl(TArray< r_8 > const & in, TArray< complex<r_8> > & out); protected: FFTWServerPlan * _p1df; FFTWServerPlan * _p1db; - FFTWServerPlan * _p2df; - FFTWServerPlan * _p2db; + FFTWServerPlan * _pndf; + FFTWServerPlan * _pndb; FFTWServerPlan * _p1drf; FFTWServerPlan * _p1drb; - FFTWServerPlan * _p2drf; - FFTWServerPlan * _p2drb; + FFTWServerPlan * _pndrf; + FFTWServerPlan * _pndrb; + + FFTArrayChecker<r_4> ckR4; + FFTArrayChecker<r_8> ckR8; }; #endif