[Bug-gsl] Memoryproblem when using fourierroutines

[Anders Wallander]

Hi!
I Sent a email about my problems to the discuss-forum and Brian Gough told me 
that you had tested the fft-routines for memoryleakage with several tools.
But after some trial-and-error-testing around I concluded that my problem lies 
in the two routines below (fourier_transform_2Ddouble_complex and 
fourier_inverse_transform_2Dcomplex_magnitude_double).

Problem:
I'm using them for a cross-correlation routine when I read in packed pictures 
one at a time into gsl_matrixes. I create a new matrix, turned 180 degrees and 
and then fouriertransform both, multiply all elements and inverse-transform 
back. This gives me a cross-correlation map. The programs works fine, except 
that these routines eates up more and more of the memory for every picture 
correlated.
For some reason gdb won't load my home-made library so I used the debugger 
valgrind which showed that something fishy was going on when allocating 
wavetable and workspace. As you see I'm very precise with free:ing all 
allocated memory before leaving the routines.
I'm using gsl-1.3.

Would appreciate if you either found the bug or point out where I have missed 
something.

Regards,
Anders Wallander

-------------------------------------------------------   
gsl_matrix_complex *
fourier_transform_2Ddouble_complex(gsl_matrix const *matrix_d1){
  size_t      row, col;
  size_t      mrow  = matrix_d1->size1;
  size_t      mcol  = matrix_d1->size2;
  double      *cpdata;                     /* array in format 
gsl_complex_packed_array */
  gsl_complex ctmp;
  gsl_fft_complex_wavetable *wave_table;   /* Pointer to wavetable in format 
complex  */
  gsl_fft_complex_workspace *work_space;   /* Pointer to workspace in format 
complex  */
  gsl_matrix_complex        *matrix_c2;
  
  if( (mcol > FOURIER_WORKSPACE/2) || (mrow > FOURIER_WORKSPACE/2)){ 
    printf("The workspace set is only enough for %d x %d matrix \n", 
FOURIER_WORKSPACE/2, FOURIER_WORKSPACE/2);
    return NULL;
  }
  if ((cpdata = malloc(2 * mcol * sizeof(double))) == NULL) {
    printf("Couldn't allocate memory for cpdata.\n");
    return NULL;
  }  
  
  matrix_c2 =  gsl_matrix_complex_alloc(mrow, mcol);
  for(row=0; row < mrow; row++){
    for(col=0; col < mcol; col++){
      cpdata[2*col]   = gsl_matrix_get(matrix_d1, row, col); 
      cpdata[2*col+1] = 0.0;                                 
    }
    wave_table = gsl_fft_complex_wavetable_alloc(mcol);      
    work_space = gsl_fft_complex_workspace_alloc(mcol); 
    gsl_fft_complex_forward(cpdata, 1, mcol, wave_table, work_space);
    for(col=0; col < mcol; col++){
      GSL_SET_COMPLEX(&ctmp, cpdata[2*col], cpdata[2*col+1]);
      gsl_matrix_complex_set(matrix_c2, row, col, ctmp); 
    }
  }
  if (mrow != mcol){                                  /* security check */
    free(cpdata);
    if ((cpdata = malloc(2*mrow* sizeof(double))) == NULL) {
      printf("Couldn't allocate memory for cpdata.\n");
      return NULL;
    }  
    gsl_fft_complex_wavetable_free(wave_table);
    gsl_fft_complex_workspace_free(work_space);
    wave_table = gsl_fft_complex_wavetable_alloc(mrow);      
    work_space = gsl_fft_complex_workspace_alloc(mrow);
    
  }
  for(col=0; col < mcol; col++){
    for(row=0; row < mrow; row++){
      cpdata[2*row]   = GSL_REAL(gsl_matrix_complex_get(matrix_c2, row, col)); 
      cpdata[2*row+1] = GSL_IMAG(gsl_matrix_complex_get(matrix_c2, row, col));
    }
    gsl_fft_complex_forward(cpdata, 1, mrow, wave_table, work_space);
    for(row=0; row < mrow; row++){
      GSL_SET_COMPLEX(&ctmp, cpdata[2*row], cpdata[2*row+1]);
      gsl_matrix_complex_set(matrix_c2, row, col, ctmp); 
    }
  } 
  gsl_fft_complex_wavetable_free(wave_table);
  gsl_fft_complex_workspace_free(work_space);
  return matrix_c2;
}

gsl_matrix *
fourier_inverse_transform_2Dcomplex_magnitude_double(gsl_matrix_complex const 
*matrix_c1){
  size_t      row, col;
  size_t      mrow = matrix_c1->size1;
  size_t      mcol = matrix_c1->size2;
  double      *cpdata;
  double      dtmp;
  gsl_complex ctmp;
  gsl_fft_complex_wavetable *wave_table;
  gsl_fft_complex_workspace *work_space;
  gsl_matrix_complex        *matrix_c2;
  gsl_matrix                *matrix_d;
  
  if( (mcol > FOURIER_WORKSPACE/2) || (mrow > FOURIER_WORKSPACE/2)){ 
    printf("The workspace set is only enough for %d x %d matrix \n", 
FOURIER_WORKSPACE/2, FOURIER_WORKSPACE/2);
    return NULL;
  }
  if ((cpdata = malloc(2 * mcol * sizeof(double))) == NULL) {
    printf("Couldn't allocate memory for cpdata.\n");
    return NULL;
  }
  
  matrix_c2 = gsl_matrix_complex_alloc(mrow, mcol);
  matrix_d  = gsl_matrix_alloc(mrow, mcol);
  
  for(row=0; row < mrow; row++){
    for(col=0; col < mcol; col++){
      cpdata[2*col]   =  GSL_REAL(gsl_matrix_complex_get(matrix_c1, row, col)); 
      cpdata[2*col+1] =  GSL_IMAG(gsl_matrix_complex_get(matrix_c1, row, col));
    }
    wave_table = gsl_fft_complex_wavetable_alloc(mcol);      
    work_space = gsl_fft_complex_workspace_alloc(mcol);
    gsl_fft_complex_inverse(cpdata, 1, mcol, wave_table, work_space);
    for(col=0; col < mcol; col++){
      GSL_SET_COMPLEX(&ctmp, cpdata[2*col], cpdata[2*col+1]);
      gsl_matrix_complex_set(matrix_c2, row, col, ctmp);
    }
  }
  
  if (mrow != mcol){                                  /* security check */
    free(cpdata);
    if ((cpdata = malloc(2 * mrow * sizeof(double))) == NULL) {
      printf("Couldn't allocate memory for cpdata.\n");
      return NULL;
    } 
    gsl_fft_complex_wavetable_free(wave_table);
    gsl_fft_complex_workspace_free(work_space);
    wave_table = gsl_fft_complex_wavetable_alloc(mrow); 
    work_space = gsl_fft_complex_workspace_alloc(mrow);
  }
  
  for(col=0; col < mcol; col++){
    for(row=0; row < mrow; row++){
      cpdata[2*row]   =  GSL_REAL(gsl_matrix_complex_get(matrix_c2, row, col)); 
      cpdata[2*row+1] =  GSL_IMAG(gsl_matrix_complex_get(matrix_c2, row, col));
    }
    gsl_fft_complex_inverse(cpdata, 1, mrow, wave_table, work_space);
    for(row=0; row < mrow; row++){
      GSL_SET_COMPLEX(&ctmp, cpdata[2*row], cpdata[2*row+1]);
      gsl_matrix_complex_set(matrix_c2, row, col, ctmp); 
    } 
  }
  gsl_fft_complex_wavetable_free(wave_table);
  gsl_fft_complex_workspace_free(work_space);
  
  for(row=0; row < mrow; row++){
    for(col=0; col < mcol; col++){
      ctmp = gsl_matrix_complex_get(matrix_c2, row, col);
      dtmp = sqrt(GSL_REAL(ctmp)*GSL_REAL(ctmp)+GSL_IMAG(ctmp)*GSL_IMAG(ctmp));
      gsl_matrix_set(matrix_d, row, col, dtmp);
    }
  }
  gsl_matrix_complex_free(matrix_c2);
  return matrix_d; 
}