src/mdlib/gmx_fft_fftw2.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  * Gromacs 4.0                         Copyright (c) 1991-2003
   5  * David van der Spoel, Erik Lindahl, University of Groningen.
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation; either version 2
  10  * of the License, or (at your option) any later version.
  11  *
  12  * To help us fund GROMACS development, we humbly ask that you cite
  13  * the research papers on the package. Check out http://www.gromacs.org
  14  *
  15  * And Hey:
  16  * Gnomes, ROck Monsters And Chili Sauce
  17  */
  18 #ifdef HAVE_CONFIG_H
  19 #include <config.h>
  20 #endif
  21
  22 #ifdef GMX_FFT_FFTW2
  23
  24 #include <errno.h>
  25 #include <string.h>
  26 #include <stdlib.h>
  27
  28
  29 #include "gmx_fft.h"
  30 #include "gmx_fatal.h"
  31
  32
  33
  34 #ifdef FFTW2_NAME_FFTW
  35 #  include<fftw.h>
  36 #  include<rfftw.h>
  37 #elif defined FFTW2_NAME_SFFTW
  38 #  include<sfftw.h>
  39 #  include<srfftw.h>
  40 #elif defined FFTW2_NAME_DFFTW
  41 #  include<dfftw.h>
  42 #  include<drfftw.h>
  43 #else
  44 #error No FFTW2 name defined - you must define one of:
  45 #error FFTW2_NAME_FFTW, FFTW2_NAME_SFFTW, FFTW2_NAME_DFFTW
  46 #endif
  47
  48
  49 /* Contents of the FFTW2 setup */
  50 struct gmx_fft
  51 {
  52     int               ndim;         /**< Number of dimensions in transform.   */
  53     int               nx;           /**< Data X dimension                     */
  54     int               ny;           /**< Data Y dimension                     */
  55     int               nz;           /**< Data Z dimension                     */
  56     /* Arrays with fftw2 plans.
  57      * First index is 0 for out-of-place, 1 for in-place transform.
  58      * Second index is 0 for backward, 1 for forward.
  59      */
  60     fftw_plan         single[2][2]; /**< Plans for 1d transforms.             */
  61     fftwnd_plan       multi[2][2];  /**< Plans for n-d transforms.            */
  62     real *            work;         /**< Avoid overwriting input for c2r ffts */
  63 };
  64
  65
  66
  67 int
  68 gmx_fft_init_1d(gmx_fft_t *        pfft,
  69                 int                nx,
  70                 enum gmx_fft_flag  flags)
  71 {
  72     int i,j;
  73     gmx_fft_t           fft;
  74     int                    fftw_flags;
  75
  76     /* FFTW2 is slow to measure, so we do not use it */
  77
  78     /* If you change this, add an #ifndef for GMX_DISABLE_FFTW_MEASURE around it! */
  79     fftw_flags = FFTW_ESTIMATE;
  80
  81     if(pfft==NULL)
  82     {
  83         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
  84         return EINVAL;
  85     }
  86     *pfft = NULL;
  87
  88     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
  89     {
  90         return ENOMEM;
  91     }
  92
  93
  94     fft->single[0][0] = fftw_create_plan(nx,FFTW_BACKWARD,FFTW_OUT_OF_PLACE|fftw_flags);
  95     fft->single[0][1] = fftw_create_plan(nx,FFTW_FORWARD,FFTW_OUT_OF_PLACE|fftw_flags);
  96     fft->single[1][0] = fftw_create_plan(nx,FFTW_BACKWARD,FFTW_IN_PLACE|fftw_flags);
  97     fft->single[1][1] = fftw_create_plan(nx,FFTW_FORWARD,FFTW_IN_PLACE|fftw_flags);
  98
  99
 100     fft->multi[0][0] = NULL;
 101     fft->multi[0][1] = NULL;
 102     fft->multi[1][0] = NULL;
 103     fft->multi[1][1] = NULL;
 104
 105     for(i=0;i<2;i++)
 106     {
 107         for(j=0;j<2;j++)
 108         {
 109             if(fft->single[i][j] == NULL)
 110             {
 111                 gmx_fatal(FARGS,"Error initializing FFTW2 plan.");
 112                 gmx_fft_destroy(fft);
 113                 return -1;
 114             }
 115         }
 116     }
 117
 118     /* No workspace needed for complex-to-complex FFTs */
 119     fft->work = NULL;
 120
 121     fft->ndim = 1;
 122     fft->nx   = nx;
 123
 124     *pfft = fft;
 125     return 0;
 126 }
 127
 128
 129
 130 int
 131 gmx_fft_init_1d_real(gmx_fft_t *        pfft,
 132                      int                nx,
 133                      enum gmx_fft_flag  flags)
 134 {
 135     int i,j;
 136     gmx_fft_t             fft;
 137     int                    fftw_flags;
 138
 139     /* FFTW2 is slow to measure, so we do not use it */
 140     /* If you change this, add an #ifndef for GMX_DISABLE_FFTW_MEASURE around it! */
 141     fftw_flags = FFTW_ESTIMATE;
 142
 143     if(pfft==NULL)
 144     {
 145         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 146         return EINVAL;
 147     }
 148     *pfft = NULL;
 149
 150     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 151     {
 152         return ENOMEM;
 153     }
 154
 155
 156     fft->single[0][0] = rfftw_create_plan(nx,FFTW_COMPLEX_TO_REAL,FFTW_OUT_OF_PLACE|fftw_flags);
 157     fft->single[0][1] = rfftw_create_plan(nx,FFTW_REAL_TO_COMPLEX,FFTW_OUT_OF_PLACE|fftw_flags);
 158     fft->single[1][0] = rfftw_create_plan(nx,FFTW_COMPLEX_TO_REAL,FFTW_IN_PLACE|fftw_flags);
 159     fft->single[1][1] = rfftw_create_plan(nx,FFTW_REAL_TO_COMPLEX,FFTW_IN_PLACE|fftw_flags);
 160
 161
 162     fft->multi[0][0] = NULL;
 163     fft->multi[0][1] = NULL;
 164     fft->multi[1][0] = NULL;
 165     fft->multi[1][1] = NULL;
 166
 167     for(i=0;i<2;i++)
 168     {
 169         for(j=0;j<2;j++)
 170         {
 171             if(fft->single[i][j] == NULL)
 172             {
 173                 gmx_fatal(FARGS,"Error initializing FFTW2 plan.");
 174                 gmx_fft_destroy(fft);
 175                 return -1;
 176             }
 177         }
 178     }
 179
 180     /* FFTW2 overwrites the input when doing out-of-place complex-to-real FFTs.
 181      * This is not acceptable for the Gromacs interface, so we define a
 182      * work array and copy the data there before doing complex-to-real FFTs.
 183      */
 184     fft->work = (real *)malloc(sizeof(real)*( (nx/2 + 1)*2) );
 185     if(fft->work == NULL)
 186     {
 187         gmx_fatal(FARGS,"Cannot allocate complex-to-real FFT workspace.");
 188         gmx_fft_destroy(fft);
 189         return ENOMEM;
 190     }
 191
 192     fft->ndim = 1;
 193     fft->nx   = nx;
 194
 195     *pfft = fft;
 196     return 0;
 197 }
 198
 199
 200
 201 int
 202 gmx_fft_init_2d(gmx_fft_t *        pfft,
 203                 int                nx,
 204                 int                ny,
 205                 enum gmx_fft_flag  flags)
 206 {
 207     int i,j;
 208     gmx_fft_t             fft;
 209     int                    fftw_flags;
 210
 211
 212     /* FFTW2 is slow to measure, so we do not use it */
 213     /* If you change this, add an #ifndef for GMX_DISABLE_FFTW_MEASURE around it! */
 214     fftw_flags = FFTW_ESTIMATE;
 215
 216     if(pfft==NULL)
 217     {
 218         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 219         return EINVAL;
 220     }
 221     *pfft = NULL;
 222
 223     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 224     {
 225         return ENOMEM;
 226     }
 227
 228     fft->single[0][0] = NULL;
 229     fft->single[0][1] = NULL;
 230     fft->single[1][0] = NULL;
 231     fft->single[1][1] = NULL;
 232
 233     fft->multi[0][0] = fftw2d_create_plan(nx,ny,FFTW_BACKWARD,FFTW_OUT_OF_PLACE|fftw_flags);
 234     fft->multi[0][1] = fftw2d_create_plan(nx,ny,FFTW_FORWARD,FFTW_OUT_OF_PLACE|fftw_flags);
 235     fft->multi[1][0] = fftw2d_create_plan(nx,ny,FFTW_BACKWARD,FFTW_IN_PLACE|fftw_flags);
 236     fft->multi[1][1] = fftw2d_create_plan(nx,ny,FFTW_FORWARD,FFTW_IN_PLACE|fftw_flags);
 237
 238     for(i=0;i<2;i++)
 239     {
 240         for(j=0;j<2;j++)
 241         {
 242             if(fft->multi[i][j] == NULL)
 243             {
 244                 gmx_fatal(FARGS,"Error initializing FFTW2 plan.");
 245                 gmx_fft_destroy(fft);
 246                 return -1;
 247             }
 248         }
 249     }
 250
 251     /* No workspace needed for complex-to-complex FFTs */
 252     fft->work = NULL;
 253
 254     fft->ndim = 2;
 255     fft->nx   = nx;
 256     fft->ny   = ny;
 257
 258     *pfft = fft;
 259     return 0;
 260 }
 261
 262
 263
 264
 265 int
 266 gmx_fft_init_2d_real(gmx_fft_t *        pfft,
 267                      int                nx,
 268                      int                ny,
 269                      enum gmx_fft_flag  flags)
 270 {
 271     int i,j;
 272     gmx_fft_t             fft;
 273     int                    fftw_flags;
 274
 275
 276     /* FFTW2 is slow to measure, so we do not use it */
 277     /* If you change this, add an #ifndef for GMX_DISABLE_FFTW_MEASURE around it! */
 278     fftw_flags = FFTW_ESTIMATE;
 279
 280     if(pfft==NULL)
 281     {
 282         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 283         return EINVAL;
 284     }
 285     *pfft = NULL;
 286
 287     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 288     {
 289         return ENOMEM;
 290     }
 291
 292     fft->single[0][0] = NULL;
 293     fft->single[0][1] = NULL;
 294     fft->single[1][0] = NULL;
 295     fft->single[1][1] = NULL;
 296
 297
 298     fft->multi[0][0] = rfftw2d_create_plan(nx,ny,FFTW_COMPLEX_TO_REAL,FFTW_OUT_OF_PLACE|fftw_flags);
 299     fft->multi[0][1] = rfftw2d_create_plan(nx,ny,FFTW_REAL_TO_COMPLEX,FFTW_OUT_OF_PLACE|fftw_flags);
 300     fft->multi[1][0] = rfftw2d_create_plan(nx,ny,FFTW_COMPLEX_TO_REAL,FFTW_IN_PLACE|fftw_flags);
 301     fft->multi[1][1] = rfftw2d_create_plan(nx,ny,FFTW_REAL_TO_COMPLEX,FFTW_IN_PLACE|fftw_flags);
 302
 303
 304     for(i=0;i<2;i++)
 305     {
 306         for(j=0;j<2;j++)
 307         {
 308             if(fft->multi[i][j] == NULL)
 309             {
 310                 gmx_fatal(FARGS,"Error initializing FFTW2 plan.");
 311                 gmx_fft_destroy(fft);
 312                 return -1;
 313             }
 314         }
 315     }
 316
 317     /* FFTW2 overwrites the input when doing out-of-place complex-to-real FFTs.
 318      * This is not acceptable for the Gromacs interface, so we define a
 319      * work array and copy the data there before doing complex-to-real FFTs.
 320      */
 321     fft->work = (real *)malloc(sizeof(real)*( nx*(ny/2 + 1)*2) );
 322     if(fft->work == NULL)
 323     {
 324         gmx_fatal(FARGS,"Cannot allocate complex-to-real FFT workspace.");
 325         gmx_fft_destroy(fft);
 326         return ENOMEM;
 327     }
 328
 329
 330     fft->ndim = 2;
 331     fft->nx   = nx;
 332     fft->ny   = ny;
 333
 334     *pfft = fft;
 335     return 0;
 336 }
 337
 338
 339 int
 340 gmx_fft_init_3d(gmx_fft_t *        pfft,
 341                 int                nx,
 342                 int                ny,
 343                 int                nz,
 344                 enum gmx_fft_flag  flags)
 345 {
 346     int i,j;
 347     gmx_fft_t             fft;
 348     int                    fftw_flags;
 349
 350
 351     /* FFTW2 is slow to measure, so we do not use it */
 352     /* If you change this, add an #ifndef for GMX_DISABLE_FFTW_MEASURE around it! */
 353     fftw_flags = FFTW_ESTIMATE;
 354
 355     if(pfft==NULL)
 356     {
 357         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 358         return EINVAL;
 359     }
 360     *pfft = NULL;
 361
 362     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 363     {
 364         return ENOMEM;
 365     }
 366
 367     fft->single[0][0] = NULL;
 368     fft->single[0][1] = NULL;
 369     fft->single[1][0] = NULL;
 370     fft->single[1][1] = NULL;
 371
 372
 373     fft->multi[0][0] = fftw3d_create_plan(nx,ny,nz,FFTW_BACKWARD,FFTW_OUT_OF_PLACE|fftw_flags);
 374     fft->multi[0][1] = fftw3d_create_plan(nx,ny,nz,FFTW_FORWARD,FFTW_OUT_OF_PLACE|fftw_flags);
 375     fft->multi[1][0] = fftw3d_create_plan(nx,ny,nz,FFTW_BACKWARD,FFTW_IN_PLACE|fftw_flags);
 376     fft->multi[1][1] = fftw3d_create_plan(nx,ny,nz,FFTW_FORWARD,FFTW_IN_PLACE|fftw_flags);
 377
 378
 379     for(i=0;i<2;i++)
 380     {
 381         for(j=0;j<2;j++)
 382         {
 383             if(fft->multi[i][j] == NULL)
 384             {
 385                 gmx_fatal(FARGS,"Error initializing FFTW2 plan.");
 386                 gmx_fft_destroy(fft);
 387                 return -1;
 388             }
 389         }
 390     }
 391
 392     /* No workspace needed for complex-to-complex FFTs */
 393     fft->work = NULL;
 394     fft->nx   = nx;
 395     fft->ny   = ny;
 396     fft->nz   = nz;
 397
 398     fft->ndim = 3;
 399
 400     *pfft = fft;
 401     return 0;
 402 }
 403
 404
 405
 406
 407 int
 408 gmx_fft_init_3d_real(gmx_fft_t *        pfft,
 409                      int                nx,
 410                      int                ny,
 411                      int                nz,
 412                      enum gmx_fft_flag  flags)
 413 {
 414     int i,j;
 415     gmx_fft_t            fft;
 416     int                    fftw_flags;
 417
 418
 419     /* FFTW2 is slow to measure, so we do not use it */
 420     /* If you change this, add an #ifndef for GMX_DISABLE_FFTW_MEASURE around it! */
 421     fftw_flags = FFTW_ESTIMATE;
 422
 423     if(pfft==NULL)
 424     {
 425         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 426         return EINVAL;
 427     }
 428     *pfft = NULL;
 429
 430     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 431     {
 432         return ENOMEM;
 433     }
 434
 435     fft->single[0][0] = NULL;
 436     fft->single[0][1] = NULL;
 437     fft->single[1][0] = NULL;
 438     fft->single[1][1] = NULL;
 439
 440
 441     fft->multi[0][0] = rfftw3d_create_plan(nx,ny,nz,FFTW_COMPLEX_TO_REAL,FFTW_OUT_OF_PLACE|fftw_flags);
 442     fft->multi[0][1] = rfftw3d_create_plan(nx,ny,nz,FFTW_REAL_TO_COMPLEX,FFTW_OUT_OF_PLACE|fftw_flags);
 443     fft->multi[1][0] = rfftw3d_create_plan(nx,ny,nz,FFTW_COMPLEX_TO_REAL,FFTW_IN_PLACE|fftw_flags);
 444     fft->multi[1][1] = rfftw3d_create_plan(nx,ny,nz,FFTW_REAL_TO_COMPLEX,FFTW_IN_PLACE|fftw_flags);
 445
 446
 447     for(i=0;i<2;i++)
 448     {
 449         for(j=0;j<2;j++)
 450         {
 451             if(fft->multi[i][j] == NULL)
 452             {
 453                 gmx_fatal(FARGS,"Error initializing FFTW2 plan.");
 454                 gmx_fft_destroy(fft);
 455                 return -1;
 456             }
 457         }
 458     }
 459
 460     /* FFTW2 overwrites the input when doing out-of-place complex-to-real FFTs.
 461      * This is not acceptable for the Gromacs interface, so we define a
 462      * work array and copy the data there before doing complex-to-real FFTs.
 463      */
 464     fft->work = (real *)malloc(sizeof(real)*( nx*ny*(nz/2 + 1)*2) );
 465     if(fft->work == NULL)
 466     {
 467         gmx_fatal(FARGS,"Cannot allocate complex-to-real FFT workspace.");
 468         gmx_fft_destroy(fft);
 469         return ENOMEM;
 470     }
 471
 472     fft->ndim = 3;
 473     fft->nx   = nx;
 474     fft->ny   = ny;
 475     fft->nz   = nz;
 476
 477     *pfft = fft;
 478     return 0;
 479 }
 480
 481
 482 int
 483 gmx_fft_1d(gmx_fft_t                  fft,
 484            enum gmx_fft_direction     dir,
 485            void *                     in_data,
 486            void *                     out_data)
 487 {
 488     int inplace   = (in_data == out_data);
 489     int isforward = (dir == GMX_FFT_FORWARD);
 490
 491     if((fft->ndim != 1) ||
 492        ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)))
 493     {
 494         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 495         return EINVAL;
 496     }
 497
 498     fftw_one(fft->single[inplace][isforward],(fftw_complex *)in_data,(fftw_complex *)out_data);
 499
 500   return 0;
 501 }
 502
 503 int
 504 gmx_fft_1d_real(gmx_fft_t                  fft,
 505                 enum gmx_fft_direction     dir,
 506                 void *                     in_data,
 507                 void *                     out_data)
 508 {
 509     /* FFTW2 1-dimensional real transforms are special.
 510      *
 511      * First, the complex data is stored in a special packed half-complex
 512      * fashion. To enable a standard common Gromacs interface this forces us
 513      * to always use out-of-place FFTs, and permute the data after
 514      * real-to-complex FFTs or before complex-to-real FFTs.
 515      *
 516      * The input is also destroyed for out-of-place complex-to-real FFTs, but
 517      * this doesn't matter since we need to permute and copy the data into
 518      * the work array first anyway.
 519      */
 520     real *     work = fft->work;
 521     t_complex *  data;
 522     int              n    = fft->nx;
 523     int              i;
 524
 525     if((fft->ndim != 1) ||
 526        ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)))
 527     {
 528         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 529         return EINVAL;
 530     }
 531
 532     if(dir==GMX_FFT_REAL_TO_COMPLEX)
 533     {
 534         rfftw_one(fft->single[0][1],(fftw_real *)in_data,(fftw_real *)work);
 535         /* permute it back into data, in standard complex format
 536          * instead of halfcomplex...
 537          */
 538         data = (t_complex *)out_data;
 539
 540         data[0].re = work[0];
 541         data[0].im = 0;
 542
 543         for(i=1;i<n/2;i++)
 544         {
 545             data[i].re = work[i];
 546             data[i].im = work[n-i];
 547         }
 548
 549         data[i].re=work[i];
 550
 551         if(2*i==n)
 552         {
 553             data[i].im=0;
 554         }
 555         else
 556         {
 557             data[i].im=work[n-i];
 558         }
 559     }
 560     else
 561     {
 562         /* Complex-to-real. First permute standard format into halfcomplex */
 563         data = (t_complex *)in_data;
 564
 565         work[0]=data[0].re;
 566
 567         for(i=1;i<n/2;i++)
 568         {
 569             work[i]  =data[i].re;
 570             work[n-i]=data[i].im;
 571         }
 572
 573         if(2*i!=n)
 574         {
 575             work[n-i]=data[i].im;
 576         }
 577
 578         rfftw_one(fft->single[0][0],(fftw_real *)work,(fftw_real *)out_data);
 579     }
 580
 581     return 0;
 582 }
 583
 584
 585 int
 586 gmx_fft_2d(gmx_fft_t                  fft,
 587            enum gmx_fft_direction     dir,
 588            void *                     in_data,
 589            void *                     out_data)
 590 {
 591     int inplace   = (in_data == out_data);
 592     int isforward = (dir == GMX_FFT_FORWARD);
 593
 594     if((fft->ndim != 2) ||
 595        ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)))
 596     {
 597         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 598         return EINVAL;
 599     }
 600
 601     fftwnd_one(fft->multi[inplace][isforward],(fftw_complex *)in_data,(fftw_complex *)out_data);
 602
 603     return 0;
 604 }
 605
 606
 607 int
 608 gmx_fft_2d_real(gmx_fft_t                  fft,
 609                 enum gmx_fft_direction     dir,
 610                 void *                     in_data,
 611                 void *                     out_data)
 612 {
 613     int inplace   = (in_data == out_data);
 614     int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
 615     int sz;
 616
 617     if((fft->ndim != 2) ||
 618        ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)))
 619     {
 620         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 621         return EINVAL;
 622     }
 623
 624     if(inplace == 0)
 625     {
 626         /* Copy data to avoid overwriting input, and redirect input ptr to work array */
 627         sz = fft->nx*(fft->ny/2 + 1)*2;
 628         memcpy(fft->work,in_data,sz*sizeof(real));
 629         in_data = fft->work;
 630     }
 631
 632     if(isforward)
 633     {
 634         rfftwnd_one_real_to_complex(fft->multi[inplace][isforward],(fftw_real *)in_data,(fftw_complex *)out_data);
 635     }
 636     else
 637     {
 638         rfftwnd_one_complex_to_real(fft->multi[inplace][isforward],(fftw_complex *)in_data,(fftw_real *)out_data);
 639     }
 640
 641     return 0;
 642 }
 643
 644
 645 int
 646 gmx_fft_3d(gmx_fft_t                  fft,
 647            enum gmx_fft_direction     dir,
 648            void *                     in_data,
 649            void *                     out_data)
 650 {
 651     int inplace   = (in_data == out_data);
 652     int isforward = (dir == GMX_FFT_FORWARD);
 653
 654     if((fft->ndim != 3) ||
 655        ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)))
 656     {
 657         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 658         return EINVAL;
 659     }
 660
 661     fftwnd_one(fft->multi[inplace][isforward],(fftw_complex *)in_data,(fftw_complex *)out_data);
 662
 663     return 0;
 664 }
 665
 666
 667 int
 668 gmx_fft_3d_real(gmx_fft_t                  fft,
 669                 enum gmx_fft_direction     dir,
 670                 void *                     in_data,
 671                 void *                     out_data)
 672 {
 673     int inplace   = (in_data == out_data);
 674     int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
 675     int sz;
 676
 677     if((fft->ndim != 3) ||
 678        ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)))
 679     {
 680         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 681         return EINVAL;
 682     }
 683
 684     if(inplace == 0)
 685     {
 686         /* Copy data to avoid overwriting input, and redirect input ptr to work array */
 687         sz = fft->nx*fft->ny*(fft->nz/2 + 1)*2;
 688         memcpy(fft->work,in_data,sz*sizeof(real));
 689         in_data = fft->work;
 690     }
 691
 692     if(isforward)
 693     {
 694         rfftwnd_one_real_to_complex(fft->multi[inplace][isforward],(fftw_real *)in_data,(fftw_complex *)out_data);
 695     }
 696     else
 697     {
 698         rfftwnd_one_complex_to_real(fft->multi[inplace][isforward],(fftw_complex *)in_data,(fftw_real *)out_data);
 699     }
 700
 701     return 0;
 702 }
 703
 704
 705
 706
 707 void
 708 gmx_fft_destroy(gmx_fft_t    fft)
 709 {
 710     int i,j;
 711
 712     if(fft != NULL)
 713     {
 714         for(i=0;i<2;i++)
 715         {
 716             for(j=0;j<2;j++)
 717             {
 718                 if(fft->single[i][j] != NULL)
 719                 {
 720                     rfftw_destroy_plan(fft->single[i][j]);
 721                     fft->single[i][j] = NULL;
 722                 }
 723                 if(fft->multi[i][j] != NULL)
 724                 {
 725                     rfftwnd_destroy_plan(fft->multi[i][j]);
 726                     fft->multi[i][j] = NULL;
 727                 }
 728             }
 729         }
 730         free(fft);
 731     }
 732 }
 733
 734
 735 #else
 736 int
 737 gmx_fft_fftw2_empty;
 738 #endif /* GMX_FFT_FFTW2 */