[Bug-gsl] Re: Memoryproblem when using fourierroutines

[Anders Wallander]

Additional info:
* 2 x 735 Mhz Pentium 3, Linux Red Hat 8.0
* GSL-1.3, it was compiled by our sysop.
* gcc version 3.2 20020903 (Red Hat Linux 8.0 3.2-7)
/Anders

On Tue, 22 Jul 2003 15:14:32 +0200
Anders Wallander <address@hidden> wrote:

> 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; 
> }