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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
| #include <fftw3.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
int main(void)
{
fftw_plan pb, pf;
double *in, *back;
fftw_complex *out;
int nx, ny;
int i, j, ij;
double xm, ym;
int n[1];
int inembed[1];
int onembed[1];
int rank;
int howmany;
int istride, ostride;
int idist, odist;
double x,y;
FILE *fpin, *fpback;
nx = 128;
ny = 64;
xm = 20.0;
ym = 20.0;
in = fftw_malloc((nx+1)*(ny+1)*sizeof(*in));
out = fftw_malloc((nx/2+1)*(ny+1)*sizeof(*out));
back = fftw_malloc((nx+1)*(ny+1)*sizeof(*back));
rank = 1;
howmany = ny+1;
istride = 1;
ostride = 1;
idist = nx+1;
odist = nx/2+1;
n[0] = nx;
inembed[0] = nx+1;
onembed[0] = nx/2+1;
fpin = fopen("in.txt" , "w");
for (i = 0; i < nx+1; i++)
{
for (j = 0; j < ny+1; j++)
{
ij = i + j*(nx+1);
x = i*xm/(double)nx;
y = j*ym/(double)ny;
in[ij] = cos(2.0*M_PI*x/xm)*sin(2.0*M_PI*y/ym);
fprintf(fpin, "%lg %lg %lg\n", x, y, in[ij]);
}
}
fclose(fpin);
pf = fftw_plan_many_dft_r2c(rank, n, howmany, in,
inembed, istride, idist, out,
onembed, ostride, odist, FFTW_MEASURE);
fftw_execute(pf);
rank = 1;
howmany = ny+1;
istride = 1;
ostride = 1;
idist = nx/2+1;
odist = nx+1;
n[0] = nx/2+1;
inembed[0] = nx/2+1;
onembed[0] = nx+1;
pb = fftw_plan_many_dft_c2r(rank, n, howmany,
out, inembed, istride,
idist, back, onembed,
ostride, odist, FFTW_MEASURE);
fftw_execute(pb);
fpback = fopen("back.txt" , "w");
for (i = 0; i < nx+1; i++)
{
for (j = 0; j < ny+1; j++)
{
ij = i + j*(nx+1);
x = i*xm/(double)nx;
y = j*ym/(double)ny;
fprintf(fpback, "%lg %lg %lg\n", x, y, back[ij]);
}
}
fclose(fpback);
fftw_destroy_plan(pf);
fftw_destroy_plan(pb);
free(in);
free(out);
free(back);
return 0;
} |
