Using fftw with column major square matrices (armadillo library)

Question

I find the armadillo C++ library very convenient for matrix computations. How does one perform a two dimensional FFT on armadillo matrices using the FFTW library?

I understand that armadillo matrix class stores data in a column major order. How do I pass this to FFTW? The fftw 3.3.3 documentation says

If you have an array stored in column-major order and wish to transform it using FFTW, it is quite easy to do. When creating the plan, simply pass the dimensions of the array to the planner in reverse order. For example, if your array is a rank three N x M x L matrix in column-major order, you should pass the dimensions of the array as if it were an L x M x N matrix (which it is, from the perspective of FFTW)

I cannot fully understand what this means, given the syntax for creating a plan is as follows.

fftw_plan fftw_plan_dft_2d(int n0, int n1,
                            fftw_complex *in, fftw_complex *out,
                            int sign, unsigned flags);

Could someone explain this?

Answer 1

It just means that you can try both

  fftw_plan plan=fftw_plan_dft_2d(4, 2,(double(*)[2])&AAA(0,0), (double(*)[2])&BBB(0,0), FFTW_FORWARD, FFTW_ESTIMATE);

and

 fftw_plan plan=fftw_plan_dft_2d(2, 4,(double(*)[2])&AAA(0,0), (double(*)[2])&BBB(0,0), FFTW_FORWARD, FFTW_ESTIMATE);

And keep the correct order.

Normally, in and out would be allocated by fftw_malloc() to keep aligned memory. Little tests show that it is already the case with cx_mat matrix from Armadillo (no horrible segmentation fault or wrong values...).

Here is the test code (and the correct order...) :

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

using namespace arma;
using namespace std;


int main(int argc, char** argv)
{
    cout << "Armadillo version: " << arma_version::as_string() << endl;


    cx_mat AAA = eye<cx_mat>(2,4);
    AAA(0,0)=0;
    AAA(0,1)=1;
    AAA(0,2)=2;
    AAA(0,3)=3;

    AAA(1,0)=0;
    AAA(1,1)=1;
    AAA(1,2)=2;
    AAA(1,3)=3;


    cx_mat BBB = eye<cx_mat>(2,4);

    fftw_plan plan=fftw_plan_dft_2d(4, 2,(double(*)[2])&AAA(0,0), (double(*)[2])&BBB(0,0), FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(plan);

    BBB.print("BBB:");
    return 0;
}

Compile it with g++ -O2 -o example1 example1.cpp -larmadillo -llapack -lblas -lfftw3 !

Bye,

Francis