/*  --------------- apcxx_cstyle.icc -----------------

   Example illustrating C++ scientific programming 
   using SOPHYA library arrays and FFT computation 
   through FFTServerInterface.   
   This example, performs the same operations as in
   apcxx.icc, through C-style programming with 
   explicit loops. This C_style should be avoided when 
   programming in C++.
                      
                                   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_cstyle.icc
     - executed within spiapp 
       Cmd> c++execfrf apcxx_cstyle.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;
int i, j;
cout << " apc_cstyle : NL= " << NL << " NC= " << NC << endl;
PrtTim("apc_cstyle_Start");

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


// Creation of the initial matrix
TMatrix< FTYP > mtx(NL, NC);
// Filling matrix with gaussian random values
// mtx = RandomSequence(RandomSequence::Gaussian, 15., 3.);
for(i=0; i<NL; i++) 
  for(j=0; j<NC; j++) mtx(i, j) = 15.+NorRand()*3.;

// Making a copy of the original matrix
//mtxs = mtx;
TMatrix< FTYP > mtxs(NL, NC);
for(i=0; i<NL; i++) 
  for(j=0; j<NC; j++) mtxs(i, j) = mtx(i, j);


// 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("apc_cstyle_AfterInit");

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

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

PrtTim("apc_cstyle_AfterFFTLoop");

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