2D FFT of Eigen matrices using FFTW library

JamieAl (236)
I'm trying to use the FFTW library with Eigen. Since Eigen's unsupported module is slower. So, I'm struggling on implementing a very basic example. Here is what I have so far:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
 
#define _USE_MATH_DEFINES
#include <cmath>
#include<math.h>
#include<stdio.h>
#include "fftw3.h"
#include <cstring>
#include <sstream>
#include <string>
#include <sys/resource.h>
#include <iostream> 
#include <vector>
#include <fstream> 
#include <iomanip>
#include <numeric>
#include <assert.h>
#include <Eigen/Dense>
#include <unsupported/Eigen/FFT>
#include <Eigen/SparseCore> 
#include <Eigen/Sparse>

using namespace std;
using namespace Eigen;
int main(){

        static const int nx = 10;
	static const int ny = 10; 
	static const int nyk = ny/2 + 1;
	static const int ncomp = 2;
	double Lx = 2*EIGEN_PI; 
        double A = (2 * EIGEN_PI)/Lx;
	
        Matrix <double, ny+1, nx> eYY; 
        eYY.setZero();

        Matrix <double, ny+1, nx> eXX; 
        eXX.setZero();

	Eigen::Matrix< double, (ny+1), (ny)> u; 
	u.setZero(); 
	for(int i = 0; i< nx; i++){
		for(int j = 0; j< ny+1; j++){
                        eXX(j + (ny+1)*i) = (i)*2*EIGEN_PI/nx; 
			eYY(j + (ny+1)*i) = -1. * cos(((j) * EIGEN_PI )/ny);
			u(j + (ny+1)*i) = pow((eYY(j + (ny+1)*i)-0.5),2) * sin(A * eXX(j + (ny+1)*i)); //test function:u 

		}
	}
	
	fftw_complex *uhk; 
	uhk = (fftw_complex*) fftw_malloc(((ny*(ny+1)))* sizeof(fftw_complex)); 
	memset(uhk, 42, ((ny*(ny+1)))* sizeof(fftw_complex)); 

	
	Eigen::MatrixXcd uh;
	uh.setZero((ny+1), (ny));
	
	
	fftw_complex *uhkOut; 
	uhkOut = (fftw_complex*) fftw_malloc(((ny*(ny+1)))* sizeof(fftw_complex)); 
	memset(uhkOut, 42, ((ny*(ny+1)))* sizeof(fftw_complex)); 

	fftw_plan r2c; 
	r2c = fftw_plan_dft_r2c_2d(ny, (ny+1), &uhk[0], &uhkOut[0], FFTW_ESTIMATE);
	fftw_execute(r2c);
	fftw_destroy_plan(r2c);
	fftw_cleanup();
}

Basically, what I am trying to do is define an Eigen matrix u and try to take the 2D fft of it using my fftw library, where I had to copy memcpy the content over. But, I am getting the error:
1
2
3
4
5
6
7
8
9
 
Test.cpp: In function ‘int main()’:
Test.cpp:455:41: error: cannot convert ‘double (*)[2]’ to ‘double*’
  455 |  r2c = fftw_plan_dft_r2c_2d(ny, (ny+1), &uhk[0], &uhkOut[0], FFTW_ESTIMATE);
      |                                         ^~~~~~~
      |                                         |
      |                                         double (*)[2]
In file included from Test.cpp:8:
/usr/include/fftw3.h:457:1: note:   initializing argument 3 of ‘fftw_plan_s* fftw_plan_dft_r2c_2d(int, int, double*, double (*)[2], unsigned int)’
  457 | FFTW_DEFINE_API(FFTW_MANGLE_DOUBLE, double, fftw_complex)

Last edited on
kigar64551 (815)
It's because you are trying to pass a fftw_complex* pointer (which is a typedef for double[2]) as "in" parameter, but in fact a double* pointer is expected here! Only the "out" array is of type fftw_complex*:
https://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html

In other words, the "input" data has to be real numbers, but the "output" data will be complex numbers.

That's why it is called "r2c", by the way ;-)
Last edited on
JamieAl (236)
@kigar64551 hey you're so right, I guess writing code at 1 am is not a good idea. Thanks though for pointing it out. I see that my size if off now since I am getting the error:
1
2
 
corrupted size vs. prev_size
Aborted

But, I will take a second closer look.
Topic archived. No new replies allowed.