Changeset 1395 in Sophya for trunk/SophyaExt

Timestamp:

Feb 12, 2001, 6:09:47 PM (25 years ago)

Author:

ansari

Message:

Corrections- adaptation interface FFTServer - Reza 12/2/2001

Location:

trunk/SophyaExt/IFFTW

Files:

• fftwserver.cc (modified) (7 diffs)
• fftwserver.h (modified) (2 diffs)

trunk/SophyaExt/IFFTW/fftwserver.cc

@@ -127 +127 @@
 FFTWServer::FFTWServer()
   : FFTServerInterface("FFTServer using FFTW package")
-  ,  ckR4(true, false) , ckR8(true, false)
+  , ckR8("FFTWServer - ", true, false)
 
 {
@@ -185 +185 @@
     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.);
+  if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
   return;
 }
@@ -211 +210 @@
     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;
 }
@@ -242 +240 @@
                                 (fftw_complex *)(out.Data()) );
   }
-  // $CHECK$ Reza 9/2/2001 , Verifier normalisation  
-  if(this->getNormalize()) out=out/complex<r_8>(sqrt((double)in.Size()),0.);
+  if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
   return;
 
@@ -252 +249 @@
 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;}
-  
-  TArray< r_8 > inTemp(size);
-  out.ReSize(size);
-
-  if (_p1drb) _p1drb->Recreate(size);
-  else _p1drb = new FFTWServerPlan(size, FFTW_COMPLEX_TO_REAL, true);
-
-  ReShapetoReal(in, inTemp);
-  rfftw_one(_p1drb->rp, (fftw_real *)(inTemp.Data()) , (fftw_real *)(out.Data()));
-  if(this->getNormalize()) out=out/pow(size,0.5);
-  */
+  
+  int rank = ckR8.CheckResize(in, out);
+  if (rank == 1) { // One dimensional transform 
+    TArray<r_8> intemp(out, false); 
+    if (_p1drb) _p1drb->Recreate(out.Size());
+    else _p1drb = new FFTWServerPlan(out.Size(), FFTW_COMPLEX_TO_REAL, true);
+
+    ReShapetoReal(in, intemp);
+    cerr << " DEBUG-FFTWServer::FFTBackward() in = \n" << in << endl;
+    cerr << " DEBUG-FFTWServer::FFTBackward() intemp = \n" << intemp << endl;
+    rfftw_one(_p1drb->rp, (fftw_real *)(intemp.Data()) , (fftw_real *)(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 (_pndrb) _pndrb->Recreate(in.NbDimensions(), sz);
+    else _pndrb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_COMPLEX_TO_REAL, true);
+
+    rfftwnd_one_complex_to_real(_pndrb->rpnd, (fftw_complex *)(in.Data()) , (fftw_real *)(out.Data()) );
+  }
+  return;
 }
 
@@ -336 +345 @@
 
 /* --Methode-- */
-void FFTWServer::ReShapetoReal( TArray< complex<r_8> > const & in, TArray< r_8 >  & out)
-{
-  int N = in.Size();
-  /*
-  int i;
-  if (in(in.Size()).imag() == 0) 
-    {
-      out(0) = in(0).real();
-      for(i=1; i<in.NElts(); i++)
-        {
-          out(i) = in(i).real();
-        }
-      for(i=1; i<in.NElts(); i++)
-        {
-          out(i+in.NElts()-1) = in(in.NElts()-i-1).imag();
-        }
-    }
-  else
-    {
-      out(0) = in(0).real();
-      for(i=1; i<in.NElts(); i++)
-        {
-          out(i) = in(i).real();
-        }
-      for(i=1; i<in.NElts(); i++)
-        {
-          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();
-  }
+void FFTWServer::ReShapetoReal(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa)
+{
+  TVector< complex<r_8> > in(ina);
+  TVector< r_8> out(outa);
+  int n = in.NElts();
+  r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
+  sa_size_t ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
+                    ncs = 2*n-1 : ncs = 2*n-2;
+
+  if (out.NElts() != ncs) {
+    cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs 
+         << " out.NElts()= " << out.NElts() << endl;
+    throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
+  }
+  //  if (ncs == 2*n-2)  {
+  out(0) = in(0).real();
+  for(int k=1; k<n; k++) {
+      out(k) = in(k).real();
+      out(ncs-k) = in(k).imag();
+  }
+  if (ncs == 2*n-2)
+    out(ncs-1) = in(n-1).real();
   
   //  for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
@@ -378 +373 @@
 
 /* --Methode-- */
-void FFTWServer::ReShapetoCompl(TArray< r_8 > const & in, TArray< complex<r_8> > & out)
-{
-  int N = in.Size();
-  //  for(int k=0; k<in.NElts(); k++) cout << "ReShapetoCompl in " << k << " "  << in(k) << endl;
-  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;
-}
-
+void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
+{
+  TVector< r_8> in(ina);
+  TVector< complex<r_8> > out(outa);
+
+  sa_size_t n = in.NElts();
+  sa_size_t ncs = n/2+1;
+  sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
+  if (out.NElts() != ncs)
+    throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
+
+  out(0) = complex<r_8> (in(0),0.);
+  for(int k=1; k<ncs; k++) 
+    out(k) = complex<r_8>(in(k), in(n-k));
+  if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);
+
+}
+

trunk/SophyaExt/IFFTW/fftwserver.h

@@ -26 +26 @@
 
 // Methodes statiques pour reordonner les donnees 
-  virtual void ReShapetoReal( TArray< complex<r_8> > const & in, TArray< 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);
 
@@ -40 +40 @@
  FFTWServerPlan * _pndrb;
 
- FFTArrayChecker<r_4> ckR4;
  FFTArrayChecker<r_8> ckR8;
 };