src/mdlib/gmx_fft_fftw3.c

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