source: Sophya/trunk/SophyaProg/Examples/apcxx.icc@ 1673

/*  ----------------- apcxx.icc -------------------

   Example illustrating C++ scientific programming 
   using SOPHYA library arrays and FFT computation 
   through FFTServerInterface.                         
                                   R. Ansari  08/2001 
   ---- Computation steps :
    > Create a matrix (NL x NC )  (mtx)
    > fill it with a gaussian distributed random values
    > make a copy of the matrix   (mtxs)
    > create a 1D filter in Fourier space 
    > Loop over matrix rows k
     >> Extract row k (fline)
     >> compute 1D Fourier transform FFTForward (fline) 
     >> apply filter in Fourier space
     >> compute backward 1D FFT 
     >> Replace matrix row with the filtered values
   

   this example code can be 
     - included in a main program
     - executed using runcxx 
       csh> runcxx -tmpdir /tmp -f apcxx.icc
     - executed within spiapp 
       Cmd> c++execfrf apcxx.icc
*/

// Select computation on float or double (r_4 r_8)
#define FTYP r_4

// Number of matrix lines and colums
int NL = 1024;
int NC = 4096;
cout << " apc_xx : NL= " << NL << " NC= " << NC << endl;
PrtTim("apcxx_Start");

// BaseArray::SetDefaultMemoryMapping(BaseArray::CMemoryMapping);
// BaseArray::SetMaxPrint(10, 3);

// Creation of the initial matrix
TMatrix< FTYP > mtx(NL, NC), mtxs;
// Filling matrix with gaussian random values
mtx = RandomSequence(RandomSequence::Gaussian, 15., 3.);
// Making a copy of the original matrix
mtxs = mtx;

// Creation and initialization of the Fourier filter filt(nu) = 1/(1+0.3*nu)
int LFFT = NC/2+1;
TVector< complex< FTYP > > filt(LFFT, BaseArray::RowVector);
filt(0) = 1.;
for(int i=1; i<filt.Size(); i++) 
  filt(i) = 1./(1+0.3*(double)i);


// Creation of the FFTServer
FFTPackServer ffts;
ffts.setNormalize(true);

PrtTim("apcxx_AfterInit");

// Vectors for FFT operations 
TVector< FTYP > fline(NC, BaseArray::RowVector);
TVector< complex< FTYP > > vfft;

//  Loop over matrix rows
for(int k=0; k<NL; k++) {
  fline = mtx.Row(k);      // Matrix row extraction 
  ffts.FFTForward(fline, vfft);    // Compute 1D forward FFT
  // Applying filter in Fourier space f(nu) =  f(nu)*filter(nu)
  vfft.MulElt(filt);             
  ffts.FFTBackward(vfft, fline);   // backward FFT
  mtx.Row(k) = fline;      // replace matrix row with filtered values
 }

PrtTim("apcxx_AfterFFTLoop");

// Macro KeepObj can be used with runcxx or within (s)piapp 
// KeepObj(mtx);
// KeepObj(mtxs);