如何使用fftw Guru界面 [英] How to use fftw Guru interface
问题描述
我曾经使用fftw_plan_dft
进行多维傅里叶变换.
I used to use fftw_plan_dft
for multi-dimensional Fourier transformation.
fftw_plan fftw_plan_dft(int rank, const int *n, fftw_complex *in,
fftw_complex *out, int sign, unsigned flags);
现在我想将64位整数传递给fftw,看来我需要使用fftw guru接口.
Now I want to pass 64 bit integer to fftw, it looks like I need to use fftw guru interface.
fftw_plan fftw_plan_guru64_dft(
int rank, const fftw_iodim64 *dims,
int howmany_rank, const fftw_iodim64 *howmany_dims,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
但是我不明白什么是howmany_rank
和howmany_dims
的意思. fftw_plan_guru_dft
的手册说:
But I do not understand what is howmany_rank
and howmany_dims
mean. The manual of fftw_plan_guru_dft
says:
这两个函数分别针对交错格式和分割格式计划复杂数据,多维DFT.变换尺寸是由尺寸(howmany_rank,howmany_dims)的多维矢量(循环)给出的(等级,暗淡). dims和howmany_dims应该分别指向fftw_iodim数组,它们的长度分别为rank和howmany_rank.
These two functions plan a complex-data, multi-dimensional DFT for the interleaved and split format, respectively. Transform dimensions are given by (rank, dims) over a multi-dimensional vector (loop) of dimensions (howmany_rank, howmany_dims). dims and howmany_dims should point to fftw_iodim arrays of length rank and howmany_rank, respectively.
我确实知道什么是维度(howmany_rank,howmany_dims)的多维向量(循环)"的意思.您能给我一个例子或解释一下如何使用此大师界面吗?
I do know know what is "multi-dimensional vector (loop) of dimensions (howmany_rank, howmany_dims)" mean. Can you give me an example or explain how to use this guru interface?
推荐答案
If the sizes and strides of your mulitdimensional arrays are larger than 2^32, the 64 bit guru interface becomes useful.
创建复杂到复杂DTF的函数的原型为:
The prototype of the function creating complex to complex DTFs is:
fftw_plan fftw_plan_guru64_dft(
int rank, const fftw_iodim64 *dims,
int howmany_rank, const fftw_iodim64 *howmany_dims,
fftw_complex *in, fftw_complex *out,
int sign, unsigned flags);
其中:
-
rank
是要执行的FFTW变换的等级,即维数. -
dims
是大小为rank
的数组.对于每个尺寸i
,dims[i].n
是行的大小,dims[i].is
是输入数组的行之间的跨度,而dims[i].os
是输出数组的行之间的跨度.例如,如果数组在内存中是连续的,则
rank
is the rank of the FFTW transform to be performed, that is the number of dimensions.dims
is an array of sizerank
. For each dimensioni
,dims[i].n
is the size of the line,dims[i].is
is the stride between lines of the input array anddims[i].os
is the stride between lines of the output array. For instance, if the array is contiguous in memory, then the documentation of the guru interface suggests to use the recurrencedims[i].is = n[i+1] * dims[i+1].is
. The number of transforms to be performed and the offsets between the starting points are given byhowmany_rank
andhowmany_dims
.howmany_rank
specifies the number of transforms featuring particular offsets.howmany_dims
is an array of sizehowmany_rank
. For each transformi
,howmany_dims[i].n
is the number of transforms to be computed, each featuring offset between inputshowmany_dims[i].is
and offset between outputhowmany_dims[i].os
.
有关这些参数的更多详细信息,请参见关于FFTW3专家界面的困惑:3个同时进行的复杂FFT
More detail about these arguments are provided in Confusion about FFTW3 guru interface: 3 simultaneous complex FFTs
以下代码调用fftw_plan_guru64_dft()
,以便它执行与fftw_plan_dft_3d()
相同的操作.可以由gcc main.c -o main -lfftw3 -lm -Wall
:
The following code calls fftw_plan_guru64_dft()
so that it performs the same thing as fftw_plan_dft_3d()
. It can be compiled by gcc main.c -o main -lfftw3 -lm -Wall
:
#include<stdlib.h>
#include<complex.h>
#include<math.h>
#include<fftw3.h>
int main(void){
fftw_plan p;
unsigned long int N = 10;
unsigned long int M = 12;
unsigned long int P = 14;
fftw_complex *in=fftw_malloc(N*M*P*sizeof(fftw_complex));
if(in==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
fftw_complex *out=fftw_malloc(N*M*P*sizeof(fftw_complex));
if(out==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
unsigned int i,j,k;
int rank=3;
fftw_iodim64 *dims=malloc(rank*sizeof(fftw_iodim64));
if(dims==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
dims[0].n=N;
dims[0].is=P*M;
dims[0].os=P*M;
dims[1].n=M;
dims[1].is=P;
dims[1].os=P;
dims[2].n=P;
dims[2].is=1;
dims[2].os=1;
int howmany_rank=1;
fftw_iodim64 *howmany_dims=malloc(howmany_rank*sizeof(fftw_iodim64));
if(howmany_dims==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
howmany_dims[0].n=1;
howmany_dims[0].is=1;
howmany_dims[0].os=1;
printf("sizeof fftw complex %ld\n",sizeof(fftw_complex));
printf("sizeof fftw_iodim64 %ld\n",sizeof(fftw_iodim64));
printf("creating the plan\n");
p=fftw_plan_guru64_dft(rank, dims,howmany_rank, howmany_dims,in, out,FFTW_FORWARD, FFTW_ESTIMATE);
if (p==NULL){fprintf(stderr,"plan creation failed\n");exit(1);}
printf("created the plan\n");
for(i=0;i<N;i++){
for(j=0;j<M;j++){
for(k=0;k<P;k++){
//printf("ijk\n");
in[(i*M+j)*P+k]=30.+12.*sin(2*3.1415926535*i/((double)N))*sin(2*3.1415926535*j/((double)M))*sin(2*3.1415926535*k/((double)P))*I;
}
}
}
fftw_execute(p);
for (i = 0; i < N; i++){
for (j = 0; j < M; j++){
for (k = 0; k < P; k++){
printf("result: %d %d %d %g %gI\n", i,j,k, creal(out[(i*M+j)*P+k]), cimag(out[(i*M+j)*P+k]));
}
}
}
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
free(dims);
free(howmany_dims);
return(0);
}
例如,宗师接口可用于计算复杂3D电场的DFT.在网格的每个点处,电场都是大小为3的向量.因此,我可以将电场存储为4D数组,最后一个维指定向量的分量.最后,可以使用guru界面一次执行三个3D DFT:
For instance, the guru interface can be used to compute the DFT of a complex 3D electric field. At each point of the grid, the electric field is a vector of size 3. Hence, I can store the electric field as a 4D array, the last dimension specifying the component of the vector. Finally, the guru interface can be used to perform the three 3D DFTs at once:
#include<stdlib.h>
#include<complex.h>
#include<math.h>
#include<fftw3.h>
int main(void){
fftw_plan p;
unsigned long int N = 10;
unsigned long int M = 12;
unsigned long int P = 14;
unsigned long int DOF = 3;
fftw_complex *in=fftw_malloc(N*M*P*DOF*sizeof(fftw_complex));
if(in==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
fftw_complex *out=fftw_malloc(N*M*P*DOF*sizeof(fftw_complex));
if(out==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
unsigned int i,j,k;
int rank=3;
fftw_iodim64 *dims=malloc(rank*sizeof(fftw_iodim64));
if(dims==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
dims[0].n=N;
dims[0].is=P*M*DOF;
dims[0].os=P*M*DOF;
dims[1].n=M;
dims[1].is=P*DOF;
dims[1].os=P*DOF;
dims[2].n=P;
dims[2].is=DOF;
dims[2].os=DOF;
int howmany_rank=1;
fftw_iodim64 *howmany_dims=malloc(howmany_rank*sizeof(fftw_iodim64));
if(howmany_dims==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
howmany_dims[0].n=DOF;
howmany_dims[0].is=1;
howmany_dims[0].os=1;
printf("sizeof fftw complex %ld\n",sizeof(fftw_complex));
printf("sizeof fftw_iodim64 %ld\n",sizeof(fftw_iodim64));
printf("creating the plan\n");
p=fftw_plan_guru64_dft(rank, dims,howmany_rank, howmany_dims,in, out,FFTW_FORWARD, FFTW_ESTIMATE);
if (p==NULL){fprintf(stderr,"plan creation failed\n");exit(1);}
printf("created the plan\n");
for(i=0;i<N;i++){
for(j=0;j<M;j++){
for(k=0;k<P;k++){
//printf("ijk\n");
in[((i*M+j)*P+k)*DOF]=30.+12.*sin(2*3.1415926535*i/((double)N))*sin(2*3.1415926535*j/((double)M))*sin(2*3.1415926535*k/((double)P))*I;
in[((i*M+j)*P+k)*DOF+1]=42.0;
in[((i*M+j)*P+k)*DOF+2]=1.0;
}
}
}
fftw_execute(p);
for (i = 0; i < N; i++){
for (j = 0; j < M; j++){
for (k = 0; k < P; k++){
printf("result: %d %d %d || %g %gI | %g %gI | %g %gI\n", i,j,k, creal(out[((i*M+j)*P+k)*DOF]), cimag(out[((i*M+j)*P+k)*DOF]),creal(out[((i*M+j)*P+k)*DOF+1]), cimag(out[((i*M+j)*P+k)*DOF+1]),creal(out[((i*M+j)*P+k)*DOF+2]), cimag(out[((i*M+j)*P+k)*DOF+2]));
}
}
}
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
free(dims);
free(howmany_dims);
return(0);
}
这篇关于如何使用fftw Guru界面的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!