[SciPy-dev] Cleaning and fixing fft in scipy ?

[David Cournapeau]

Travis Oliphant wrote:
> Travis Oliphant wrote:
>> Could you be more clear about where the problem in your eyes lies?   
>> There are multiple sources for the fft (original fftpack files + 
>> interface files to other fft libraries if the user has those installed). 
>>
>> The ifdefs that I could find are just in the interface files to the 
>> "other" fft libraries that a user might have installed.  I don't care if 
>> that is redone so that the setup.py file just uses different sources as 
>> opposed to defining pre-processor variables, but you might check with 
>> Pearu since he is the author of those interfaces.
>>
>> But, I also don't really see the problem with the way it is done, now.
>>   
>
> I don't mean to discourage effort here.  I'm just hesitant to start 
> changing things without some kind of feedback from the original author 
> (Pearu in this case).
I agree, and that's exactly why I send this email.
>
> I just want to make sure he is aware of and agrees with changes that are 
> being made.   I know Pearu spent a fair bit of time creating the current 
> fft interfaces.  They may not be perfect but they are quite flexible and 
> allow people to use fft's from multiple libraries.  
I don't intend to change anything to the API. The reasoning is that 
scipy implementation of fft is suboptimal when using fftw3: it was 
really slow for some time, and I applied a quick fix, which is neither 
efficient or elegant (see ticket1 of scipy trac: 
http://projects.scipy.org/scipy/scipy/ticket/1). This patch took me a 
good time to get it right because of the fft sources organization. I 
understand that style is a matter of preferences more than anything 
else, and I am certainly not a reference, but I really think that with 
the kind of following code, it is really difficult to understand what's 
going on:

  int i;
  complex_double *ptr = inout;
#ifndef WITH_MKL
#if defined(WITH_FFTW)
  fftw_plan plan = NULL;
#endif
#endif
#if defined WITH_MKL
  DFTI_DESCRIPTOR_HANDLE desc_handle;
#else
  double* wsave = NULL;
#endif
#ifdef WITH_FFTWORK
  coef_dbl* coef = NULL;
#endif
#ifndef WITH_MKL
#ifdef WITH_DJBFFT
  int j;
  complex_double *ptrc = NULL;
  unsigned int *f = NULL;
#endif
#endif
#ifdef WITH_FFTWORK
  if (ispow2le2e30(n)) {
    i = get_cache_id_zfftwork(n);
    coef = caches_zfftwork[i].coef;
  } else
#endif
#ifndef WITH_MKL
#ifdef WITH_DJBFFT
  switch (n) {
  case 2:;case 4:;case 8:;case 16:;case 32:;case 64:;case 128:;case 256:;
  case 512:;case 1024:;case 2048:;case 4096:;case 8192:
    i = get_cache_id_zdjbfft(n);
    f = caches_zdjbfft[i].f;
    ptrc = (complex_double*)caches_zdjbfft[i].ptr;
  }
  if (f==0)
#endif
#endif
#ifdef WITH_MKL
  desc_handle = caches_zmklfft[get_cache_id_zmklfft(n)].desc_handle;
#elif defined WITH_FFTW
    plan = caches_zfftw[get_cache_id_zfftw(n,direction)].plan;
#else
    wsave = caches_zfftpack[get_cache_id_zfftpack(n)].wsave;
#endif
 
I think it is much more readable to have one file with (fftw3)

"""
complex_double  *ptr    = inout;
fftw_complex    *ptrm   = NULL;
fftw_plan       plan    = NULL;

int i;

plan = caches_zfftw[get_cache_id_zfftw(n, dir)].plan;
"""

One with (fftw case)

"""
int i;
complex_double *ptr = inout;
fftw_plan plan = NULL;
plan = caches_zfftw[get_cache_id_zfftw(n, direction)].plan;
"""

Etc... I guess the #ifdef solution was OK with one ot two 
implementation, but now, I fail to see how anyone can understand what 
the function does, which variables it is using, etc...

To improve the fftw3 implementation, I think we must have a better 
caching system, and I don't see how I could do that with the way the 
source is organized right now.
>
> Most of the fft behavior comes from those multiple libraries.   So, I 
> just want to be clear about which fft implementation and interface is 
> being modified and exactly what the problems are.
No interface would be changed; actually, I wouldn't be surprised if the 
compiled object code would be exactly the same with my suggestion for 
cleaning up (To be sure that I don't introduce bugs, I actually compare 
each new file with the preprocessed original file with gcc -E).

David