最近在写经典密度泛函方面的程序代码， 为了提高程序的运行速度， 卷积部分用FFT来计算。

花了半天时间研究了一下FFTW的用法，FFTW在对长度为N的实数序列做变换时， 结果只保留长度为(N/2+1)的复数数组， 其余的值可以根据Y_k = Y^*_N-k来获得，这在编程的时候要小心处理。

另外，FFTW在做逆变换时，没有除以周期长度，这点也要注意。

自己写的一个小练习，上传上来作为记录。

#include<complex.h>
#include<fftw3.h>
#include<stdio.h>
#include<stdlib.h>
#define N 80
int main()
{
int i;
double *g,*f,*c;
fftw_complex *f_out,*g_out,*c_out;
fftw_plan p1, p2, p3;
f = (double*)malloc(sizeof(double) * N);
g = (double*)malloc(sizeof(double) * N);
c = (double*)malloc(sizeof(double) * N);
f_out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
g_out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
c_out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

if((f==NULL)||(g==NULL)||(c==NULL))
{
printf("Error:insufficient available memoryn");
}
else
{
for(i=0; i<N; i++)
{
f[i] = 1;
g[i] = 2;
}
}

p1 = fftw_plan_dft_r2c_1d(N, f, f_out,FFTW_ESTIMATE);
p2 = fftw_plan_dft_r2c_1d(N, g, g_out,FFTW_ESTIMATE);
p3 = fftw_plan_dft_c2r_1d(N, c_out, c,FFTW_ESTIMATE);

fftw_execute(p1); /* repeat as needed */
fftw_execute(p2); /* repeat as needed */

for(i=0; i<N; i++){
if(i<= (N/2 + 1)){
c_out[i] = f_out[i] * g_out[i];
}else{
c_out[i] = conj(f_out[N-i]) * conj(g_out[N-i]);
}
}


fftw_execute(p3); /* repeat as needed */
fftw_destroy_plan(p1);
fftw_destroy_plan(p2);
fftw_destroy_plan(p3);
fftw_cleanup();
for(i=0;i<N;i++)/*OUTPUT*/
{
printf("%f,%f,%fn",f[i],g[i],c[i]/N);
}

free(f);
free(g);
free(c);
fftw_free(f_out);
fftw_free(g_out);
fftw_free(c_out);
return 0;
    }