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