fftw Guruインターフェイスの使い方

Question

私は多次元フーリエ変換にfftw_plan_dftを使用しました。

fftw_plan fftw_plan_dft(int rank, const int *n, fftw_complex *in, 
         fftw_complex *out, int sign, unsigned flags); 

は今、私はFFTWの第一人者・インタフェースを使用する必要があるように見える、FFTWする64ビット整数を渡したいです。

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のマニュアルは言う：

これら二つの機能は、それぞれ、インターリーブ及び分割形式のための複雑なデータ、マルチ次元DFTを計画します。変形寸法は、寸法（howmany_rank、h​​owmany_dims）の多次元ベクトル（ループ）にわたって（ランク、ディム）によって与えられる。 dimm、howmany_dimsは、それぞれlength rankとhowmany_rankのfftw_iodim配列を指すはずです。

"次元の多次元ベクトル（howmany_rank、h​​owmany_dims）"が意味することは分かっています。あなたは私に例を与えたり、このグル・インターフェースの使い方を説明できますか？

Answer 1

多次元配列のサイズとストライドが2^32より大きい場合、64 bit 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); 

：

複合DTFSに複合体を作成関数のプロトタイプである

• rankがFFTW実行される変換の順位であり、それは次元数です。
• dimsは、サイズrankの配列です。各次元についてi、dims[i].nは行のサイズ、dims[i].isは入力配列の行間のストライド、dims[i].osは出力配列の行間のストライドです。たとえば、配列がメモリ内で連続している場合は、the documentation of the guru interfaceは再帰dims[i].is = n[i+1] * dims[i+1].isを使用することを示します。 実行する変換の数と開始点間のオフセットは、howmany_rankとhowmany_dimsで与えられます。
• howmany_rankは、特定のオフセットを特徴とする変換の数を指定します。
• howmany_dimsは、サイズhowmany_rankの配列です。 iの各変換では、howmany_dims[i].nは計算する変換の数であり、それぞれの入力間のオフセットはhowmany_dims[i].isで、出力間のオフセットはhowmany_dims[i].osです。

次のコードは、fftw_plan_dft_3d()と同じことを実行するようにfftw_plan_guru64_dft()を呼び出します。 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配列として保存できます。最後の次元はベクトルの成分を指定します。最後に、グルインターフェースは一度に3次元DFTを実行するために使用することができます。

#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=3; 
    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)*3]=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)*3+1]=42.0; 
       in[((i*M+j)*P+k)*3+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)*3]), cimag(out[((i*M+j)*P+k)*3]),creal(out[((i*M+j)*P+k)*3+1]), cimag(out[((i*M+j)*P+k)*3+1]),creal(out[((i*M+j)*P+k)*3+2]), cimag(out[((i*M+j)*P+k)*3+2])); 
      } 
     } 
    } 


    fftw_destroy_plan(p); 
    fftw_free(in); 
    fftw_free(out); 

    free(dims); 
    free(howmany_dims); 

    return(0); 
}