Complex Multiplication and convolution with FFTW

Hi!
I am currently use FFTW to compute a two dimentional DFT of real data, I
readed the FFTW documentation an proced as follow:

/*** The fft.h file ***

#include <fftw3.h>
using namespace std;
class fft {

 public:  
  fft();
  fft(int rows, int cols);
  ~fft();
  fftw_complex * forward(double *);
  double * backward(fftw_complex *);
  fftw_complex * multiply(fftw_complex *, fftw_complex *);

 private: 
  fftw_plan plan_backward;
  fftw_plan plan_forward;
  int nx;
  int ny;
  int n;
  double *out;  
  fftw_complex * forward_result;
  fftw_complex* multiply_result;
};

/*** the fft.C file ***/

#include <iostream>
#include <rfftw.h>
#include "fft.h"

fft::fft(int rows, int cols) {
  nx = rows;
  ny = cols;
  n = nx*ny; 
  
}

fft::~fft() {
  fftw_destroy_plan (plan_backward );
  fftw_destroy_plan (plan_forward );
  delete[](out);
}
fftw_complex * fft::forward(double * in) {
  int nc = (n/2)+1;

  forward_result = (fftw_complex *) fftw_malloc ( sizeof ( fftw_complex )
* nx *(ny/2+1) );

  plan_forward = fftw_plan_dft_r2c_2d( nx, ny, in, forward_result,
FFTW_ESTIMATE);
  fftw_execute ( plan_forward );
  return(forward_result);
}
double * fft::backward(fftw_complex * in) {
  out = new double[nx*(ny/2+1)];
  plan_backward = fftw_plan_dft_c2r_2d(nx, ny, in, out, FFTW_ESTIMATE);
  fftw_execute ( plan_backward );
  /*  normalization */
  for(int i=0; i < nx*ny; i++)
    out[i] /= nx*ny;
  return(out);
}

fftw_complex * fft::multiply(fftw_complex * A, fftw_complex * B) {
  fftw_complex * R;
  double scale = 1.0 / (nx * ny);

  R = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) *nx * ny);
  for(int i= 0; i < nx; i++)
    for(int j= 0; j < ny/2+1; j++){
      int ij = i*(ny/2+1) + j;
      R[i][j].re = (A[ij].re * B[ij].re   
		    - A[ij].im * B[ij].im) * scale;  
      R[i][j].im = (A[ij].re * B[ij].im   
		    + A[ij].im * B[ij].re) * scale;   
    }     
  return(R);
}

/** The matrix.h file **/
#include <iostream>
#include <fstream>

//#include "fft.h"

using namespace std;

//class fft;

class Cmatrix
{
 private:
 protected:
  
  double * arr;
  int m_nRows;
  int nCols;
 
 public:
  
  Cmatrix();
  Cmatrix::Cmatrix(const Cmatrix& foo) {
    arr = foo.arr;
    m_nRows = foo.m_nRows;
    nCols = foo.nCols;
  }

  Cmatrix( int rows,int cols);
  Cmatrix(int rows, int cols, double * a);
  
  void setValue(int row,int col,double value);
  const double getValue( int row,  int col)const; 
  const int getRows() const ;
  const int getCols() const ;
  double * convolve(Cmatrix);
  Cmatrix * convolve(Cmatrix * pCmatrix_r);
  double * Cmatrix::getArray();
};
#endif

/** matrix.C file **/

#include "matrix.h"
#include <assert.h>

//constructors, destructors

Cmatrix::Cmatrix() {
  arr = NULL;
  m_nRows = 0;
  m_nCols = 0;

  arr = new double [rows*cols];
  assert(arr!=0);
  /* Initialization */
  for(int i=0;i<(rows*cols);i++)    {
    arr[i]=0.0;
  }
}

Cmatrix::~Cmatrix()
{
  if(arr!=NULL)    {
      delete[](arr);
      arr=NULL;
  }
}

void Cmatrix::setValue(int row, int col, double value)
{
  // check if access is legal 
  if( (row>m_nRows) || (col>m_nCols) || (row<=0) || (col<=0))    {
   cerr<<"ERROR: MATRIX OUT OF BOUND ACCESS, VALUE NOT SET !"<<endl;
     return;
  }
  int l = (row-1)*m_nCols + col-1;
  arr[l] = value;
}

const double Cmatrix::getValue( int row,  int col) const
{
 /* check if access is legal */
  const int ZERO=0;
  if( (row>m_nRows) || (col>m_nCols) || (row<=0) || (col<=0))
    {
       cerr<<"ERROR: MATRIX OUT OF BOUND ACCESS, ZERO RETURNED !"<<endl;
       return 0;
    }
  return arr[(row-1)*m_nCols+col-1];
}

/*inline*/ const int Cmatrix::getRows() const {

  return(m_nRows);
}
/*inline*/ const int Cmatrix::getCols() const {
  
  return(nCols);
}
/*inline*/ double * Cmatrix::getArray() {
  return(arr);
}

double * Cmatrix::convolve(Cmatrix Y) {

     fftw_real a[M][2*(N/2+1)], b[M][2*(N/2+1)], c[M][N];
     fftw_complex *A, *B, C[M][N/2+1];
     fftw_plan p, pinv, rhs;
     fftw_real scale = 1.0 / (M * N);
     int i, j;

  int M, N, nRowsZ, nColsZ;
  double * arrZ, * arrY;
  //h-- doublecomplex * pW;

  arrY = Y.getArray();
  M = Y.getRows();
  N = Y.getCols();
   
  nRowsZ = M+m_nRows-1;
  nColsZ = N+m_nCols-1;
  arrZ = new double [nRowsZ*nColsZ];

/* aliases for accessing complex transform outputs: */
     A = (fftw_complex*) &a[0][0];
     B = (fftw_complex*) &b[0][0];
     for (i = 0; i < M; ++i)
          for (j = 0; j < N; ++j) {
               a[i][j] = ... ;
               b[i][j] = ... ;
          }

  p = fftw_plan_dft_r2c_2d(M,N,&a[0][0],NULL,FFTW_ESTIMATE);
  pinv = fftw_plan_dft_r2c_2d(M,N,&b[0][0],NULL,FFTW_ESTIMATE);
  fftw_execute (p);
  fftw_execute (pinv);

     for (i = 0; i < M; ++i)
          for (j = 0; j < N/2+1; ++j) {
               int ij = i*(N/2+1) + j;
               C[i][j].re = (A[ij].re * B[ij].re
                             - A[ij].im * B[ij].im) * scale;
               C[i][j].im = (A[ij].re * B[ij].im
                             + A[ij].im * B[ij].re) * scale;
          }

  /* inverse transform to get c, the convolution of a and b;
        this has the side effect of overwriting C */

  rhs = fftw_plan_dft_c2r_2d(M,N,&C[0][0],&c[0][0],FFTW_ESTIMATE);
  fftw_execute (rhs);

  for(int a=0; a < nRowsY; a++)
    for(int b=0; b < nColsY; b++)
      arrZ[a*nColsY+b] = c[a][b];
  return(arrZ);
}

 Cmatrix * Cmatrix::convolve(Cmatrix * pCmatrix_rhs)
{
  
  fft F1(this->getRows(),this->getCols());
  // perform forward FFT transform 
  fftw_complex * f1 = F1.forward(this->getArray());
  fftw_complex * gx = F1.forward(pCmatrix_rhs->getArray()); 
  
  double * gy=F1.forward(m_pMatrixGreenY->getArray());
  // multiply in complex variable domain 
  fftw_complex * px  = F1.multiply(f1,gx);
  fftw_complex * py = F1.multiply(f1,gy);
  // delete to avoid memory leakage 
  delete []f1;
  delete []gx; 
  delete []gy;
  //perform backward FFT transform 
  double * rx  = F1.backward(px);
  double * ry  = F1.backward(py);
  // delete 
  delete []px;
  // result matrices generated from result double arrays 
Cmatrix * pCmatrix_result =new
Cmatrix(this->getRows(),this->getCols(),rx);
  return pCmatrix_result;
}


as you see my procedure convolve(Cmatrix Y) is imcomplete because I don't
understand the convolve with FFTW, and I am not sure about the procedure
multiply(fftw_complex * A, fftw_complex * B). pleased help me to correct
my faults.
Thanks

Hi!
little corrections

/*** The fft_l.h file ***

#include <fftw3.h>
using namespace std;
class fft {

 public:  
  fft();
  fft(int rows, int cols);
  ~fft();
  fftw_complex * forward(double *);
  double * backward(fftw_complex *);
  fftw_complex * multiply(fftw_complex *, fftw_complex *);

 private:
  fftw_plan plan_backward;
  fftw_plan plan_forward;
  int nx;
  int ny;
  int n;
  double *out;  
  fftw_complex * forward_result;
  fftw_complex* multiply_result;
};

/*** the fft_l.C file ***/

#include <iostream>
#include <rfftw.h>
#include "fft.h"

fft::fft(int rows, int cols) {
  nx = rows;
  ny = cols;
  n = nx*ny;
  
}

fft::~fft() {
  fftw_destroy_plan (plan_backward );
  fftw_destroy_plan (plan_forward );
  delete[](out);
}
fftw_complex * fft::forward(double * in) {
  int nc = (n/2)+1;

  forward_result = (fftw_complex *) fftw_malloc ( sizeof ( fftw_complex )
* nx *(ny/2+1) );

  plan_forward = fftw_plan_dft_r2c_2d( nx, ny, in, forward_result,
FFTW_ESTIMATE);
  fftw_execute ( plan_forward );
  return(forward_result);
}
double * fft::backward(fftw_complex * in) {
  out = new double[nx*(ny/2+1)];
  plan_backward = fftw_plan_dft_c2r_2d(nx, ny, in, out, FFTW_ESTIMATE);
  fftw_execute ( plan_backward );
  /*  normalization */
  for(int i=0; i < nx*ny; i++)
    out[i] /= nx*ny;
  return(out);
}

fftw_complex * fft::multiply(fftw_complex * A, fftw_complex * B) {
  fftw_complex * R;
  double scale = 1.0 / (nx * ny);

  R = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) *nx * ny);
  for(int i= 0; i < nx; i++)
    for(int j= 0; j < ny/2+1; j++){
      int ij = i*(ny/2+1) + j;
      R[i][j].re = (A[ij].re * B[ij].re   
            - A[ij].im * B[ij].im) * scale;  
      R[i][j].im = (A[ij].re * B[ij].im   
            + A[ij].im * B[ij].re) * scale;   
    }     
  return(R);
}

/** The matrix.h file **/
#include <iostream>
#include <fstream>
#include "fft_l.h"

using namespace std;
class fft;

class Cmatrix
{
 private:
 protected:
  
  double * arr;
  int m_nRows;
  int nCols;
 
 public:
  
  Cmatrix();
  Cmatrix::Cmatrix(const Cmatrix& foo) {
    arr = foo.arr;
    m_nRows = foo.m_nRows;
    m_nCols = foo.m_nCols;
  }

  Cmatrix( int rows,int cols);
  Cmatrix(int rows, int cols, double * a);
  
  void setValue(int row,int col,double value);
  const double getValue( int row,  int col)const;
  const int getRows() const ;
  const int getCols() const ;
  double * convolve(Cmatrix);
  Cmatrix * convolve(Cmatrix * pCmatrix_r);
  double * Cmatrix::getArray();
};
#endif

/** matrix.C file **/

#include "matrix.h"
#include <assert.h>

//constructors, destructors

Cmatrix::Cmatrix() {
  arr = NULL;
  m_nRows = 0;
  m_nCols = 0;

  arr = new double [rows*cols];
  assert(arr!=0);
  /* Initialization */
  for(int i=0;i<(rows*cols);i++)    {
    arr[i]=0.0;
  }
}

Cmatrix::~Cmatrix()
{
  if(arr!=NULL)    {
      delete[](arr);
      arr=NULL;
  }
}

void Cmatrix::setValue(int row, int col, double value)
{
  // check if access is legal
  if( (row>m_nRows) || (col>m_nCols) || (row<=0) || (col<=0))    {
   cerr<<"ERROR: MATRIX OUT OF BOUND ACCESS, VALUE NOT SET !"<<endl;
     return;
  }
  int l = (row-1)*m_nCols + col-1;
  arr[l] = value;
}

const double Cmatrix::getValue( int row,  int col) const
{
 /* check if access is legal */
  const int ZERO=0;
  if( (row>m_nRows) || (col>m_nCols) || (row<=0) || (col<=0))
    {
       cerr<<"ERROR: MATRIX OUT OF BOUND ACCESS, ZERO RETURNED !"<<endl;
       return 0;
    }
  return arr[(row-1)*m_nCols+col-1];
}

/*inline*/ const int Cmatrix::getRows() const {

  return(m_nRows);
}
/*inline*/ const int Cmatrix::getCols() const {
  
  return(nCols);
}
/*inline*/ double * Cmatrix::getArray() {
  return(arr);
}

double * Cmatrix::convolve(Cmatrix Y) {

     fftw_real a[M][2*(N/2+1)], b[M][2*(N/2+1)], c[M][N];
     fftw_complex *A, *B, C[M][N/2+1];
     fftw_plan p, pinv, rhs;
     fftw_real scale = 1.0 / (M * N);
     int i, j;

  int M, N, nRowsZ, nColsZ;
  double * arrZ, * arrY;
  //h-- doublecomplex * pW;

  arrY = Y.getArray();
  M = Y.getRows();
  N = Y.getCols();
   
  nRowsZ = M+m_nRows-1;
  nColsZ = N+m_nCols-1;
  arrZ = new double [nRowsZ*nColsZ];

/* aliases for accessing complex transform outputs: */
     A = (fftw_complex*) &a[0][0];
     B = (fftw_complex*) &b[0][0];
     for (i = 0; i < M; ++i)
          for (j = 0; j < N; ++j) {
               a[i][j] = ?? ;
               b[i][j] = ?? ;
          }

  p = fftw_plan_dft_r2c_2d(M,N,&a[0][0],NULL,FFTW_ESTIMATE);
  pinv = fftw_plan_dft_r2c_2d(M,N,&b[0][0],NULL,FFTW_ESTIMATE);
  fftw_execute (p);
  fftw_execute (pinv);

     for (i = 0; i < M; ++i)
          for (j = 0; j < N/2+1; ++j) {
               int ij = i*(N/2+1) + j;
               C[i][j].re = (A[ij].re * B[ij].re
                             - A[ij].im * B[ij].im) * scale;
               C[i][j].im = (A[ij].re * B[ij].im
                             + A[ij].im * B[ij].re) * scale;
          }

  /* inverse transform to get c, the convolution of a and b;
        this has the side effect of overwriting C */

  rhs = fftw_plan_dft_c2r_2d(M,N,&C[0][0],&c[0][0],FFTW_ESTIMATE);
  fftw_execute (rhs);

  for(int a=0; a < nRowsY; a++)
    for(int b=0; b < nColsY; b++)
      arrZ[a*nColsY+b] = c[a][b];
  return(arrZ);
}

 Cmatrix * Cmatrix::convolve(Cmatrix * pCmatrix_rhs)
{
  
  fft F1(this->getRows(),this->getCols());
  // perform forward FFT transform
  fftw_complex * f1 = F1.forward(this->getArray());
  fftw_complex * gx = F1.forward(pCmatrix_rhs->getArray());
  
  double * gy=F1.forward(m_pMatrixGreenY->getArray());
  // multiply in complex variable domain
  fftw_complex * px  = F1.multiply(f1,gx);
  fftw_complex * py = F1.multiply(f1,gy);
  // delete to avoid memory leakage
  fftw_free(f1);
  fftw_free(gx);
  fftw_freeg(y);
  //perform backward FFT transform
  double * rx  = F1.backward(px);
  double * ry  = F1.backward(py);
  // delete
  fftw_free(px);
  // result matrices generated from result double arrays
Cmatrix * pCmatrix_result =new
Cmatrix(this->getRows(),this->getCols(),rx);
  return pCmatrix_result;
}

The function multiply(fftw_complex * A, fftw_complex * B) to Perform
pointwise multiplication of A and B give me lot of faults like:
fft_l.C:118: request for member `re' in `*((((i * 2) + j) * 8) + R)',
which is
   of non-aggregate type `double'
fft_l.C:118: request for member `re' in `*(A + (+(ij * 16)))', which is
of
   non-aggregate type `double[2]'
fft_l.C:119: request for member `re' in `*(B + (+(ij * 16)))', which is
of
   non-aggregate type `double[2]'
fft_l.C:119: request for member `im' in `*(A + (+(ij * 16)))', which is
of
   non-aggregate type `double[2]'

the second point is to perform the two dimentional convolution of a Matrix
with the function convolve(Cmatrix Y), I have learned that is possible to
perfom it using the two dimentional fft routines to compute the
two-dimensional convolutions of the two-dimensional arrays A and B as
follow:
    1. Compute the two-dimensional FFT of A.
    2. Compute the two-dimensional FFT of B.
    3. Perform pointwise multiplication of A and B.
    4. Compute the inverse two-dimensional FFT of the previous result.
my problem is to define the two dimentional arrays A and B throught the
matrix and used the methode in the documantation of FFTW
http://www.fftw.org/fftw2_doc/fftw_2.html as I tried to do it in function
convolve(Cmatrix Y)
I'll like to know how can I compute the convolution of a matrix. 
Pleased help me.
Thanks