Changeset 3000 in Sophya for trunk/SophyaExt/IFFTW/fftwserver.cc


Ignore:
Timestamp:
Jul 3, 2006, 1:02:46 PM (19 years ago)
Author:
ansari
Message:

Passage a FFTW3 + suppression const des arguments FFTForward/Backward cmv+Reza 03/07/2006

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaExt/IFFTW/fftwserver.cc

    r2994 r3000  
    11#include "sopnamsp.h"
    22#include "fftwserver.h"
    3 #include "FFTW/fftw.h"
    4 #include "FFTW/rfftw.h"
    53
    64/*!
    … …  
    3634#define MAXND_FFTW 5
    3735
    38 namespace SOPHYA {
     36#ifdef FFTW_V2_EXTSOP
     37#include "fftw2server.cc"
     38#else
     39#include "fftw3server.cc"
     40#endif
    3941
    40 class FFTWServerPlan {
    41 public:
    42   FFTWServerPlan(int n, fftw_direction dir, bool fgreal=false);
    43   FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal=false);
    44   ~FFTWServerPlan();
    45   void Recreate(int n);
    46   void Recreate(int nd, int * nxyz);
    47 
    48   static void FillSizes(const BaseArray & in, int * sz);
    49 
    50   int _n;  // Array dimension for 1-d arrays
    51   int _nd; // Nb of dimensions for n-d arrays
    52   int _nxyz[MAXND_FFTW]; // Array dimensions for n-d arrays
    53   fftw_direction _dir;
    54 
    55   fftw_plan p;
    56   rfftw_plan rp;
    57   fftwnd_plan pnd;
    58   rfftwnd_plan rpnd;
    59    
    60 };
    61 
    62 } // Fin du namespace
    63 
    64 FFTWServerPlan::FFTWServerPlan(int n, fftw_direction dir, bool fgreal)
    65 {
    66   if (n < 1)
    67     throw ParmError("FFTWServerPlan: Array size <= 0 !");
    68   p = NULL;
    69   rp = NULL;
    70   pnd = NULL;
    71   rpnd = NULL;
    72   _nd = -1;
    73   for(int k=0; k<MAXND_FFTW; k++) _nxyz[k] = -10;
    74   _n = n; 
    75   _dir = dir;
    76   if (fgreal) {
    77     rp = rfftw_create_plan(n, dir, FFTW_ESTIMATE);
    78     if (rp == NULL)
    79       throw AllocationError("FFTWServerPlan: failed to create plan (rp) !");
    80   }
    81   else {
    82     p = fftw_create_plan(n, dir, FFTW_ESTIMATE);
    83     if (p == NULL)
    84       throw AllocationError("FFTWServerPlan: failed to create plan (p) !");
    85   }
    86    
    87 }
    88 
    89 FFTWServerPlan::FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal)
    90 {
    91   int k;
    92   if (nd > MAXND_FFTW)
    93     throw ParmError("FFTWServerPlan: Array rank (nd) > MAXND_FFTW !");
    94   p = NULL;
    95   rp = NULL;
    96   pnd = NULL;
    97   rpnd = NULL;
    98   _n = -10;
    99   _nd = nd;
    100   for(k=0; k<nd; k++) {
    101     if (nxyz[k] < 1)
    102        throw ParmError("FFTWServerPlan: One of the Array size <= 0 !");   
    103     _nxyz[k] = nxyz[k];
    104   }
    105   for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
    106   _dir = dir;
    107   if (fgreal) {
    108     rpnd = rfftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
    109     if (rpnd == NULL)
    110       throw AllocationError("FFTWServerPlan: failed to create plan (rpnd) !");
    111   }
    112   else {
    113     pnd = fftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
    114     if (pnd == NULL)
    115       throw AllocationError("FFTWServerPlan: failed to create plan (pnd) !");
    116   }
    117 
    118 }
    119 
    120 FFTWServerPlan::~FFTWServerPlan()
    121 {
    122   if (p) fftw_destroy_plan(p);
    123   if (rp) rfftw_destroy_plan(rp);
    124   if (pnd) fftwnd_destroy_plan(pnd);
    125   if (rpnd) rfftwnd_destroy_plan(rpnd);
    126 }
    127 
    128 void
    129 FFTWServerPlan::Recreate(int n)
    130 {
    131   if (n < 1)
    132    throw ParmError("FFTWServerPlan::Recreate(n) n < 0 !");
    133   if (_nd > 0) 
    134    throw ParmError("FFTWServerPlan::Recreate(n) Multi-dimensional plan ! > 0 !");
    135   if (n == _n) return;
    136   _n = n; 
    137   if (p) {
    138     fftw_destroy_plan(p);
    139     p = fftw_create_plan(n, _dir, FFTW_ESTIMATE);
    140     if (p == NULL)
    141       throw AllocationError("FFTWServerPlan::Recreate failed to create plan (p) !");
    142   }
    143   else {
    144     rfftw_destroy_plan(rp);
    145     rp = rfftw_create_plan(n, _dir, FFTW_ESTIMATE);
    146     if (rp == NULL)
    147       throw AllocationError("FFTWServerPlan::Recreate failed to create plan (rp) !");
    148   }
    149 
    150 }
    151 
    152 void
    153 FFTWServerPlan::Recreate(int nd, int * nxyz)
    154 {
    155   if (_n > 0)
    156     throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) 1-dimensional plan !");
    157   int k;
    158   if (nd == _nd) {
    159     bool samepl = true;
    160     for (int k=0; k<nd; k++)
    161       if (nxyz[k] != _nxyz[k]) samepl = false;
    162     if (samepl) return;
    163   }
    164   if (nd > MAXND_FFTW)
    165     throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) Array rank (nd) > MAXND_FFTW !");
    166   _nd = nd;
    167   for(k=0; k<nd; k++) {
    168     if (nxyz[k] < 1)
    169        throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) One of the Array size <= 0 !");   
    170     _nxyz[k] = nxyz[k];
    171   }
    172   for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
    173   if (pnd) {
    174     fftwnd_destroy_plan(pnd);
    175     pnd = fftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
    176     if (pnd == NULL)
    177       throw AllocationError("FFTWServerPlan::Recreate failed to create plan (pnd) !");
    178   }
    179   else {
    180     rfftwnd_destroy_plan(rpnd);
    181     rpnd = rfftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
    182     if (rpnd == NULL)
    183       throw AllocationError("FFTWServerPlan::Recreate failed to create plan (rpnd) !");
    184   }
    185 
    186 }
    187 
    188 void
    189 FFTWServerPlan::FillSizes(const BaseArray & in, int * sz)
    190 {
    191   int k1 = 0;
    192   int k2 = 0;
    193   for(k1=in.NbDimensions()-1; k1>=0; k1--) {
    194     sz[k2] = in.Size(k1); k2++;
    195   }
    196 }
    197 
    198 /* --Methode-- */
    199 FFTWServer::FFTWServer()
    200   : FFTServerInterface("FFTServer using FFTW package")
    201   , ckR8("FFTWServer - ", true, false)
    202 
    203 {
    204   _p1df = NULL;
    205   _p1db = NULL;
    206   _pndf = NULL;
    207   _pndb = NULL;
    208 
    209   _p1drf = NULL;
    210   _p1drb = NULL;
    211   _pndrf = NULL;
    212   _pndrb = NULL;
    213 }
    214 
    215 
    216 /* --Methode-- */
    217 FFTWServer::~FFTWServer()
    218 {
    219   if (_p1df) delete _p1df ;
    220   if (_p1db) delete _p1db ;
    221   if (_pndf) delete _pndf ;
    222   if (_pndb) delete _pndb ;
    223 
    224   if (_p1drf) delete _p1drf ;
    225   if (_p1drb) delete _p1drb ;
    226   if (_pndrf) delete _pndrf ;
    227   if (_pndrb) delete _pndrb ;
    228 }
    229 
    230 /* --Methode-- */
    231 FFTServerInterface * FFTWServer::Clone()
    232 {
    233   return (new FFTWServer) ;
    234 }
    235 
    236 /* --Methode-- */
    237 void
    238 FFTWServer::FFTForward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out)
    239 {
    240   int rank = ckR8.CheckResize(in, out);
    241   if (rank == 1) { // One dimensional transform
    242     if (_p1df) _p1df->Recreate(in.Size());
    243     else _p1df = new FFTWServerPlan(in.Size(), FFTW_FORWARD, false);
    244     fftw_one(_p1df->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
    245   }
    246   else {   // Multi dimensional
    247     if (in.NbDimensions() > MAXND_FFTW)
    248       throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
    249     int sz[MAXND_FFTW];
    250     FFTWServerPlan::FillSizes(in, sz);
    251     //    int k1 = 0;
    252     //    int k2 = 0;
    253     //    for(k1=in.NbDimensions()-1; k1>=0; k1--) {
    254     //      sz[k2] = in.Size(k1); k2++;
    255     //    }
    256     if (_pndf) _pndf->Recreate(in.NbDimensions(), sz);
    257     else _pndf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_FORWARD, false);
    258     fftwnd_one(_pndf->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
    259   } 
    260   if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
    261   return;
    262 }
    263 
    264 /* --Methode-- */
    265 void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out)
    266 {
    267   int rank = ckR8.CheckResize(in, out);
    268   if (rank == 1) { // One dimensional transform
    269     if (_p1db) _p1db->Recreate(in.Size());
    270     else _p1db = new FFTWServerPlan(in.Size(), FFTW_BACKWARD, false);
    271     fftw_one(_p1db->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
    272   }
    273   else {   // Multi dimensional
    274     if (in.NbDimensions() > MAXND_FFTW)
    275       throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
    276     int sz[MAXND_FFTW];
    277     FFTWServerPlan::FillSizes(in, sz);
    278     //    int k1 = 0;
    279     //    int k2 = 0;
    280     //    for(k1=in.NbDimensions()-1; k1>=0; k1--) {
    281     //      sz[k2] = in.Size(k1); k2++;
    282     //    }
    283     if (_pndb) _pndb->Recreate(in.NbDimensions(), sz);
    284     else _pndb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_BACKWARD, false);
    285     fftwnd_one(_pndb->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
    286   }
    287 
    288   return;
    289 }
    290 
    291 
    292 void FFTWServer::FFTForward(TArray< r_8 > const & in, TArray< complex<r_8> > & out)
    293 { 
    294   int rank = ckR8.CheckResize(in, out);
    295   if (rank == 1) { // One dimensional transform
    296     if (_p1drf) _p1drf->Recreate(in.Size());
    297     else _p1drf = new FFTWServerPlan(in.Size(), FFTW_REAL_TO_COMPLEX, true);
    298     TArray<r_8> outtemp;
    299     outtemp.ReSize(in);
    300     rfftw_one(_p1drf->rp, (fftw_real *)(in.Data()) , (fftw_real *)(outtemp.Data()));
    301     ReShapetoCompl(outtemp, out);
    302   }
    303   else {   // Multi dimensional
    304     if (in.NbDimensions() > MAXND_FFTW)
    305       throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
    306     int sz[MAXND_FFTW];
    307     FFTWServerPlan::FillSizes(in, sz);
    308     //    int k1 = 0;
    309     //    int k2 = 0;
    310     //   for(k1=in.NbDimensions()-1; k1>=0; k1--) {
    311     //      sz[k2] = in.Size(k1); k2++;
    312     //    }
    313     if (_pndrf) _pndrf->Recreate(in.NbDimensions(), sz);
    314     else _pndrf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_REAL_TO_COMPLEX, true);
    315     rfftwnd_one_real_to_complex(_pndrf->rpnd, (fftw_real *)(in.Data()) ,
    316                                 (fftw_complex *)(out.Data()) );
    317   }
    318   if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
    319   return;
    320 
    321 }
    322 
    323 
    324 
    325 void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< r_8 > & out,
    326                              bool usoutsz)
    327 {
    328  
    329   int rank = ckR8.CheckResize(in, out, usoutsz);
    330   if (rank == 1) { // One dimensional transform
    331     TArray<r_8> intemp;
    332     intemp.ReSize(out);
    333     if (_p1drb) _p1drb->Recreate(out.Size());
    334     else _p1drb = new FFTWServerPlan(out.Size(), FFTW_COMPLEX_TO_REAL, true);
    335 
    336     ReShapetoReal(in, intemp, usoutsz);
    337     //    cerr << " DEBUG-FFTWServer::FFTBackward() in = \n" << in << endl;
    338     //    cerr << " DEBUG-FFTWServer::FFTBackward() intemp = \n" << intemp << endl;
    339     rfftw_one(_p1drb->rp, (fftw_real *)(intemp.Data()) , (fftw_real *)(out.Data()));
    340 
    341   }
    342   else {   // Multi dimensional
    343     if (in.NbDimensions() > MAXND_FFTW)
    344       throw ParmError("FFTWServer::FFTForward( complex<r_8>, r_8 ) rank > MAXND_FFTW !");
    345     int sz[MAXND_FFTW];
    346     FFTWServerPlan::FillSizes(out, sz);
    347     //    int k1 = 0;
    348     //    int k2 = 0;
    349     //    for(k1=out.NbDimensions()-1; k1>=0; k1--) {
    350     //      sz[k2] = out.Size(k1); k2++;
    351     //    }
    352     if (_pndrb) _pndrb->Recreate(in.NbDimensions(), sz);
    353     else _pndrb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_COMPLEX_TO_REAL, true);
    354 
    355     rfftwnd_one_complex_to_real(_pndrb->rpnd, (fftw_complex *)(in.Data()) , (fftw_real *)(out.Data()) );
    356   }
    357   return;
    358 }
    359 
    360 
    361 
    362 /* --Methode-- */
    363 void FFTWServer::ReShapetoReal(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa,
    364                                bool usz)
    365 {
    366   TVector< complex<r_8> > in(ina);
    367   TVector< r_8> out(outa);
    368   int n = in.NElts();
    369   r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
    370   sa_size_t ncs;
    371   if (usz) {
    372     if ( (out.NElts() != 2*n-1) && (out.NElts() != 2*n-2) )
    373       throw SzMismatchError("FFTWServer::ReShapetoReal(..., true) - Wrong output array size ");
    374     ncs = out.NElts();
    375   }
    376   else {
    377     ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
    378                     ncs = 2*n-1 : ncs = 2*n-2;
    379 
    380     if (out.NElts() != ncs) {
    381       cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
    382            << " out.NElts()= " << out.NElts() << endl;
    383       throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
    384     }
    385   }
    386   //  if (ncs == 2*n-2)  {
    387   out(0) = in(0).real();
    388   for(int k=1; k<n; k++) {
    389       out(k) = in(k).real();
    390       out(ncs-k) = in(k).imag();
    391   }
    392   if (ncs == 2*n-2)
    393     out(n-1) = in(n-1).real();
    394  
    395   //  for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
    396 }
    397 
    398 
    399 /* --Methode-- */
    400 void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
    401 {
    402   TVector< r_8> in(ina);
    403   TVector< complex<r_8> > out(outa);
    404 
    405   sa_size_t n = in.NElts();
    406   sa_size_t ncs = n/2+1;
    407   sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
    408   if (out.NElts() != ncs)
    409     throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
    410 
    411   out(0) = complex<r_8> (in(0),0.);
    412   for(int k=1; k<ncs; k++)
    413     out(k) = complex<r_8>(in(k), in(n-k));
    414   if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);
    415 
    416 }
    417 
Note: See TracChangeset for help on using the changeset viewer.