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