[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Re: Memoryproblem when using fourierroutines
From: |
Anders Wallander |
Subject: |
[Bug-gsl] Re: Memoryproblem when using fourierroutines |
Date: |
Tue, 22 Jul 2003 15:45:38 +0200 |
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;
> }