src/mdlib/gmx_fft_fftw3.c

   1 /*
   2  *
   3  * Gromacs 4.0                         Copyright (c) 1991-2003
   4  * David van der Spoel, Erik Lindahl, University of Groningen.
   5  *
   6  * This program is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU General Public License
   8  * as published by the Free Software Foundation; either version 2
   9  * of the License, or (at your option) any later version.
  10  *
  11  * To help us fund GROMACS development, we humbly ask that you cite
  12  * the research papers on the package. Check out http://www.gromacs.org
  13  *
  14  * And Hey:
  15  * Gnomes, ROck Monsters And Chili Sauce
  16  */
  17 #ifdef HAVE_CONFIG_H
  18 #include <config.h>
  19 #endif
  20
  21 #ifdef GMX_FFT_FFTW3
  22
  23 #include <errno.h>
  24 #include <stdlib.h>
  25
  26 #include <fftw3.h>
  27
  28 #include "gmx_fft.h"
  29 #include "gmx_fatal.h"
  30
  31
  32 #ifdef GMX_DOUBLE
  33 #define FFTWPREFIX(name) fftw_ ## name
  34 #else
  35 #define FFTWPREFIX(name) fftwf_ ## name
  36 #endif
  37
  38
  39
  40 #ifdef GMX_THREAD_MPI
  41 /* none of the fftw3 calls, except execute(), are thread-safe, so
  42    we need to serialize them with this mutex. */
  43 static tMPI_Thread_mutex_t big_fftw_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
  44 static gmx_bool gmx_fft_threads_initialized=FALSE;
  45 #define FFTW_LOCK tMPI_Thread_mutex_lock(&big_fftw_mutex)
  46 #define FFTW_UNLOCK tMPI_Thread_mutex_unlock(&big_fftw_mutex)
  47 #else /* GMX_THREAD_MPI */
  48 #define FFTW_LOCK
  49 #define FFTW_UNLOCK
  50 #endif /* GMX_THREAD_MPI */
  51
  52 /* We assume here that aligned memory starts at multiple of 16 bytes and unaligned memory starts at multiple of 8 bytes. The later is guranteed for all malloc implementation.
  53    Consequesences:
  54    - It is not allowed to use these FFT plans from memory which doesn't have a starting address as a multiple of 8 bytes.
  55      This is OK as long as the memory directly comes from malloc and is not some subarray within alloated memory.
  56    - This has to be fixed if any future architecute requires memory to be aligned to multiples of 32 bytes.
  57 */
  58
  59 struct gmx_fft
  60 {
  61     /* Three alternatives (unaligned/aligned, out-of-place/in-place, forward/backward)
  62      * results in 8 different FFTW plans. Keep track of them with 3 array indices:
  63      * first index:   0=unaligned, 1=aligned
  64      * second index:  0=out-of-place, 1=in-place
  65      * third index:   0=backward, 1=forward
  66      */
  67     FFTWPREFIX(plan)         plan[2][2][2];
  68     /* Catch user mistakes */
  69     int                      real_transform;
  70     int                      ndim;
  71 };
  72
  73
  74
  75 int
  76 gmx_fft_init_1d(gmx_fft_t *        pfft,
  77                 int                nx,
  78                 gmx_fft_flag       flags)
  79 {
  80     return gmx_fft_init_many_1d(pfft,nx,1,flags);
  81 }
  82
  83
  84 int
  85 gmx_fft_init_many_1d(gmx_fft_t *        pfft,
  86                      int                nx,
  87                      int                howmany,
  88                      gmx_fft_flag       flags)
  89 {
  90     gmx_fft_t              fft;
  91     FFTWPREFIX(complex)   *p1,*p2,*up1,*up2;
  92     size_t                 pc;
  93     int                    i,j,k;
  94     int                    fftw_flags;
  95
  96 #ifdef GMX_DISABLE_FFTW_MEASURE
  97     flags |= GMX_FFT_FLAG_CONSERVATIVE;
  98 #endif
  99
 100     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 101
 102     if(pfft==NULL)
 103     {
 104         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 105         return EINVAL;
 106     }
 107     *pfft = NULL;
 108
 109     FFTW_LOCK;
 110     if( (fft = (gmx_fft_t)FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
 111     {
 112         FFTW_UNLOCK;
 113         return ENOMEM;
 114     }
 115
 116     /* allocate aligned, and extra memory to make it unaligned */
 117     p1  = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
 118     if(p1==NULL)
 119     {
 120         FFTWPREFIX(free)(fft);
 121         FFTW_UNLOCK;
 122         return ENOMEM;
 123     }
 124
 125     p2  = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
 126     if(p2==NULL)
 127     {
 128         FFTWPREFIX(free)(p1);
 129         FFTWPREFIX(free)(fft);
 130         FFTW_UNLOCK;
 131         return ENOMEM;
 132     }
 133
 134     /* make unaligned pointers.
 135      * In double precision the actual complex datatype will be 16 bytes,
 136      * so go to a char pointer and force an offset of 8 bytes instead.
 137      */
 138     pc = (size_t)p1;
 139     pc += 8;
 140     up1 = (FFTWPREFIX(complex) *)pc;
 141
 142     pc = (size_t)p2;
 143     pc += 8;
 144     up2 = (FFTWPREFIX(complex) *)pc;
 145
 146     /*                            int rank, const int *n, int howmany,
 147                                   fftw_complex *in, const int *inembed,
 148                                   int istride, int idist,
 149                                   fftw_complex *out, const int *onembed,
 150                                   int ostride, int odist,
 151                                   int sign, unsigned flags */
 152     fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
 153     fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
 154     fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
 155     fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
 156     fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
 157     fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
 158     fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
 159     fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
 160
 161     for(i=0;i<2;i++)
 162     {
 163         for(j=0;j<2;j++)
 164         {
 165             for(k=0;k<2;k++)
 166             {
 167                 if(fft->plan[i][j][k] == NULL)
 168                 {
 169                     gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
 170                     FFTW_UNLOCK;
 171                     gmx_fft_destroy(fft);
 172                     FFTW_LOCK;
 173                     FFTWPREFIX(free)(p1);
 174                     FFTWPREFIX(free)(p2);
 175                     FFTW_UNLOCK;
 176                     return -1;
 177                 }
 178             }
 179         }
 180     }
 181
 182     FFTWPREFIX(free)(p1);
 183     FFTWPREFIX(free)(p2);
 184
 185     fft->real_transform = 0;
 186     fft->ndim           = 1;
 187
 188     *pfft = fft;
 189     FFTW_UNLOCK;
 190     return 0;
 191 }
 192
 193
 194 int
 195 gmx_fft_init_1d_real(gmx_fft_t *        pfft,
 196                      int                nx,
 197                      gmx_fft_flag       flags)
 198 {
 199     return gmx_fft_init_many_1d_real(pfft, nx, 1, flags);
 200 }
 201
 202 int
 203 gmx_fft_init_many_1d_real(gmx_fft_t *        pfft,
 204                      int                nx,
 205                      int                howmany,
 206                      gmx_fft_flag       flags)
 207 {
 208     gmx_fft_t              fft;
 209     real            *p1,*p2,*up1,*up2;
 210     size_t                pc;
 211     int                   i,j,k;
 212     int                    fftw_flags;
 213
 214 #ifdef GMX_DISABLE_FFTW_MEASURE
 215     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 216 #endif
 217
 218     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 219
 220     if(pfft==NULL)
 221     {
 222         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 223         return EINVAL;
 224     }
 225     *pfft = NULL;
 226
 227     FFTW_LOCK;
 228     if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
 229     {
 230         FFTW_UNLOCK;
 231         return ENOMEM;
 232     }
 233
 234     /* allocate aligned, and extra memory to make it unaligned */
 235     p1  = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
 236     if(p1==NULL)
 237     {
 238         FFTWPREFIX(free)(fft);
 239         FFTW_UNLOCK;
 240         return ENOMEM;
 241     }
 242
 243     p2  = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
 244     if(p2==NULL)
 245     {
 246         FFTWPREFIX(free)(p1);
 247         FFTWPREFIX(free)(fft);
 248         FFTW_UNLOCK;
 249         return ENOMEM;
 250     }
 251
 252     /* make unaligned pointers.
 253      * In double precision the actual complex datatype will be 16 bytes,
 254      * so go to a char pointer and force an offset of 8 bytes instead.
 255      */
 256     pc = (size_t)p1;
 257     pc += 8;
 258     up1 = (real *)pc;
 259
 260     pc = (size_t)p2;
 261     pc += 8;
 262     up2 = (real *)pc;
 263
 264     /*                                int rank, const int *n, int howmany,
 265                                       double *in, const int *inembed,
 266                                       int istride, int idist,
 267                                       fftw_complex *out, const int *onembed,
 268                                       int ostride, int odist,
 269                                       unsigned flag    */
 270     fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft_r2c)(1,&nx,howmany,up1,0,1,(nx/2+1)*2,(FFTWPREFIX(complex) *)up2,0,1,(nx/2+1),fftw_flags);
 271     fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft_r2c)(1,&nx,howmany,up1,0,1,(nx/2+1)*2,(FFTWPREFIX(complex) *)up1,0,1,(nx/2+1),fftw_flags);
 272     fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft_r2c)(1,&nx,howmany, p1,0,1,(nx/2+1)*2,(FFTWPREFIX(complex) *)p2 ,0,1,(nx/2+1),fftw_flags);
 273     fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft_r2c)(1,&nx,howmany, p1,0,1,(nx/2+1)*2,(FFTWPREFIX(complex) *)p1 ,0,1,(nx/2+1),fftw_flags);
 274
 275     fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft_c2r)(1,&nx,howmany,(FFTWPREFIX(complex) *)up1,0,1,(nx/2+1),up2,0,1,(nx/2+1)*2,fftw_flags);
 276     fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft_c2r)(1,&nx,howmany,(FFTWPREFIX(complex) *)up1,0,1,(nx/2+1),up1,0,1,(nx/2+1)*2,fftw_flags);
 277     fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft_c2r)(1,&nx,howmany,(FFTWPREFIX(complex) *) p1,0,1,(nx/2+1), p2,0,1,(nx/2+1)*2,fftw_flags);
 278     fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft_c2r)(1,&nx,howmany,(FFTWPREFIX(complex) *) p1,0,1,(nx/2+1), p1,0,1,(nx/2+1)*2,fftw_flags);
 279
 280     for(i=0;i<2;i++)
 281     {
 282         for(j=0;j<2;j++)
 283         {
 284             for(k=0;k<2;k++)
 285             {
 286                 if(fft->plan[i][j][k] == NULL)
 287                 {
 288                     gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
 289                     FFTW_UNLOCK;
 290                     gmx_fft_destroy(fft);
 291                     FFTW_LOCK;
 292                     FFTWPREFIX(free)(p1);
 293                     FFTWPREFIX(free)(p2);
 294                     FFTW_UNLOCK;
 295                     return -1;
 296                 }
 297             }
 298         }
 299     }
 300
 301     FFTWPREFIX(free)(p1);
 302     FFTWPREFIX(free)(p2);
 303
 304     fft->real_transform = 1;
 305     fft->ndim           = 1;
 306
 307     *pfft = fft;
 308     FFTW_UNLOCK;
 309     return 0;
 310 }
 311
 312
 313
 314 int
 315 gmx_fft_init_2d(gmx_fft_t *        pfft,
 316                 int                nx,
 317                 int                ny,
 318                 gmx_fft_flag       flags)
 319 {
 320     gmx_fft_t              fft;
 321     FFTWPREFIX(complex)   *p1,*p2,*up1,*up2;
 322     size_t                 pc;
 323     int                   i,j,k;
 324     int                    fftw_flags;
 325
 326 #ifdef GMX_DISABLE_FFTW_MEASURE
 327     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 328 #endif
 329
 330     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 331
 332     if(pfft==NULL)
 333     {
 334         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 335         return EINVAL;
 336     }
 337     *pfft = NULL;
 338
 339     FFTW_LOCK;
 340     if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
 341     {
 342         FFTW_UNLOCK;
 343         return ENOMEM;
 344     }
 345
 346     /* allocate aligned, and extra memory to make it unaligned */
 347     p1  = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny+2));
 348     if(p1==NULL)
 349     {
 350         FFTWPREFIX(free)(fft);
 351         FFTW_UNLOCK;
 352         return ENOMEM;
 353     }
 354
 355     p2  = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny+2));
 356     if(p2==NULL)
 357     {
 358         FFTWPREFIX(free)(p1);
 359         FFTWPREFIX(free)(fft);
 360         FFTW_UNLOCK;
 361         return ENOMEM;
 362     }
 363
 364     /* make unaligned pointers.
 365      * In double precision the actual complex datatype will be 16 bytes,
 366      * so go to a char pointer and force an offset of 8 bytes instead.
 367      */
 368     pc = (size_t)p1;
 369     pc += 8;
 370     up1 = (FFTWPREFIX(complex) *)pc;
 371
 372     pc = (size_t)p2;
 373     pc += 8;
 374     up2 = (FFTWPREFIX(complex) *)pc;
 375
 376
 377     fft->plan[0][0][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up2,FFTW_BACKWARD,fftw_flags);
 378     fft->plan[0][0][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up2,FFTW_FORWARD,fftw_flags);
 379     fft->plan[0][1][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up1,FFTW_BACKWARD,fftw_flags);
 380     fft->plan[0][1][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up1,FFTW_FORWARD,fftw_flags);
 381
 382     fft->plan[1][0][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p2,FFTW_BACKWARD,fftw_flags);
 383     fft->plan[1][0][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p2,FFTW_FORWARD,fftw_flags);
 384     fft->plan[1][1][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p1,FFTW_BACKWARD,fftw_flags);
 385     fft->plan[1][1][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p1,FFTW_FORWARD,fftw_flags);
 386
 387
 388     for(i=0;i<2;i++)
 389     {
 390         for(j=0;j<2;j++)
 391         {
 392             for(k=0;k<2;k++)
 393             {
 394                 if(fft->plan[i][j][k] == NULL)
 395                 {
 396                     gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
 397                     FFTW_UNLOCK;
 398                     gmx_fft_destroy(fft);
 399                     FFTW_LOCK;
 400                     FFTWPREFIX(free)(p1);
 401                     FFTWPREFIX(free)(p2);
 402                     FFTW_UNLOCK;
 403                     return -1;
 404                 }
 405             }
 406         }
 407     }
 408
 409     FFTWPREFIX(free)(p1);
 410     FFTWPREFIX(free)(p2);
 411
 412     fft->real_transform = 0;
 413     fft->ndim           = 2;
 414
 415     *pfft = fft;
 416     FFTW_UNLOCK;
 417     return 0;
 418 }
 419
 420
 421
 422 int
 423 gmx_fft_init_2d_real(gmx_fft_t *        pfft,
 424                      int                nx,
 425                      int                ny,
 426                      gmx_fft_flag       flags)
 427 {
 428     gmx_fft_t              fft;
 429     real            *p1,*p2,*up1,*up2;
 430     size_t                pc;
 431     int                   i,j,k;
 432     int                    fftw_flags;
 433
 434 #ifdef GMX_DISABLE_FFTW_MEASURE
 435     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 436 #endif
 437
 438     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 439
 440     if(pfft==NULL)
 441     {
 442         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 443         return EINVAL;
 444     }
 445     *pfft = NULL;
 446
 447     FFTW_LOCK;
 448     if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
 449     {
 450         FFTW_UNLOCK;
 451         return ENOMEM;
 452     }
 453
 454     /* allocate aligned, and extra memory to make it unaligned */
 455     p1  = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*(ny/2+1)*2 + 2) );
 456     if(p1==NULL)
 457     {
 458         FFTWPREFIX(free)(fft);
 459         FFTW_UNLOCK;
 460         return ENOMEM;
 461     }
 462
 463     p2  = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*(ny/2+1)*2 + 2) );
 464     if(p2==NULL)
 465     {
 466         FFTWPREFIX(free)(p1);
 467         FFTWPREFIX(free)(fft);
 468         FFTW_UNLOCK;
 469         return ENOMEM;
 470     }
 471
 472     /* make unaligned pointers.
 473      * In double precision the actual complex datatype will be 16 bytes,
 474      * so go to a char pointer and force an offset of 8 bytes instead.
 475      */
 476     pc = (size_t)p1;
 477     pc += 8;
 478     up1 = (real *)pc;
 479
 480     pc = (size_t)p2;
 481     pc += 8;
 482     up2 = (real *)pc;
 483
 484
 485     fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)up1,up2,fftw_flags);
 486     fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,up1,(FFTWPREFIX(complex) *)up2,fftw_flags);
 487     fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)up1,up1,fftw_flags);
 488     fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,up1,(FFTWPREFIX(complex) *)up1,fftw_flags);
 489
 490     fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)p1,p2,fftw_flags);
 491     fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,p1,(FFTWPREFIX(complex) *)p2,fftw_flags);
 492     fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)p1,p1,fftw_flags);
 493     fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,p1,(FFTWPREFIX(complex) *)p1,fftw_flags);
 494
 495
 496     for(i=0;i<2;i++)
 497     {
 498         for(j=0;j<2;j++)
 499         {
 500             for(k=0;k<2;k++)
 501             {
 502                 if(fft->plan[i][j][k] == NULL)
 503                 {
 504                     gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
 505                     FFTW_UNLOCK;
 506                     gmx_fft_destroy(fft);
 507                     FFTW_LOCK;
 508                     FFTWPREFIX(free)(p1);
 509                     FFTWPREFIX(free)(p2);
 510                     FFTW_UNLOCK;
 511                     return -1;
 512                 }
 513             }
 514         }
 515     }
 516
 517     FFTWPREFIX(free)(p1);
 518     FFTWPREFIX(free)(p2);
 519
 520     fft->real_transform = 1;
 521     fft->ndim           = 2;
 522
 523     *pfft = fft;
 524     FFTW_UNLOCK;
 525     return 0;
 526 }
 527
 528
 529
 530 int
 531 gmx_fft_init_3d(gmx_fft_t *        pfft,
 532                 int                nx,
 533                 int                ny,
 534                 int                nz,
 535                 gmx_fft_flag       flags)
 536 {
 537     gmx_fft_t              fft;
 538     FFTWPREFIX(complex)   *p1,*p2,*up1,*up2;
 539     size_t                 pc;
 540     int                   i,j,k;
 541     int                    fftw_flags;
 542
 543 #ifdef GMX_DISABLE_FFTW_MEASURE
 544     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 545 #endif
 546
 547     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 548
 549     if(pfft==NULL)
 550     {
 551         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 552         return EINVAL;
 553     }
 554     *pfft = NULL;
 555
 556     FFTW_LOCK;
 557     if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
 558     {
 559         FFTW_UNLOCK;
 560         return ENOMEM;
 561     }
 562
 563     /* allocate aligned, and extra memory to make it unaligned */
 564     p1  = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny*nz+2));
 565     if(p1==NULL)
 566     {
 567         FFTWPREFIX(free)(fft);
 568         FFTW_UNLOCK;
 569         return ENOMEM;
 570     }
 571
 572     p2  = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny*nz+2));
 573     if(p2==NULL)
 574     {
 575         FFTWPREFIX(free)(p1);
 576         FFTWPREFIX(free)(fft);
 577         FFTW_UNLOCK;
 578         return ENOMEM;
 579     }
 580
 581     /* make unaligned pointers.
 582         * In double precision the actual complex datatype will be 16 bytes,
 583         * so go to a char pointer and force an offset of 8 bytes instead.
 584         */
 585     pc = (size_t)p1;
 586     pc += 8;
 587     up1 = (FFTWPREFIX(complex) *)pc;
 588
 589     pc = (size_t)p2;
 590     pc += 8;
 591     up2 = (FFTWPREFIX(complex) *)pc;
 592
 593
 594     fft->plan[0][0][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up2,FFTW_BACKWARD,fftw_flags);
 595     fft->plan[0][0][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up2,FFTW_FORWARD,fftw_flags);
 596     fft->plan[0][1][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up1,FFTW_BACKWARD,fftw_flags);
 597     fft->plan[0][1][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up1,FFTW_FORWARD,fftw_flags);
 598
 599     fft->plan[1][0][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p2,FFTW_BACKWARD,fftw_flags);
 600     fft->plan[1][0][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p2,FFTW_FORWARD,fftw_flags);
 601     fft->plan[1][1][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p1,FFTW_BACKWARD,fftw_flags);
 602     fft->plan[1][1][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p1,FFTW_FORWARD,fftw_flags);
 603
 604
 605     for(i=0;i<2;i++)
 606     {
 607         for(j=0;j<2;j++)
 608         {
 609             for(k=0;k<2;k++)
 610             {
 611                 if(fft->plan[i][j][k] == NULL)
 612                 {
 613                     gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
 614                     FFTW_UNLOCK;
 615                     gmx_fft_destroy(fft);
 616                     FFTW_LOCK;
 617                     FFTWPREFIX(free)(p1);
 618                     FFTWPREFIX(free)(p2);
 619                     FFTW_UNLOCK;
 620                     return -1;
 621                 }
 622             }
 623         }
 624     }
 625
 626     FFTWPREFIX(free)(p1);
 627     FFTWPREFIX(free)(p2);
 628
 629     fft->real_transform = 0;
 630     fft->ndim           = 3;
 631
 632     *pfft = fft;
 633     FFTW_UNLOCK;
 634     return 0;
 635 }
 636
 637
 638
 639 int
 640 gmx_fft_init_3d_real(gmx_fft_t *        pfft,
 641                      int                nx,
 642                      int                ny,
 643                      int                nz,
 644                      gmx_fft_flag       flags)
 645 {
 646     gmx_fft_t             fft;
 647     real            *p1,*p2,*up1,*up2;
 648     size_t                pc;
 649     int                   i,j,k;
 650     int                    fftw_flags;
 651
 652 #ifdef GMX_DISABLE_FFTW_MEASURE
 653     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 654 #endif
 655
 656     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 657
 658     if(pfft==NULL)
 659     {
 660         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 661         return EINVAL;
 662     }
 663     *pfft = NULL;
 664
 665     FFTW_LOCK;
 666     if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
 667     {
 668         FFTW_UNLOCK;
 669         return ENOMEM;
 670     }
 671
 672     /* allocate aligned, and extra memory to make it unaligned */
 673     p1  = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*ny*(nz/2+1)*2 + 2) );
 674     if(p1==NULL)
 675     {
 676         FFTWPREFIX(free)(fft);
 677         FFTW_UNLOCK;
 678         return ENOMEM;
 679     }
 680
 681     p2  = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*ny*(nz/2+1)*2 + 2) );
 682     if(p2==NULL)
 683     {
 684         FFTWPREFIX(free)(p1);
 685         FFTWPREFIX(free)(fft);
 686         FFTW_UNLOCK;
 687         return ENOMEM;
 688     }
 689
 690     /* make unaligned pointers.
 691      * In double precision the actual complex datatype will be 16 bytes,
 692      * so go to a void pointer and force an offset of 8 bytes instead.
 693      */
 694     pc = (size_t)p1;
 695     pc += 8;
 696     up1 = (real *)pc;
 697
 698     pc = (size_t)p2;
 699     pc += 8;
 700     up2 = (real *)pc;
 701
 702
 703     fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)up1,up2,fftw_flags);
 704     fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,up1,(FFTWPREFIX(complex) *)up2,fftw_flags);
 705     fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)up1,up1,fftw_flags);
 706     fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,up1,(FFTWPREFIX(complex) *)up1,fftw_flags);
 707
 708     fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)p1,p2,fftw_flags);
 709     fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,p1,(FFTWPREFIX(complex) *)p2,fftw_flags);
 710     fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)p1,p1,fftw_flags);
 711     fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,p1,(FFTWPREFIX(complex) *)p1,fftw_flags);
 712
 713
 714     for(i=0;i<2;i++)
 715     {
 716         for(j=0;j<2;j++)
 717         {
 718             for(k=0;k<2;k++)
 719             {
 720                 if(fft->plan[i][j][k] == NULL)
 721                 {
 722                     gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
 723                     FFTW_UNLOCK;
 724                     gmx_fft_destroy(fft);
 725                     FFTW_LOCK;
 726                     FFTWPREFIX(free)(p1);
 727                     FFTWPREFIX(free)(p2);
 728                     FFTW_UNLOCK;
 729                     return -1;
 730                 }
 731             }
 732         }
 733     }
 734
 735     FFTWPREFIX(free)(p1);
 736     FFTWPREFIX(free)(p2);
 737
 738     fft->real_transform = 1;
 739     fft->ndim           = 3;
 740
 741     *pfft = fft;
 742     FFTW_UNLOCK;
 743     return 0;
 744 }
 745
 746
 747 int
 748 gmx_fft_1d               (gmx_fft_t                  fft,
 749                           enum gmx_fft_direction     dir,
 750                           void *                     in_data,
 751                           void *                     out_data)
 752 {
 753     int           aligned   = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
 754     int           inplace   = (in_data == out_data);
 755     int           isforward = (dir == GMX_FFT_FORWARD);
 756
 757     /* Some checks */
 758     if( (fft->real_transform == 1) || (fft->ndim != 1) ||
 759         ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
 760     {
 761         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 762         return EINVAL;
 763     }
 764
 765     FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
 766                             (FFTWPREFIX(complex) *)in_data,
 767                             (FFTWPREFIX(complex) *)out_data);
 768
 769     return 0;
 770 }
 771
 772 int
 773 gmx_fft_many_1d               (gmx_fft_t                  fft,
 774                                enum gmx_fft_direction     dir,
 775                                void *                     in_data,
 776                                void *                     out_data)
 777 {
 778     return gmx_fft_1d(fft,dir,in_data,out_data);
 779 }
 780
 781 int
 782 gmx_fft_1d_real          (gmx_fft_t                  fft,
 783                           enum gmx_fft_direction     dir,
 784                           void *                     in_data,
 785                           void *                     out_data)
 786 {
 787     int           aligned   = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
 788     int           inplace   = (in_data == out_data);
 789     int           isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
 790
 791     /* Some checks */
 792     if( (fft->real_transform != 1) || (fft->ndim != 1) ||
 793         ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 794     {
 795         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 796         return EINVAL;
 797     }
 798
 799     if(isforward)
 800     {
 801         FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
 802                                     (real *)in_data,(FFTWPREFIX(complex) *)out_data);
 803     }
 804     else
 805     {
 806         FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
 807                                     (FFTWPREFIX(complex) *)in_data,(real *)out_data);
 808     }
 809
 810     return 0;
 811 }
 812
 813 int
 814 gmx_fft_many_1d_real     (gmx_fft_t                  fft,
 815                           enum gmx_fft_direction     dir,
 816                           void *                     in_data,
 817                           void *                     out_data)
 818 {
 819     return gmx_fft_1d_real(fft,dir,in_data,out_data);
 820 }
 821
 822 int
 823 gmx_fft_2d               (gmx_fft_t                  fft,
 824                           enum gmx_fft_direction     dir,
 825                           void *                     in_data,
 826                           void *                     out_data)
 827 {
 828     int           aligned   = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
 829     int           inplace   = (in_data == out_data);
 830     int           isforward = (dir == GMX_FFT_FORWARD);
 831
 832     /* Some checks */
 833     if( (fft->real_transform == 1) || (fft->ndim != 2) ||
 834         ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
 835     {
 836         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 837         return EINVAL;
 838     }
 839
 840     FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
 841                             (FFTWPREFIX(complex) *)in_data,
 842                             (FFTWPREFIX(complex) *)out_data);
 843
 844     return 0;
 845 }
 846
 847
 848 int
 849 gmx_fft_2d_real          (gmx_fft_t                  fft,
 850                           enum gmx_fft_direction     dir,
 851                           void *                     in_data,
 852                           void *                     out_data)
 853 {
 854     int           aligned   = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
 855     int           inplace   = (in_data == out_data);
 856     int           isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
 857
 858     /* Some checks */
 859     if( (fft->real_transform != 1) || (fft->ndim != 2) ||
 860         ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 861     {
 862         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 863         return EINVAL;
 864     }
 865
 866     if(isforward)
 867     {
 868         FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
 869                                     (real *)in_data,
 870                                     (FFTWPREFIX(complex) *)out_data);
 871     }
 872     else
 873     {
 874         FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
 875                                     (FFTWPREFIX(complex) *)in_data,
 876                                     (real *)out_data);
 877     }
 878
 879
 880     return 0;
 881 }
 882
 883
 884 int
 885 gmx_fft_3d               (gmx_fft_t                  fft,
 886                           enum gmx_fft_direction     dir,
 887                           void *                     in_data,
 888                           void *                     out_data)
 889 {
 890     int           aligned   = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
 891     int           inplace   = (in_data == out_data);
 892     int           isforward = (dir == GMX_FFT_FORWARD);
 893
 894     /* Some checks */
 895     if( (fft->real_transform == 1) || (fft->ndim != 3) ||
 896         ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
 897     {
 898         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 899         return EINVAL;
 900     }
 901
 902     FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
 903                             (FFTWPREFIX(complex) *)in_data,
 904                             (FFTWPREFIX(complex) *)out_data);
 905
 906     return 0;
 907 }
 908
 909
 910 int
 911 gmx_fft_3d_real          (gmx_fft_t                  fft,
 912                           enum gmx_fft_direction     dir,
 913                           void *                     in_data,
 914                           void *                     out_data)
 915 {
 916     int           aligned   = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
 917     int           inplace   = (in_data == out_data);
 918     int           isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
 919
 920     /* Some checks */
 921     if( (fft->real_transform != 1) || (fft->ndim != 3) ||
 922         ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 923     {
 924         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 925         return EINVAL;
 926     }
 927
 928     if(isforward)
 929     {
 930         FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
 931                                     (real *)in_data,
 932                                     (FFTWPREFIX(complex) *)out_data);
 933     }
 934     else
 935     {
 936         FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
 937                                     (FFTWPREFIX(complex) *)in_data,
 938                                     (real *)out_data);
 939     }
 940
 941
 942     return 0;
 943 }
 944
 945
 946 void
 947 gmx_fft_destroy(gmx_fft_t      fft)
 948 {
 949     int                   i,j,k;
 950
 951     if(fft != NULL)
 952     {
 953         for(i=0;i<2;i++)
 954         {
 955             for(j=0;j<2;j++)
 956             {
 957                 for(k=0;k<2;k++)
 958                 {
 959                     if(fft->plan[i][j][k] != NULL)
 960                     {
 961                         FFTW_LOCK;
 962                         FFTWPREFIX(destroy_plan)(fft->plan[i][j][k]);
 963                         FFTW_UNLOCK;
 964                         fft->plan[i][j][k] = NULL;
 965                     }
 966                 }
 967             }
 968         }
 969         FFTW_LOCK;
 970         FFTWPREFIX(free)(fft);
 971         FFTW_UNLOCK;
 972     }
 973
 974 }
 975
 976 void
 977 gmx_many_fft_destroy(gmx_fft_t    fft)
 978 {
 979     gmx_fft_destroy(fft);
 980 }
 981
 982 #else
 983 int
 984 gmx_fft_fftw3_empty;
 985 #endif /* GMX_FFT_FFTW3 */