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