How do I use fftw_plan_many_dft on a transposed array of data?

Question

I have a 2D array of data stored in column-major (Fortran-style) format, and I'd like to take the FFT of each row. I would like to avoid transposing the array (it is not square). For example, my array

fftw_complex* data = new fftw_complex[21*256];

contains entries [r0_val0, r1_val0,..., r20_val0, r0_val1,...,r20_val255].

I can use fftw_plan_many_dft to make a plan to solve each of the 21 FFTs in-place in the data array if it is row-major, e.g. [r0_val0, r0_val1,..., r0_val255, r1_val0,...,r20_val255]:

int main() {
  int N = 256;
  int howmany = 21;
  fftw_complex* data = new fftw_complex[N*howmany];
  fftw_plan p;

  // this plan is OK
  p = fftw_plan_many_dft(1,&N,howmany,data,NULL,1,N,data,NULL,1,N,FFTW_FORWARD,FFTW_MEASURE);
  // do stuff...

  return 0;
}

According to the documentation (section 4.4.1 of the FFTW manual), the signature for the function is

fftw_plan fftw_plan_many_dft(int rank, const int *n, int howmany,
                              fftw_complex *in, const int *inembed,
                              int istride, int idist,
                              fftw_complex *out, const int *onembed,
                              int ostride, int odist,
                              int sign, unsigned flags);

and I should be able to use the stride and dist parameters to set the indexing. From what I can understand from the documentation, the entries in the array to be transformed are indexed as in + j*istride + k*idist where j=0..n-1 and k=0..howmany-1. (My arrays are 1D and there are howmany of them). However, the following code results in a seg. fault (edit: the stride length is wrong, see update below):

int main() {
  int N = 256;
  int howmany = 21;
  fftw_complex* data = new fftw_complex[N*howmany];
  fftw_plan p;

  // this call results in a seg. fault.
  p = fftw_plan_many_dft(1,&N,howmany,data,NULL,N,1,data,NULL,N,1,FFTW_FORWARD,FFTW_MEASURE);

  return 0;
}

Update:

I made an error choosing the stride length. The correct call is (the correct stride length is howmany, not N):

int main() {
  int N = 256;
  int howmany = 21;
  fftw_complex* data = new fftw_complex[N*howmany];
  fftw_plan p;

  // OK
  p = fftw_plan_many_dft(1,&N,howmany,data,NULL,howmany,1,data,NULL,howmany,1,FFTW_FORWARD,FFTW_MEASURE);
  // do stuff  

  return 0;
}

Larder

Answer 1

The function works as documented. I made an error with the stride length, which should actually be howmany in this case. I have updated the question to reflect this.

I find the documentation for FFTW is somewhat difficult to comprehend without examples (I could just be illiterate...), so I'm posting a more detailed example below, comparing the usual use of fftw_plan_dft_1d with fftw_plan_many_dft. To recap, in the case of howmany arrays with length N that are stored in a contiguous block of memory referenced as in, the array elements j for each transform k are

*(in + j*istride + k*idist)

The following two pieces of code are equivalent. In the first, the conversion from some 2D array is done explicitly, and in the second the fftw_plan_many_dft call is used to do everything in-place.

Explicit Copy

int N, howmany;
// ...
fftw_complex* data = (fftw_complex*) fftw_malloc(N*howmany*sizeof(fftw_complex));
// ... load data with howmany arrays of length N 
int istride, idist;
// ... if data is column-major, set istride=howmany, idist=1
//    if data is row-major, set istride=1, idist=N

fftw_complex* in = (fftw_complex*) fftw_malloc(N*sizeof(fftw_complex));
fftw_complex* out = (fftw_complex*) fftw_malloc(N*sizeof(fftw_complex));
fftw_plan p = fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_MEASURE);
for (int k = 0; k < howmany; ++k) {
  for (int j = 0; j < N; ++j) {
    in[j] = data[j*istride + k*idist];
  }
  fftw_execute(p);
  // do something with out
}

Plan Many

int N, howmany;
// ...
fftw_complex* data = (fftw_complex*) fftw_malloc(N*howmany*sizeof(fftw_complex));
// ... load data with howmany arrays of length N 
int istride, idist;
// ... if data is column-major, set istride=howmany, idist=1
//    if data is row-major, set istride=1, idist=N

fftw_plan p = fftw_plan_many_dft(1,&N,howmany,data,NULL,howmany,1,data,NULL,howmany,1,FFTW_FORWARD,FFTW_MEASURE);
fftw_execute(p);

Larder

I am not sure this is correct @Dylan. If you change your data array in then you have to run fftw_plan_dft_1d a second time using your new input array. Although I can't be sure it will give you an error as written, I can tell you the documentation explicitly warns against it (page 4 of the pdf manual to version 3.2.2). – Breland

@Breland I think you may have misunderstood the manual. Are you referring to the sentence that begins “If you want to transform a different array of the same size”? That is in reference to using a different array, rather than to changing the contents of the existing one. – Fernandina