CODEWITHSUNDEEP
c++fftwarmadillo
user1285419
user1285419

Reputation: 2225

Interfacing Armadillo with FFTW

I am trying to use FFTW in the matrix of Armadillo package. I need to perform N independent FFT on a 2D matrix. Following the manual of FFTW and other resource online, I have the following code:

#include <armadillo>
#include <iostream>
#include "fftw3.h"

using namespace std;
using namespace arma;

fftw_plan fplan;

int main(int argc, char *argv[])
{
  cx_mat psi;
  int NCOL=2, NROW=4;
  int frank=1;
  int howmany=NCOL;
  int n[]={NROW};         
  int idist=NROW, odist=NROW; 
  int istride=1, ostride=1;
  int *inembed=n, *onembed=n;

  psi.resize(NROW, NCOL);
  psi(0,0)=cx_double(1,1); psi(0,1)=cx_double(1,1);
  psi(1,0)=cx_double(2,1); psi(1,1)=cx_double(1,2);
  psi(2,0)=cx_double(3,1); psi(2,1)=cx_double(1,3);
  psi(3,0)=cx_double(4,1); psi(3,1)=cx_double(1,4);
  cout << psi << endl; // the output is correct at here

  fftw_complex* in = reinterpret_cast<fftw_complex*>(psi.memptr());

  fplan = fftw_plan_many_dft(frank, n, howmany,
                             in, inembed, istride, idist,
                             in, onembed, ostride, odist,
                             FFTW_FORWARD, FFTW_MEASURE);

  cout << endl << "after fftw plan: " << endl << psi << endl; // all zeros

  // fftw_execute(fplan); // output zeros again if execute
  cx_double* A_mem=psi.memptr();
  cout << endl << "by reference: " << endl << *A_mem << " " << *(A_mem+1) << endl;

  return 0;
}

The code compiles without any error. But if I run the code, after fftw_plan_many_dft, it outputs all zeroes instead. If I remove fftw_plan_many_dft, the output is correct. I don't even execute the fftw plan so why it clean my data by just setting up the plan?

Upvotes: 1

Views: 1380

Answers (1)

Paul R
Paul R

Reputation: 213220

With the FFTW_MEASURE flag the buffers are actually used to tune the FFT algorithm. You should create the plan first, and then initialise the input data.

Note that FFTW_ESTIMATE does not exhibit this behaviour, so you could also just change this flag in your current code to fix the problem, but then your code would probably run slower.

Upvotes: 3

Related Questions