C++

Forum

 
2D FFT of Eigen matrices using FFTW library

JamieAl (234)
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 (784)
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 (234)
@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.