Commit 9d1991bf authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

remove support for obsolete legacy FFT libraries and point -DFFT_FFTW to FFTW3

parent 997142a4

doc/src/Section_start.txt

+4 −5
Original line number Diff line number Diff line
@@ -413,7 +413,7 @@ uses (for performing 1d FFTs) when running the particle-particle 
particle-mesh (PPPM) option for long-range Coulombics via the

"kspace_style"_kspace_style.html command.



LAMMPS supports various open-source or vendor-supplied FFT libraries

LAMMPS supports common open-source or vendor-supplied FFT libraries

for this purpose.  If you leave these 3 variables blank, LAMMPS will

use the open-source "KISS FFT library"_http://kissfft.sf.net, which is

included in the LAMMPS distribution.  This library is portable to all
@@ -423,10 +423,9 @@ package in your build, you can also leave the 3 variables blank. 


Otherwise, select which kinds of FFTs to use as part of the FFT_INC

setting by a switch of the form -DFFT_XXX.  Recommended values for XXX

are: MKL, SCSL, FFTW2, and FFTW3.  Legacy options are: INTEL, SGI,

ACML, and T3E.  For backward compatability, using -DFFT_FFTW will use

the FFTW2 library.  Using -DFFT_NONE will use the KISS library

described above.

are: MKL or FFTW3.  FFTW2 and NONE are supported as legacy options.

Selecting -DFFT_FFTW will use the FFTW3 library and -DFFT_NONE will

use the KISS library described above.



You may also need to set the FFT_INC, FFT_PATH, and FFT_LIB variables,

so the compiler and linker can find the needed FFT header and library

src/KSPACE/fft3d.cpp

+16 −431
Original line number Diff line number Diff line
@@ -72,20 +72,7 @@ void fft_3d(FFT_DATA *in, FFT_DATA *out, int flag, struct fft_plan_3d *plan) 


  // system specific constants



#if defined(FFT_SCSL)

  int isys = 0;

  FFT_PREC scalef = 1.0;

#elif defined(FFT_DEC)

  char c = 'C';

  char f = 'F';

  char b = 'B';

  int one = 1;

#elif defined(FFT_T3E)

  int isys = 0;

  double scalef = 1.0;

#elif defined(FFT_ACML)

  int info;

#elif defined(FFT_FFTW3)

#if defined(FFT_FFTW3)

  FFTW_API(plan) theplan;

#else

  // nothing to do for other FFTs.
@@ -109,35 +96,11 @@ void fft_3d(FFT_DATA *in, FFT_DATA *out, int flag, struct fft_plan_3d *plan) 
  total = plan->total1;

  length = plan->length1;



#if defined(FFT_SGI)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(flag,length,&data[offset],1,plan->coeff1);

#elif defined(FFT_SCSL)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(flag,length,scalef,&data[offset],&data[offset],plan->coeff1,

           plan->work1,&isys);

#elif defined(FFT_ACML)

  num=total/length;

  FFT_1D(&flag,&num,&length,data,plan->coeff1,&info);

#elif defined(FFT_INTEL)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(&data[offset],&length,&flag,plan->coeff1);

#elif defined(FFT_MKL)

#if defined(FFT_MKL)

  if (flag == -1)

    DftiComputeForward(plan->handle_fast,data);

  else

    DftiComputeBackward(plan->handle_fast,data);

#elif defined(FFT_DEC)

  if (flag == -1)

    for (offset = 0; offset < total; offset += length)

      FFT_1D(&c,&c,&f,&data[offset],&data[offset],&length,&one);

  else

    for (offset = 0; offset < total; offset += length)

      FFT_1D(&c,&c,&b,&data[offset],&data[offset],&length,&one);

#elif defined(FFT_T3E)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(&flag,&length,&scalef,&data[offset],&data[offset],plan->coeff1,

           plan->work1,&isys);

#elif defined(FFT_FFTW2)

  if (flag == -1)

    fftw(plan->plan_fast_forward,total/length,data,1,length,NULL,0,0);
@@ -172,35 +135,11 @@ void fft_3d(FFT_DATA *in, FFT_DATA *out, int flag, struct fft_plan_3d *plan) 
  total = plan->total2;

  length = plan->length2;



#if defined(FFT_SGI)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(flag,length,&data[offset],1,plan->coeff2);

#elif defined(FFT_SCSL)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(flag,length,scalef,&data[offset],&data[offset],plan->coeff2,

           plan->work2,&isys);

#elif defined(FFT_ACML)

  num=total/length;

  FFT_1D(&flag,&num,&length,data,plan->coeff2,&info);

#elif defined(FFT_INTEL)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(&data[offset],&length,&flag,plan->coeff2);

#elif defined(FFT_MKL)

#if defined(FFT_MKL)

  if (flag == -1)

    DftiComputeForward(plan->handle_mid,data);

  else

    DftiComputeBackward(plan->handle_mid,data);

#elif defined(FFT_DEC)

  if (flag == -1)

    for (offset = 0; offset < total; offset += length)

      FFT_1D(&c,&c,&f,&data[offset],&data[offset],&length,&one);

  else

    for (offset = 0; offset < total; offset += length)

      FFT_1D(&c,&c,&b,&data[offset],&data[offset],&length,&one);

#elif defined(FFT_T3E)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(&flag,&length,&scalef,&data[offset],&data[offset],plan->coeff2,

           plan->work2,&isys);

#elif defined(FFT_FFTW2)

  if (flag == -1)

    fftw(plan->plan_mid_forward,total/length,data,1,length,NULL,0,0);
@@ -235,35 +174,11 @@ void fft_3d(FFT_DATA *in, FFT_DATA *out, int flag, struct fft_plan_3d *plan) 
  total = plan->total3;

  length = plan->length3;



#if defined(FFT_SGI)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(flag,length,&data[offset],1,plan->coeff3);

#elif defined(FFT_SCSL)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(flag,length,scalef,&data[offset],&data[offset],plan->coeff3,

           plan->work3,&isys);

#elif defined(FFT_ACML)

  num=total/length;

  FFT_1D(&flag,&num,&length,data,plan->coeff3,&info);

#elif defined(FFT_INTEL)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(&data[offset],&length,&flag,plan->coeff3);

#elif defined(FFT_MKL)

#if defined(FFT_MKL)

  if (flag == -1)

    DftiComputeForward(plan->handle_slow,data);

  else

    DftiComputeBackward(plan->handle_slow,data);

#elif defined(FFT_DEC)

  if (flag == -1)

    for (offset = 0; offset < total; offset += length)

      FFT_1D(&c,&c,&f,&data[offset],&data[offset],&length,&one);

  else

    for (offset = 0; offset < total; offset += length)

      FFT_1D(&c,&c,&b,&data[offset],&data[offset],&length,&one);

#elif defined(FFT_T3E)

  for (offset = 0; offset < total; offset += length)

    FFT_1D(&flag,&length,&scalef,&data[offset],&data[offset],plan->coeff3,

           plan->work3,&isys);

#elif defined(FFT_FFTW2)

  if (flag == -1)

    fftw(plan->plan_slow_forward,total/length,data,1,length,NULL,0,0);
@@ -292,7 +207,6 @@ void fft_3d(FFT_DATA *in, FFT_DATA *out, int flag, struct fft_plan_3d *plan) 
             plan->post_plan);



  // scaling if required

#if !defined(FFT_T3E) && !defined(FFT_ACML)

  if (flag == 1 && plan->scaled) {

    norm = plan->norm;

    num = plan->normnum;
@@ -309,25 +223,6 @@ void fft_3d(FFT_DATA *in, FFT_DATA *out, int flag, struct fft_plan_3d *plan) 
#endif

    }

  }

#endif



#ifdef FFT_T3E

  if (flag == 1 && plan->scaled) {

    norm = plan->norm;

    num = plan->normnum;

    for (i = 0; i < num; i++) out[i] *= (norm,norm);

  }

#endif



#ifdef FFT_ACML

  norm = plan->norm;

  num = plan->normnum;

  for (i = 0; i < num; i++) {

    out[i].re *= norm;

    out[i].im *= norm;

  }

#endif



}



/* ----------------------------------------------------------------------
@@ -369,23 +264,6 @@ struct fft_plan_3d *fft_3d_create_plan( 
  int np1,np2,ip1,ip2;

  int list[50];



  // system specific variables



#ifdef FFT_SCSL

  FFT_DATA dummy_d[5];

  FFT_PREC dummy_p[5];

  int isign,isys;

  FFT_PREC scalef;

#endif

#ifdef FFT_INTEL

  FFT_DATA dummy;

#endif

#ifdef FFT_T3E

  FFT_DATA dummy[5];

  int isign,isys;

  double scalef;

#endif



  // query MPI info



  MPI_Comm_rank(comm,&me);
@@ -421,8 +299,7 @@ struct fft_plan_3d *fft_3d_create_plan( 
    first_klo = in_klo;

    first_khi = in_khi;

    plan->pre_plan = NULL;

  }

  else {

  } else {

    first_ilo = 0;

    first_ihi = nfast - 1;

    first_jlo = ip1*nmid/np1;
@@ -483,8 +360,7 @@ struct fft_plan_3d *fft_3d_create_plan( 
    third_jhi = out_jhi;

    third_klo = out_klo;

    third_khi = out_khi;

  }

  else {

  } else {

    third_ilo = ip1*nfast/np1;

    third_ihi = (ip1+1)*nfast/np1 - 1;

    third_jlo = ip2*nmid/np2;
@@ -601,136 +477,7 @@ struct fft_plan_3d *fft_3d_create_plan( 
  // system specific pre-computation of 1d FFT coeffs

  // and scaling normalization



#if defined(FFT_SGI)



  plan->coeff1 = (FFT_DATA *) malloc((nfast+15)*sizeof(FFT_DATA));

  plan->coeff2 = (FFT_DATA *) malloc((nmid+15)*sizeof(FFT_DATA));

  plan->coeff3 = (FFT_DATA *) malloc((nslow+15)*sizeof(FFT_DATA));



  if (plan->coeff1 == NULL || plan->coeff2 == NULL ||

      plan->coeff3 == NULL) return NULL;



  FFT_1D_INIT(nfast,plan->coeff1);

  FFT_1D_INIT(nmid,plan->coeff2);

  FFT_1D_INIT(nslow,plan->coeff3);



  if (scaled == 0)

    plan->scaled = 0;

  else {

    plan->scaled = 1;

    plan->norm = 1.0/(nfast*nmid*nslow);

    plan->normnum = (out_ihi-out_ilo+1) * (out_jhi-out_jlo+1) *

      (out_khi-out_klo+1);

  }



#elif defined(FFT_SCSL)



  plan->coeff1 = (FFT_PREC *) malloc((2*nfast+30)*sizeof(FFT_PREC));

  plan->coeff2 = (FFT_PREC *) malloc((2*nmid+30)*sizeof(FFT_PREC));

  plan->coeff3 = (FFT_PREC *) malloc((2*nslow+30)*sizeof(FFT_PREC));



  if (plan->coeff1 == NULL || plan->coeff2 == NULL ||

      plan->coeff3 == NULL) return NULL;



  plan->work1 = (FFT_PREC *) malloc((2*nfast)*sizeof(FFT_PREC));

  plan->work2 = (FFT_PREC *) malloc((2*nmid)*sizeof(FFT_PREC));

  plan->work3 = (FFT_PREC *) malloc((2*nslow)*sizeof(FFT_PREC));



  if (plan->work1 == NULL || plan->work2 == NULL ||

      plan->work3 == NULL) return NULL;



  isign = 0;

  scalef = 1.0;

  isys = 0;



  FFT_1D_INIT(isign,nfast,scalef,dummy_d,dummy_d,plan->coeff1,dummy_p,&isys);

  FFT_1D_INIT(isign,nmid,scalef,dummy_d,dummy_d,plan->coeff2,dummy_p,&isys);

  FFT_1D_INIT(isign,nslow,scalef,dummy_d,dummy_d,plan->coeff3,dummy_p,&isys);



  if (scaled == 0)

    plan->scaled = 0;

  else {

    plan->scaled = 1;

    plan->norm = 1.0/(nfast*nmid*nslow);

    plan->normnum = (out_ihi-out_ilo+1) * (out_jhi-out_jlo+1) *

      (out_khi-out_klo+1);

  }



#elif defined(FFT_ACML)



  plan->coeff1 = (FFT_DATA *) malloc((3*nfast+100)*sizeof(FFT_DATA));

  plan->coeff2 = (FFT_DATA *) malloc((3*nmid+100)*sizeof(FFT_DATA));

  plan->coeff3 = (FFT_DATA *) malloc((3*nslow+100)*sizeof(FFT_DATA));



  if (plan->coeff1 == NULL || plan->coeff2 == NULL ||

      plan->coeff3 == NULL) return NULL;



  int isign = 100;

  int isys = 1;

  int info = 0;

  FFT_DATA *dummy = NULL;



  FFT_1D(&isign,&isys,&nfast,dummy,plan->coeff1,&info);

  FFT_1D(&isign,&isys,&nmid,dummy,plan->coeff2,&info);

  FFT_1D(&isign,&isys,&nslow,dummy,plan->coeff3,&info);



  if (scaled == 0) {

    plan->scaled = 0;

    plan->norm = sqrt(nfast*nmid*nslow);

    plan->normnum = (out_ihi-out_ilo+1) * (out_jhi-out_jlo+1) *

      (out_khi-out_klo+1);

  } else {

    plan->scaled = 1;

    plan->norm = sqrt(nfast*nmid*nslow);

    plan->normnum = (out_ihi-out_ilo+1) * (out_jhi-out_jlo+1) *

      (out_khi-out_klo+1);

  }



#elif defined(FFT_INTEL)



  flag = 0;



  num = 0;

  factor(nfast,&num,list);

  for (i = 0; i < num; i++)

    if (list[i] != 2 && list[i] != 3 && list[i] != 5) flag = 1;

  num = 0;

  factor(nmid,&num,list);

  for (i = 0; i < num; i++)

    if (list[i] != 2 && list[i] != 3 && list[i] != 5) flag = 1;

  num = 0;

  factor(nslow,&num,list);

  for (i = 0; i < num; i++)

    if (list[i] != 2 && list[i] != 3 && list[i] != 5) flag = 1;



  MPI_Allreduce(&flag,&fftflag,1,MPI_INT,MPI_MAX,comm);

  if (fftflag) {

    if (me == 0) printf("ERROR: FFTs are not power of 2,3,5\n");

    return NULL;

  }



  plan->coeff1 = (FFT_DATA *) malloc((3*nfast/2+1)*sizeof(FFT_DATA));

  plan->coeff2 = (FFT_DATA *) malloc((3*nmid/2+1)*sizeof(FFT_DATA));

  plan->coeff3 = (FFT_DATA *) malloc((3*nslow/2+1)*sizeof(FFT_DATA));



  if (plan->coeff1 == NULL || plan->coeff2 == NULL ||

      plan->coeff3 == NULL) return NULL;



  flag = 0;

  FFT_1D_INIT(&dummy,&nfast,&flag,plan->coeff1);

  FFT_1D_INIT(&dummy,&nmid,&flag,plan->coeff2);

  FFT_1D_INIT(&dummy,&nslow,&flag,plan->coeff3);



  if (scaled == 0) {

    plan->scaled = 1;

    plan->norm = nfast*nmid*nslow;

    plan->normnum = (out_ihi-out_ilo+1) * (out_jhi-out_jlo+1) *

      (out_khi-out_klo+1);

  }

  else

    plan->scaled = 0;



#elif defined(FFT_MKL)

#if defined(FFT_MKL)

  DftiCreateDescriptor( &(plan->handle_fast), FFT_MKL_PREC, DFTI_COMPLEX, 1, (MKL_LONG)nfast);

  DftiSetValue(plan->handle_fast, DFTI_NUMBER_OF_TRANSFORMS, (MKL_LONG)plan->total1/nfast);

  DftiSetValue(plan->handle_fast, DFTI_PLACEMENT,DFTI_INPLACE);
@@ -761,50 +508,6 @@ struct fft_plan_3d *fft_3d_create_plan( 
      (out_khi-out_klo+1);

  }



#elif defined(FFT_DEC)



  if (scaled == 0) {

    plan->scaled = 1;

    plan->norm = nfast*nmid*nslow;

    plan->normnum = (out_ihi-out_ilo+1) * (out_jhi-out_jlo+1) *

      (out_khi-out_klo+1);

  }

  else

    plan->scaled = 0;



#elif defined(FFT_T3E)



  plan->coeff1 = (double *) malloc((12*nfast)*sizeof(double));

  plan->coeff2 = (double *) malloc((12*nmid)*sizeof(double));

  plan->coeff3 = (double *) malloc((12*nslow)*sizeof(double));



  if (plan->coeff1 == NULL || plan->coeff2 == NULL ||

      plan->coeff3 == NULL) return NULL;



  plan->work1 = (double *) malloc((8*nfast)*sizeof(double));

  plan->work2 = (double *) malloc((8*nmid)*sizeof(double));

  plan->work3 = (double *) malloc((8*nslow)*sizeof(double));



  if (plan->work1 == NULL || plan->work2 == NULL ||

      plan->work3 == NULL) return NULL;



  isign = 0;

  scalef = 1.0;

  isys = 0;



  FFT_1D_INIT(&isign,&nfast,&scalef,dummy,dummy,plan->coeff1,dummy,&isys);

  FFT_1D_INIT(&isign,&nmid,&scalef,dummy,dummy,plan->coeff2,dummy,&isys);

  FFT_1D_INIT(&isign,&nslow,&scalef,dummy,dummy,plan->coeff3,dummy,&isys);



  if (scaled == 0)

    plan->scaled = 0;

  else {

    plan->scaled = 1;

    plan->norm = 1.0/(nfast*nmid*nslow);

    plan->normnum = (out_ihi-out_ilo+1) * (out_jhi-out_jlo+1) *

      (out_khi-out_klo+1);

  }



#elif defined(FFT_FFTW2)



  plan->plan_fast_forward =
@@ -941,36 +644,10 @@ void fft_3d_destroy_plan(struct fft_plan_3d *plan) 
  if (plan->copy) free(plan->copy);

  if (plan->scratch) free(plan->scratch);



#if defined(FFT_SGI)

  free(plan->coeff1);

  free(plan->coeff2);

  free(plan->coeff3);

#elif defined(FFT_SCSL)

  free(plan->coeff1);

  free(plan->coeff2);

  free(plan->coeff3);

  free(plan->work1);

  free(plan->work2);

  free(plan->work3);

#elif defined(FFT_ACML)

  free(plan->coeff1);

  free(plan->coeff2);

  free(plan->coeff3);

#elif defined(FFT_INTEL)

  free(plan->coeff1);

  free(plan->coeff2);

  free(plan->coeff3);

#elif defined(FFT_MKL)

#if defined(FFT_MKL)

  DftiFreeDescriptor(&(plan->handle_fast));

  DftiFreeDescriptor(&(plan->handle_mid));

  DftiFreeDescriptor(&(plan->handle_slow));

#elif defined(FFT_T3E)

  free(plan->coeff1);

  free(plan->coeff2);

  free(plan->coeff3);

  free(plan->work1);

  free(plan->work2);

  free(plan->work3);

#elif defined(FFT_FFTW2)

  if (plan->plan_slow_forward != plan->plan_fast_forward &&

      plan->plan_slow_forward != plan->plan_mid_forward) {
@@ -1015,38 +692,31 @@ void factor(int n, int *num, int *list) 
{

  if (n == 1) {

    return;

  }

  else if (n % 2 == 0) {

  } else if (n % 2 == 0) {

    *list = 2;

    (*num)++;

    factor(n/2,num,list+1);

  }

  else if (n % 3 == 0) {

  } else if (n % 3 == 0) {

    *list = 3;

    (*num)++;

    factor(n/3,num,list+1);

  }

  else if (n % 5 == 0) {

  } else if (n % 5 == 0) {

    *list = 5;

    (*num)++;

    factor(n/5,num,list+1);

  }

  else if (n % 7 == 0) {

  } else if (n % 7 == 0) {

    *list = 7;

    (*num)++;

    factor(n/7,num,list+1);

  }

  else if (n % 11 == 0) {

  } else if (n % 11 == 0) {

    *list = 11;

    (*num)++;

    factor(n/11,num,list+1);

  }

  else if (n % 13 == 0) {

  } else if (n % 13 == 0) {

    *list = 13;

    (*num)++;

    factor(n/13,num,list+1);

  }

  else {

  } else {

    *list = n;

    (*num)++;

    return;
@@ -1089,23 +759,6 @@ void fft_1d_only(FFT_DATA *data, int nsize, int flag, struct fft_plan_3d *plan) 
  int i,total,length,offset,num;

  FFT_SCALAR norm, *data_ptr;



  // system specific constants



#ifdef FFT_SCSL

  int isys = 0;

  FFT_PREC scalef = 1.0;

#endif

#ifdef FFT_DEC

  char c = 'C';

  char f = 'F';

  char b = 'B';

  int one = 1;

#endif

#ifdef FFT_T3E

  int isys = 0;

  double scalef = 1.0;

#endif



  // total = size of data needed in each dim

  // length = length of 1d FFT in each dim

  // total/length = # of 1d FFTs in each dim
@@ -1131,39 +784,7 @@ void fft_1d_only(FFT_DATA *data, int nsize, int flag, struct fft_plan_3d *plan) 
  // perform 1d FFTs in each of 3 dimensions

  // data is just an array of 0.0



#ifdef FFT_SGI

  for (offset = 0; offset < total1; offset += length1)

    FFT_1D(flag,length1,&data[offset],1,plan->coeff1);

  for (offset = 0; offset < total2; offset += length2)

    FFT_1D(flag,length2,&data[offset],1,plan->coeff2);

  for (offset = 0; offset < total3; offset += length3)

    FFT_1D(flag,length3,&data[offset],1,plan->coeff3);

#elif defined(FFT_SCSL)

  for (offset = 0; offset < total1; offset += length1)

    FFT_1D(flag,length1,scalef,&data[offset],&data[offset],plan->coeff1,

           plan->work1,&isys);

  for (offset = 0; offset < total2; offset += length2)

    FFT_1D(flag,length2,scalef,&data[offset],&data[offset],plan->coeff2,

           plan->work2,&isys);

  for (offset = 0; offset < total3; offset += length3)

    FFT_1D(flag,length3,scalef,&data[offset],&data[offset],plan->coeff3,

           plan->work3,&isys);

#elif defined(FFT_ACML)

  int info=0;

  num=total1/length1;

  FFT_1D(&flag,&num,&length1,data,plan->coeff1,&info);

  num=total2/length2;

  FFT_1D(&flag,&num,&length2,data,plan->coeff2,&info);

  num=total3/length3;

  FFT_1D(&flag,&num,&length3,data,plan->coeff3,&info);

#elif defined(FFT_INTEL)

  for (offset = 0; offset < total1; offset += length1)

    FFT_1D(&data[offset],&length1,&flag,plan->coeff1);

  for (offset = 0; offset < total2; offset += length2)

    FFT_1D(&data[offset],&length2,&flag,plan->coeff2);

  for (offset = 0; offset < total3; offset += length3)

    FFT_1D(&data[offset],&length3,&flag,plan->coeff3);

#elif defined(FFT_MKL)

#if defined(FFT_MKL)

  if (flag == -1) {

    DftiComputeForward(plan->handle_fast,data);

    DftiComputeForward(plan->handle_mid,data);
@@ -1173,32 +794,6 @@ void fft_1d_only(FFT_DATA *data, int nsize, int flag, struct fft_plan_3d *plan) 
    DftiComputeBackward(plan->handle_mid,data);

    DftiComputeBackward(plan->handle_slow,data);

  }

#elif defined(FFT_DEC)

  if (flag == -1) {

    for (offset = 0; offset < total1; offset += length1)

      FFT_1D(&c,&c,&f,&data[offset],&data[offset],&length1,&one);

    for (offset = 0; offset < total2; offset += length2)

      FFT_1D(&c,&c,&f,&data[offset],&data[offset],&length2,&one);

    for (offset = 0; offset < total3; offset += length3)

      FFT_1D(&c,&c,&f,&data[offset],&data[offset],&length3,&one);

  } else {

    for (offset = 0; offset < total1; offset += length1)

      FFT_1D(&c,&c,&b,&data[offset],&data[offset],&length1,&one);

    for (offset = 0; offset < total2; offset += length2)

      FFT_1D(&c,&c,&b,&data[offset],&data[offset],&length2,&one);

    for (offset = 0; offset < total3; offset += length3)

      FFT_1D(&c,&c,&b,&data[offset],&data[offset],&length3,&one);

  }

#elif defined(FFT_T3E)

  for (offset = 0; offset < total1; offset += length1)

    FFT_1D(&flag,&length1,&scalef,&data[offset],&data[offset],plan->coeff1,

           plan->work1,&isys);

  for (offset = 0; offset < total2; offset += length2)

    FFT_1D(&flag,&length2,&scalef,&data[offset],&data[offset],plan->coeff2,

           plan->work2,&isys);

  for (offset = 0; offset < total3; offset += length3)

    FFT_1D(&flag,&length3,&scalef,&data[offset],&data[offset],plan->coeff3,

           plan->work3,&isys);

#elif defined(FFT_FFTW2)

  if (flag == -1) {

    fftw(plan->plan_fast_forward,total1/length1,data,1,0,NULL,0,0);
@@ -1247,7 +842,6 @@ void fft_1d_only(FFT_DATA *data, int nsize, int flag, struct fft_plan_3d *plan) 
  // scaling if required

  // limit num to size of data



#ifndef FFT_T3E

  if (flag == 1 && plan->scaled) {

    norm = plan->norm;

    num = MIN(plan->normnum,nsize);
@@ -1264,13 +858,4 @@ void fft_1d_only(FFT_DATA *data, int nsize, int flag, struct fft_plan_3d *plan) 
#endif

    }

  }

#endif



#ifdef FFT_T3E

  if (flag == 1 && plan->scaled) {

    norm = plan->norm;

    num = MIN(plan->normnum,nsize);

    for (i = 0; i < num; i++) data[i] *= (norm,norm);

  }

#endif

}

src/KSPACE/fft3d.h

+4 −154
Original line number Diff line number Diff line
@@ -27,7 +27,7 @@ typedef double FFT_SCALAR; 


// set default fftw library. switch to FFT_FFTW3 when convenient.

#ifdef FFT_FFTW

#define FFT_FFTW2

#define FFT_FFTW3

#endif



// -------------------------------------------------------------------------
@@ -36,73 +36,11 @@ typedef double FFT_SCALAR; 


#if FFT_PRECISION == 1



#if defined(FFT_SGI)

#include "fft.h"

typedef complex FFT_DATA;

#define FFT_1D cfft1d

#define FFT_1D_INIT cfft1di

extern "C" {

  int cfft1d(int, int, FFT_DATA *, int, FFT_DATA *);

  FFT_DATA *cfft1di(int, FFT_DATA *);

}



#elif defined(FFT_SCSL)

#include <scsl_fft.h>

typedef scsl_complex FFT_DATA;

typedef float FFT_PREC;

#define FFT_1D ccfft

#define FFT_1D_INIT ccfft

extern "C" {

  int ccfft(int, int, FFT_PREC, FFT_DATA *, FFT_DATA *,

                      FFT_PREC *, FFT_PREC *, int *);

}



#elif defined(FFT_ACML)

typedef struct {

  float re;

  float im;

} FFT_DATA;

#define FFT_1D cfft1m_

extern "C" {

  void cfft1m_(int *, int *, int *, FFT_DATA *, FFT_DATA *, int *);

}



#elif defined(FFT_INTEL)

typedef struct {

  float re;

  float im;

} FFT_DATA;

#define FFT_1D cfft1d_

#define FFT_1D_INIT cfft1d_

extern "C" {

  void cfft1d_(FFT_DATA *, int *, int *, FFT_DATA *);

}



#elif defined(FFT_MKL)

#if defined(FFT_MKL)

#include "mkl_dfti.h"

typedef float _Complex FFT_DATA;

#define FFT_MKL_PREC DFTI_SINGLE



#elif defined(FFT_DEC)

typedef struct {

  float re;

  float im;

} FFT_DATA;

#define FFT_1D cfft_

extern "C" {

  void cfft_(char *, char *, char *, FFT_DATA *, FFT_DATA *, int *, int *);

}



#elif defined(FFT_T3E)

#include <complex.h>

typedef complex single FFT_DATA;

#define FFT_1D GGFFT

#define FFT_1D_INIT GGFFT

extern "C" {

  void GGFFT(int *, int *, double *, FFT_DATA *, FFT_DATA *,

             double *, double *, int *);

}



#elif defined(FFT_FFTW2)

#if defined(FFTW_SIZE)

#include "sfftw.h"
@@ -138,73 +76,11 @@ typedef struct kiss_fft_state* kiss_fft_cfg; 


#elif FFT_PRECISION == 2



#if defined(FFT_SGI)

#include "fft.h"

typedef zomplex FFT_DATA;

#define FFT_1D zfft1d

#define FFT_1D_INIT zfft1di

extern "C" {

  int zfft1d(int, int, FFT_DATA *, int, FFT_DATA *);

  FFT_DATA *zfft1di(int, FFT_DATA *);

}



#elif defined(FFT_SCSL)

#include <scsl_fft.h>

typedef scsl_zomplex FFT_DATA;

typedef double FFT_PREC;

#define FFT_1D zzfft

#define FFT_1D_INIT zzfft

extern "C" {

  int zzfft(int, int, FFT_PREC, FFT_DATA *, FFT_DATA *,

                      FFT_PREC *, FFT_PREC *, int *);

}



#elif defined(FFT_ACML)

typedef struct {

  double re;

  double im;

} FFT_DATA;

#define FFT_1D zfft1m_

extern "C" {

  void zfft1m_(int *, int *, int *, FFT_DATA *, FFT_DATA *, int *);

}



#elif defined(FFT_INTEL)

typedef struct {

  double re;

  double im;

} FFT_DATA;

#define FFT_1D zfft1d_

#define FFT_1D_INIT zfft1d_

extern "C" {

  void zfft1d_(FFT_DATA *, int *, int *, FFT_DATA *);

}



#elif defined(FFT_MKL)

#if defined(FFT_MKL)

#include "mkl_dfti.h"

typedef double _Complex FFT_DATA;

#define FFT_MKL_PREC DFTI_DOUBLE



#elif defined(FFT_DEC)

typedef struct {

  double re;

  double im;

} FFT_DATA;

#define FFT_1D zfft_

extern "C" {

  void zfft_(char *, char *, char *, FFT_DATA *, FFT_DATA *, int *, int *);

}



#elif defined(FFT_T3E)

#include <complex.h>

typedef complex double FFT_DATA;

#define FFT_1D CCFFT

#define FFT_1D_INIT CCFFT

extern "C" {

  void CCFFT(int *, int *, double *, FFT_DATA *, FFT_DATA *,

             double *, double *, int *);

}



#elif defined(FFT_FFTW2)

#if defined(FFTW_SIZE)

#include "dfftw.h"
@@ -258,36 +134,10 @@ struct fft_plan_3d { 
  double norm;                      // normalization factor for rescaling



                                    // system specific 1d FFT info

#if defined(FFT_SGI)

  FFT_DATA *coeff1;

  FFT_DATA *coeff2;

  FFT_DATA *coeff3;

#elif defined(FFT_SCSL)

  FFT_PREC *coeff1;

  FFT_PREC *coeff2;

  FFT_PREC *coeff3;

  FFT_PREC *work1;

  FFT_PREC *work2;

  FFT_PREC *work3;

#elif defined(FFT_ACML)

  FFT_DATA *coeff1;

  FFT_DATA *coeff2;

  FFT_DATA *coeff3;

#elif defined(FFT_INTEL)

  FFT_DATA *coeff1;

  FFT_DATA *coeff2;

  FFT_DATA *coeff3;

#elif defined(FFT_MKL)

#if defined(FFT_MKL)

  DFTI_DESCRIPTOR *handle_fast;

  DFTI_DESCRIPTOR *handle_mid;

  DFTI_DESCRIPTOR *handle_slow;

#elif defined(FFT_T3E)

  double *coeff1;

  double *coeff2;

  double *coeff3;

  double *work1;

  double *work2;

  double *work3;

#elif defined(FFT_FFTW2)

  fftw_plan plan_fast_forward;

  fftw_plan plan_fast_backward;

src/MAKE/MACHINES/Makefile.bgl

+1 −1
Original line number Diff line number Diff line
@@ -56,7 +56,7 @@ MPI_LIB = -lmpich.rts 
# PATH = path for FFT library

# LIB = name of FFT library



FFT_INC =       -DFFT_FFTW

FFT_INC =       -DFFT_FFTW2

FFT_PATH = 

FFT_LIB =	-lfftw

src/MAKE/MACHINES/Makefile.bgq

+1 −1
Original line number Diff line number Diff line
@@ -8,7 +8,7 @@ SHELL = /bin/bash 
# do not edit this section

# select which compiler by editing Makefile.bgq.details



include ../MAKE/Makefile.bgq.details

include ../MAKE/MACHINES/Makefile.bgq.details



include Makefile.package.settings

include Makefile.package

CRA Git | Maintained and supported by SUSTech CRA and CCSE