src/mdlib/gmx_fft_mkl.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_MKL
  42
  43 #include <errno.h>
  44 #include <stdlib.h>
  45
  46 #include <mkl_dfti.h>
  47
  48
  49 #include "gmx_fft.h"
  50 #include "gmx_fatal.h"
  51
  52
  53 /* For MKL version (<10.0), we should define MKL_LONG. */
  54 #ifndef MKL_LONG
  55 #define MKL_LONG long int
  56 #endif
  57
  58
  59 #ifdef GMX_DOUBLE
  60 #define GMX_DFTI_PREC  DFTI_DOUBLE
  61 #else
  62 #define GMX_DFTI_PREC  DFTI_SINGLE
  63 #endif
  64
  65 /* Contents of the Intel MKL FFT fft datatype.
  66  *
  67  * Note that this is one of several possible implementations of gmx_fft_t.
  68  *
  69  *  The MKL _API_ supports 1D,2D, and 3D transforms, including real-to-complex.
  70  *  Unfortunately the actual library implementation does not support 3D real
  71  *  transforms as of version 7.2, and versions before 7.0 don't support 2D real
  72  *  either. In addition, the multi-dimensional storage format for real data
  73  *  is not compatible with our padding.
  74  *
  75  *  To work around this we roll our own 2D and 3D real-to-complex transforms,
  76  *  using separate X/Y/Z handles defined to perform (ny*nz), (nx*nz), and
  77  *  (nx*ny) transforms at once when necessary. To perform strided multiple
  78  *  transforms out-of-place (i.e., without padding in the last dimension)
  79  *  on the fly we also need to separate the forward and backward
  80  *  handles for real-to-complex/complex-to-real data permutation.
  81  *
  82  *  This makes it necessary to define 3 handles for in-place FFTs, and 4 for
  83  *  the out-of-place transforms. Still, whenever possible we try to use
  84  *  a single 3D-transform handle instead.
  85  *
  86  *  So, the handles are enumerated as follows:
  87  *
  88  *  1D FFT (real too):    Index 0 is the handle for the entire FFT
  89  *  2D complex FFT:       Index 0 is the handle for the entire FFT
  90  *  3D complex FFT:       Index 0 is the handle for the entire FFT
  91  *  2D, inplace real FFT: 0=FFTx, 1=FFTy handle
  92  *  2D, ooplace real FFT: 0=FFTx, 1=real-to-complex FFTy, 2=complex-to-real FFTy
  93  *  3D, inplace real FFT: 0=FFTx, 1=FFTy, 2=FFTz handle
  94  *  3D, ooplace real FFT: 0=FFTx, 1=FFTy, 2=r2c FFTz, 3=c2r FFTz
  95  *
  96  *  Intel people reading this: Learn from FFTW what a good interface looks like :-)
  97  *
  98  */
  99 struct gmx_fft
 100 {
 101     int                ndim;              /**< Number of dimensions in FFT  */
 102     int                nx;                /**< Length of X transform        */
 103     int                ny;                /**< Length of Y transform        */
 104     int                nz;                /**< Length of Z transform        */
 105     int                real_fft;          /**< 1 if real FFT, otherwise 0   */
 106     DFTI_DESCRIPTOR *  inplace[3];        /**< in-place FFT                 */
 107     DFTI_DESCRIPTOR *  ooplace[4];        /**< out-of-place FFT             */
 108     t_complex *    work;              /**< Enable out-of-place c2r FFT  */
 109 };
 110
 111
 112
 113 int
 114 gmx_fft_init_1d(gmx_fft_t *        pfft,
 115                 int                nx,
 116                 gmx_fft_flag  flags)
 117 {
 118     gmx_fft_t      fft;
 119     int            d;
 120     int            status;
 121
 122     if(pfft==NULL)
 123     {
 124         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 125         return EINVAL;
 126     }
 127     *pfft = NULL;
 128
 129     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 130     {
 131         return ENOMEM;
 132     }
 133
 134     /* Mark all handles invalid */
 135     for(d=0;d<3;d++)
 136     {
 137         fft->inplace[d] = fft->ooplace[d] = NULL;
 138     }
 139     fft->ooplace[3] = NULL;
 140
 141
 142     status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
 143
 144     if( status == 0 )
 145         status = DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE);
 146
 147     if( status == 0 )
 148         status = DftiCommitDescriptor(fft->inplace[0]);
 149
 150
 151     if( status == 0 )
 152         status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
 153
 154     if( status == 0)
 155         DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE);
 156
 157     if( status == 0)
 158         DftiCommitDescriptor(fft->ooplace[0]);
 159
 160
 161     if( status != 0 )
 162     {
 163         gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
 164         gmx_fft_destroy(fft);
 165         return status;
 166     }
 167
 168     fft->ndim     = 1;
 169     fft->nx       = nx;
 170     fft->real_fft = 0;
 171     fft->work     = NULL;
 172
 173     *pfft = fft;
 174     return 0;
 175 }
 176
 177
 178
 179 int
 180 gmx_fft_init_1d_real(gmx_fft_t *        pfft,
 181                      int                nx,
 182                      gmx_fft_flag  flags)
 183 {
 184     gmx_fft_t      fft;
 185     int            d;
 186     int            status;
 187
 188     if(pfft==NULL)
 189     {
 190         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 191         return EINVAL;
 192     }
 193     *pfft = NULL;
 194
 195     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 196     {
 197         return ENOMEM;
 198     }
 199
 200     /* Mark all handles invalid */
 201     for(d=0;d<3;d++)
 202     {
 203         fft->inplace[d] = fft->ooplace[d] = NULL;
 204     }
 205     fft->ooplace[3] = NULL;
 206
 207     status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nx);
 208
 209     if( status == 0 )
 210         status = DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE);
 211
 212     if( status == 0 )
 213         status = DftiCommitDescriptor(fft->inplace[0]);
 214
 215
 216     if( status == 0 )
 217         status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nx);
 218
 219     if( status == 0 )
 220         status = DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE);
 221
 222     if( status == 0 )
 223         status = DftiCommitDescriptor(fft->ooplace[0]);
 224
 225
 226     if(status == DFTI_UNIMPLEMENTED)
 227     {
 228         gmx_fatal(FARGS,
 229                   "The linked Intel MKL version (<6.0?) cannot do real FFTs.");
 230         gmx_fft_destroy(fft);
 231         return status;
 232     }
 233
 234
 235     if( status != 0 )
 236     {
 237         gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
 238         gmx_fft_destroy(fft);
 239         return status;
 240     }
 241
 242     fft->ndim     = 1;
 243     fft->nx       = nx;
 244     fft->real_fft = 1;
 245     fft->work     = NULL;
 246
 247     *pfft = fft;
 248     return 0;
 249 }
 250
 251
 252
 253 int
 254 gmx_fft_init_2d(gmx_fft_t *        pfft,
 255                 int                nx,
 256                 int                ny,
 257                 gmx_fft_flag  flags)
 258 {
 259     gmx_fft_t      fft;
 260     int            d;
 261     int            status;
 262     MKL_LONG       length[2];
 263
 264     if(pfft==NULL)
 265     {
 266         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 267         return EINVAL;
 268     }
 269     *pfft = NULL;
 270
 271     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 272     {
 273         return ENOMEM;
 274     }
 275
 276     /* Mark all handles invalid */
 277     for(d=0;d<3;d++)
 278     {
 279         fft->inplace[d] = fft->ooplace[d] = NULL;
 280     }
 281     fft->ooplace[3] = NULL;
 282
 283     length[0] = nx;
 284     length[1] = ny;
 285
 286     status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,2,length);
 287
 288     if( status == 0 )
 289         status = DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE);
 290
 291     if( status == 0 )
 292         status = DftiCommitDescriptor(fft->inplace[0]);
 293
 294
 295     if( status == 0 )
 296         status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,2,length);
 297
 298     if( status == 0 )
 299         status = DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE);
 300
 301     if( status == 0 )
 302         status = DftiCommitDescriptor(fft->ooplace[0]);
 303
 304
 305     if( status != 0 )
 306     {
 307         gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
 308         gmx_fft_destroy(fft);
 309         return status;
 310     }
 311
 312     fft->ndim     = 2;
 313     fft->nx       = nx;
 314     fft->ny       = ny;
 315     fft->real_fft = 0;
 316     fft->work     = NULL;
 317
 318     *pfft = fft;
 319     return 0;
 320 }
 321
 322
 323
 324 int
 325 gmx_fft_init_2d_real(gmx_fft_t *        pfft,
 326                      int                nx,
 327                      int                ny,
 328                      gmx_fft_flag  flags)
 329 {
 330     gmx_fft_t      fft;
 331     int            d;
 332     int            status;
 333     MKL_LONG       stride[2];
 334     MKL_LONG       nyc;
 335
 336     if(pfft==NULL)
 337     {
 338         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 339         return EINVAL;
 340     }
 341     *pfft = NULL;
 342
 343     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 344     {
 345         return ENOMEM;
 346     }
 347
 348     nyc = (ny/2 + 1);
 349
 350     /* Mark all handles invalid */
 351     for(d=0;d<3;d++)
 352     {
 353         fft->inplace[d] = fft->ooplace[d] = NULL;
 354     }
 355     fft->ooplace[3] = NULL;
 356
 357     /* Roll our own 2D real transform using multiple transforms in MKL,
 358      * since the current MKL versions does not support our storage format,
 359      * and all but the most recent don't even have 2D real FFTs.
 360      */
 361
 362     /* In-place X FFT */
 363     status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
 364
 365     if ( status == 0 )
 366     {
 367         stride[0]  = 0;
 368         stride[1]  = nyc;
 369
 370         status =
 371             (DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE)    ||
 372              DftiSetValue(fft->inplace[0],DFTI_NUMBER_OF_TRANSFORMS,nyc)  ||
 373              DftiSetValue(fft->inplace[0],DFTI_INPUT_DISTANCE,1)          ||
 374              DftiSetValue(fft->inplace[0],DFTI_INPUT_STRIDES,stride)      ||
 375              DftiSetValue(fft->inplace[0],DFTI_OUTPUT_DISTANCE,1)         ||
 376              DftiSetValue(fft->inplace[0],DFTI_OUTPUT_STRIDES,stride));
 377     }
 378
 379     if( status == 0 )
 380         status = DftiCommitDescriptor(fft->inplace[0]);
 381
 382     /* Out-of-place X FFT */
 383     if( status == 0 )
 384         status = DftiCreateDescriptor(&(fft->ooplace[0]),GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
 385
 386     if( status == 0 )
 387     {
 388         stride[0] = 0;
 389         stride[1] = nyc;
 390
 391         status =
 392             (DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE) ||
 393              DftiSetValue(fft->ooplace[0],DFTI_NUMBER_OF_TRANSFORMS,nyc)   ||
 394              DftiSetValue(fft->ooplace[0],DFTI_INPUT_DISTANCE,1)           ||
 395              DftiSetValue(fft->ooplace[0],DFTI_INPUT_STRIDES,stride)       ||
 396              DftiSetValue(fft->ooplace[0],DFTI_OUTPUT_DISTANCE,1)          ||
 397              DftiSetValue(fft->ooplace[0],DFTI_OUTPUT_STRIDES,stride));
 398     }
 399
 400     if( status == 0 )
 401         status = DftiCommitDescriptor(fft->ooplace[0]);
 402
 403
 404     /* In-place Y FFT  */
 405     if( status == 0 )
 406         status = DftiCreateDescriptor(&fft->inplace[1],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)ny);
 407
 408     if( status == 0 )
 409     {
 410         stride[0] = 0;
 411         stride[1] = 1;
 412
 413         status =
 414             (DftiSetValue(fft->inplace[1],DFTI_PLACEMENT,DFTI_INPLACE)             ||
 415              DftiSetValue(fft->inplace[1],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx)  ||
 416              DftiSetValue(fft->inplace[1],DFTI_INPUT_DISTANCE,2*nyc)               ||
 417              DftiSetValue(fft->inplace[1],DFTI_INPUT_STRIDES,stride)               ||
 418              DftiSetValue(fft->inplace[1],DFTI_OUTPUT_DISTANCE,2*nyc)              ||
 419              DftiSetValue(fft->inplace[1],DFTI_OUTPUT_STRIDES,stride)              ||
 420              DftiCommitDescriptor(fft->inplace[1]));
 421     }
 422
 423
 424     /* Out-of-place real-to-complex (affects output distance) Y FFT */
 425     if( status == 0 )
 426         status = DftiCreateDescriptor(&fft->ooplace[1],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)ny);
 427
 428     if( status == 0 )
 429     {
 430         stride[0] = 0;
 431         stride[1] = 1;
 432
 433         status =
 434             (DftiSetValue(fft->ooplace[1],DFTI_PLACEMENT,DFTI_NOT_INPLACE)           ||
 435              DftiSetValue(fft->ooplace[1],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx)    ||
 436              DftiSetValue(fft->ooplace[1],DFTI_INPUT_DISTANCE,(MKL_LONG)ny)          ||
 437              DftiSetValue(fft->ooplace[1],DFTI_INPUT_STRIDES,stride)                 ||
 438              DftiSetValue(fft->ooplace[1],DFTI_OUTPUT_DISTANCE,2*nyc)                ||
 439              DftiSetValue(fft->ooplace[1],DFTI_OUTPUT_STRIDES,stride)                ||
 440              DftiCommitDescriptor(fft->ooplace[1]));
 441     }
 442
 443
 444     /* Out-of-place complex-to-real (affects output distance) Y FFT */
 445     if( status == 0 )
 446         status = DftiCreateDescriptor(&fft->ooplace[2],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)ny);
 447
 448     if( status == 0 )
 449     {
 450         stride[0] = 0;
 451         stride[1] = 1;
 452
 453         status =
 454             (DftiSetValue(fft->ooplace[2],DFTI_PLACEMENT,DFTI_NOT_INPLACE)           ||
 455              DftiSetValue(fft->ooplace[2],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx)    ||
 456              DftiSetValue(fft->ooplace[2],DFTI_INPUT_DISTANCE,2*nyc)                 ||
 457              DftiSetValue(fft->ooplace[2],DFTI_INPUT_STRIDES,stride)                 ||
 458              DftiSetValue(fft->ooplace[2],DFTI_OUTPUT_DISTANCE,(MKL_LONG)ny)         ||
 459              DftiSetValue(fft->ooplace[2],DFTI_OUTPUT_STRIDES,stride)                ||
 460              DftiCommitDescriptor(fft->ooplace[2]));
 461     }
 462
 463
 464     if ( status == 0 )
 465     {
 466         if ((fft->work = (t_complex *)malloc(sizeof(t_complex)*(nx*(ny/2+1)))) == NULL)
 467         {
 468             status = ENOMEM;
 469         }
 470     }
 471
 472     if( status != 0 )
 473     {
 474         gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
 475         gmx_fft_destroy(fft);
 476         return status;
 477     }
 478
 479     fft->ndim     = 2;
 480     fft->nx       = nx;
 481     fft->ny       = ny;
 482     fft->real_fft = 1;
 483
 484     *pfft = fft;
 485     return 0;
 486 }
 487
 488
 489
 490 int
 491 gmx_fft_init_3d(gmx_fft_t *        pfft,
 492                 int                nx,
 493                 int                ny,
 494                 int                nz,
 495                 gmx_fft_flag  flags)
 496 {
 497     gmx_fft_t      fft;
 498     int            d;
 499     MKL_LONG       length[3];
 500     int            status;
 501
 502     if(pfft==NULL)
 503     {
 504         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 505         return EINVAL;
 506     }
 507     *pfft = NULL;
 508
 509     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 510     {
 511         return ENOMEM;
 512     }
 513
 514     /* Mark all handles invalid */
 515     for(d=0;d<3;d++)
 516     {
 517         fft->inplace[d] = fft->ooplace[d] = NULL;
 518     }
 519     fft->ooplace[3] = NULL;
 520
 521     length[0] = nx;
 522     length[1] = ny;
 523     length[2] = nz;
 524
 525     status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,(MKL_LONG)3,length);
 526
 527     if( status == 0 )
 528         status = DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE);
 529
 530     if( status == 0 )
 531         status = DftiCommitDescriptor(fft->inplace[0]);
 532
 533
 534     if( status == 0 )
 535         status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,(MKL_LONG)3,length);
 536
 537     if( status == 0 )
 538         status = DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE);
 539
 540     if( status == 0 )
 541         status = DftiCommitDescriptor(fft->ooplace[0]);
 542
 543
 544     if( status != 0 )
 545     {
 546         gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
 547         gmx_fft_destroy(fft);
 548         return status;
 549     }
 550
 551
 552     fft->ndim     = 3;
 553     fft->nx       = nx;
 554     fft->ny       = ny;
 555     fft->nz       = nz;
 556     fft->real_fft = 0;
 557     fft->work     = NULL;
 558
 559     *pfft = fft;
 560     return 0;
 561 }
 562
 563
 564
 565
 566 int
 567 gmx_fft_init_3d_real(gmx_fft_t *        pfft,
 568                      int                nx,
 569                      int                ny,
 570                      int                nz,
 571                      gmx_fft_flag  flags)
 572 {
 573     gmx_fft_t      fft;
 574     int            d;
 575     int            status;
 576     MKL_LONG       stride[2];
 577     int            nzc;
 578
 579     if(pfft==NULL)
 580     {
 581         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 582         return EINVAL;
 583     }
 584     *pfft = NULL;
 585
 586     nzc = (nz/2 + 1);
 587
 588     if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 589     {
 590         return ENOMEM;
 591     }
 592
 593     /* Mark all handles invalid */
 594     for(d=0;d<3;d++)
 595     {
 596         fft->inplace[d] = fft->ooplace[d] = NULL;
 597     }
 598     fft->ooplace[3] = NULL;
 599
 600     /* Roll our own 3D real transform using multiple transforms in MKL,
 601      * since the current MKL versions does not support our storage format
 602      * or 3D real transforms.
 603      */
 604
 605     /* In-place X FFT.
 606      * ny*nzc complex-to-complex transforms, length nx
 607      * transform distance: 1
 608      * element strides: ny*nzc
 609      */
 610     status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
 611
 612     if ( status == 0)
 613     {
 614         stride[0] = 0;
 615         stride[1] = ny*nzc;
 616
 617         status =
 618         (DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE)                ||
 619          DftiSetValue(fft->inplace[0],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)ny*nzc) ||
 620          DftiSetValue(fft->inplace[0],DFTI_INPUT_DISTANCE,1)                      ||
 621          DftiSetValue(fft->inplace[0],DFTI_INPUT_STRIDES,stride)                  ||
 622          DftiSetValue(fft->inplace[0],DFTI_OUTPUT_DISTANCE,1)                     ||
 623          DftiSetValue(fft->inplace[0],DFTI_OUTPUT_STRIDES,stride)                 ||
 624          DftiCommitDescriptor(fft->inplace[0]));
 625     }
 626
 627     /* Out-of-place X FFT:
 628      * ny*nzc complex-to-complex transforms, length nx
 629      * transform distance: 1
 630      * element strides: ny*nzc
 631      */
 632     if( status == 0 )
 633         status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
 634
 635     if( status == 0 )
 636     {
 637         stride[0] = 0;
 638         stride[1] = ny*nzc;
 639
 640         status =
 641         (DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE)              ||
 642          DftiSetValue(fft->ooplace[0],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)ny*nzc)   ||
 643          DftiSetValue(fft->ooplace[0],DFTI_INPUT_DISTANCE,1)                        ||
 644          DftiSetValue(fft->ooplace[0],DFTI_INPUT_STRIDES,stride)                    ||
 645          DftiSetValue(fft->ooplace[0],DFTI_OUTPUT_DISTANCE,1)                       ||
 646          DftiSetValue(fft->ooplace[0],DFTI_OUTPUT_STRIDES,stride)                   ||
 647          DftiCommitDescriptor(fft->ooplace[0]));
 648     }
 649
 650
 651     /* In-place Y FFT.
 652      * We cannot do all NX*NZC transforms at once, so define a handle to do
 653      * NZC transforms, and then execute it NX times.
 654      * nzc complex-to-complex transforms, length ny
 655      * transform distance: 1
 656      * element strides: nzc
 657      */
 658     if( status == 0 )
 659         status = DftiCreateDescriptor(&fft->inplace[1],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)ny);
 660
 661     if( status == 0 )
 662     {
 663         stride[0] = 0;
 664         stride[1] = nzc;
 665
 666         status =
 667         (DftiSetValue(fft->inplace[1],DFTI_PLACEMENT,DFTI_INPLACE)                ||
 668          DftiSetValue(fft->inplace[1],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nzc)    ||
 669          DftiSetValue(fft->inplace[1],DFTI_INPUT_DISTANCE,1)                      ||
 670          DftiSetValue(fft->inplace[1],DFTI_INPUT_STRIDES,stride)                  ||
 671          DftiSetValue(fft->inplace[1],DFTI_OUTPUT_DISTANCE,1)                     ||
 672          DftiSetValue(fft->inplace[1],DFTI_OUTPUT_STRIDES,stride)                 ||
 673          DftiCommitDescriptor(fft->inplace[1]));
 674     }
 675
 676
 677     /* Out-of-place Y FFT:
 678      * We cannot do all NX*NZC transforms at once, so define a handle to do
 679      * NZC transforms, and then execute it NX times.
 680      * nzc complex-to-complex transforms, length ny
 681      * transform distance: 1
 682      * element strides: nzc
 683      */
 684     if( status == 0 )
 685         status = DftiCreateDescriptor(&fft->ooplace[1],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)ny);
 686
 687     if( status == 0 )
 688     {
 689         stride[0] = 0;
 690         stride[1] = nzc;
 691
 692         status =
 693         (DftiSetValue(fft->ooplace[1],DFTI_PLACEMENT,DFTI_NOT_INPLACE)            ||
 694          DftiSetValue(fft->ooplace[1],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nzc)    ||
 695          DftiSetValue(fft->ooplace[1],DFTI_INPUT_DISTANCE,1)                      ||
 696          DftiSetValue(fft->ooplace[1],DFTI_INPUT_STRIDES,stride)                  ||
 697          DftiSetValue(fft->ooplace[1],DFTI_OUTPUT_DISTANCE,1)                     ||
 698          DftiSetValue(fft->ooplace[1],DFTI_OUTPUT_STRIDES,stride)                 ||
 699          DftiCommitDescriptor(fft->ooplace[1]));
 700     }
 701
 702     /* In-place Z FFT:
 703      * nx*ny real-to-complex transforms, length nz
 704      * transform distance: nzc*2 -> nzc*2
 705      * element strides: 1
 706      */
 707     if( status == 0 )
 708         status = DftiCreateDescriptor(&fft->inplace[2],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nz);
 709
 710     if( status == 0 )
 711     {
 712         stride[0] = 0;
 713         stride[1] = 1;
 714
 715         status =
 716         (DftiSetValue(fft->inplace[2],DFTI_PLACEMENT,DFTI_INPLACE)               ||
 717          DftiSetValue(fft->inplace[2],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx*ny) ||
 718          DftiSetValue(fft->inplace[2],DFTI_INPUT_DISTANCE,(MKL_LONG)nzc*2)       ||
 719          DftiSetValue(fft->inplace[2],DFTI_INPUT_STRIDES,stride)                 ||
 720          DftiSetValue(fft->inplace[2],DFTI_OUTPUT_DISTANCE,(MKL_LONG)nzc*2)      ||
 721          DftiSetValue(fft->inplace[2],DFTI_OUTPUT_STRIDES,stride)                ||
 722          DftiCommitDescriptor(fft->inplace[2]));
 723     }
 724
 725
 726     /* Out-of-place real-to-complex (affects distance) Z FFT:
 727      * nx*ny real-to-complex transforms, length nz
 728      * transform distance: nz -> nzc*2
 729      * element STRIDES: 1
 730      */
 731     if( status == 0 )
 732         status = DftiCreateDescriptor(&fft->ooplace[2],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nz);
 733
 734     if( status == 0 )
 735     {
 736         stride[0] = 0;
 737         stride[1] = 1;
 738
 739         status =
 740         (DftiSetValue(fft->ooplace[2],DFTI_PLACEMENT,DFTI_NOT_INPLACE)           ||
 741          DftiSetValue(fft->ooplace[2],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx*ny) ||
 742          DftiSetValue(fft->ooplace[2],DFTI_INPUT_DISTANCE,(MKL_LONG)nz)          ||
 743          DftiSetValue(fft->ooplace[2],DFTI_INPUT_STRIDES,stride)                 ||
 744          DftiSetValue(fft->ooplace[2],DFTI_OUTPUT_DISTANCE,(MKL_LONG)nzc*2)      ||
 745          DftiSetValue(fft->ooplace[2],DFTI_OUTPUT_STRIDES,stride)                ||
 746          DftiCommitDescriptor(fft->ooplace[2]));
 747     }
 748
 749
 750     /* Out-of-place complex-to-real (affects distance) Z FFT:
 751      * nx*ny real-to-complex transforms, length nz
 752      * transform distance: nzc*2 -> nz
 753      * element STRIDES: 1
 754      */
 755     if( status == 0 )
 756         status = DftiCreateDescriptor(&fft->ooplace[3],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nz);
 757
 758     if( status == 0 )
 759     {
 760         stride[0] = 0;
 761         stride[1] = 1;
 762
 763         status =
 764             (DftiSetValue(fft->ooplace[3],DFTI_PLACEMENT,DFTI_NOT_INPLACE)           ||
 765              DftiSetValue(fft->ooplace[3],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx*ny) ||
 766              DftiSetValue(fft->ooplace[3],DFTI_INPUT_DISTANCE,(MKL_LONG)nzc*2)       ||
 767              DftiSetValue(fft->ooplace[3],DFTI_INPUT_STRIDES,stride)                 ||
 768              DftiSetValue(fft->ooplace[3],DFTI_OUTPUT_DISTANCE,(MKL_LONG)nz)         ||
 769              DftiSetValue(fft->ooplace[3],DFTI_OUTPUT_STRIDES,stride)                ||
 770              DftiCommitDescriptor(fft->ooplace[3]));
 771     }
 772
 773
 774     if ( status == 0 )
 775     {
 776         if ((fft->work = (t_complex *)malloc(sizeof(t_complex)*(nx*ny*(nz/2+1)))) == NULL)
 777         {
 778             status = ENOMEM;
 779         }
 780     }
 781
 782
 783     if( status != 0 )
 784     {
 785         gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
 786         gmx_fft_destroy(fft);
 787         return status;
 788     }
 789
 790
 791     fft->ndim     = 3;
 792     fft->nx       = nx;
 793     fft->ny       = ny;
 794     fft->nz       = nz;
 795     fft->real_fft = 1;
 796
 797     *pfft = fft;
 798     return 0;
 799 }
 800
 801
 802
 803
 804 int
 805 gmx_fft_1d(gmx_fft_t                  fft,
 806            enum gmx_fft_direction     dir,
 807            void *                     in_data,
 808            void *                     out_data)
 809 {
 810     int inplace = (in_data == out_data);
 811     int status = 0;
 812
 813     if( (fft->real_fft == 1) || (fft->ndim != 1) ||
 814         ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
 815     {
 816         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 817         return EINVAL;
 818     }
 819
 820     if(dir==GMX_FFT_FORWARD)
 821     {
 822         if(inplace)
 823         {
 824             status = DftiComputeForward(fft->inplace[0],in_data);
 825         }
 826         else
 827         {
 828             status = DftiComputeForward(fft->ooplace[0],in_data,out_data);
 829         }
 830     }
 831     else
 832     {
 833         if(inplace)
 834         {
 835             status = DftiComputeBackward(fft->inplace[0],in_data);
 836         }
 837         else
 838         {
 839             status = DftiComputeBackward(fft->ooplace[0],in_data,out_data);
 840         }
 841     }
 842
 843     if( status != 0 )
 844     {
 845         gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
 846         status = -1;
 847     }
 848
 849     return status;
 850 }
 851
 852
 853
 854 int
 855 gmx_fft_1d_real(gmx_fft_t                  fft,
 856                 enum gmx_fft_direction     dir,
 857                 void *                     in_data,
 858                 void *                     out_data)
 859 {
 860     int inplace = (in_data == out_data);
 861     int status = 0;
 862
 863     if( (fft->real_fft != 1) || (fft->ndim != 1) ||
 864         ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 865     {
 866         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 867         return EINVAL;
 868     }
 869
 870     if(dir==GMX_FFT_REAL_TO_COMPLEX)
 871     {
 872         if(inplace)
 873         {
 874             status = DftiComputeForward(fft->inplace[0],in_data);
 875         }
 876         else
 877         {
 878             status = DftiComputeForward(fft->ooplace[0],in_data,out_data);
 879         }
 880     }
 881     else
 882     {
 883         if(inplace)
 884         {
 885             status = DftiComputeBackward(fft->inplace[0],in_data);
 886         }
 887         else
 888         {
 889             status = DftiComputeBackward(fft->ooplace[0],in_data,out_data);
 890         }
 891     }
 892
 893     if( status != 0 )
 894     {
 895         gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
 896         status = -1;
 897     }
 898
 899     return status;
 900 }
 901
 902
 903 int
 904 gmx_fft_2d(gmx_fft_t                  fft,
 905            enum gmx_fft_direction     dir,
 906            void *                     in_data,
 907            void *                     out_data)
 908 {
 909     int inplace = (in_data == out_data);
 910     int status = 0;
 911
 912     if( (fft->real_fft == 1) || (fft->ndim != 2) ||
 913         ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
 914     {
 915         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 916         return EINVAL;
 917     }
 918
 919     if(dir==GMX_FFT_FORWARD)
 920     {
 921         if(inplace)
 922         {
 923             status = DftiComputeForward(fft->inplace[0],in_data);
 924         }
 925         else
 926         {
 927             status = DftiComputeForward(fft->ooplace[0],in_data,out_data);
 928         }
 929     }
 930     else
 931     {
 932         if(inplace)
 933         {
 934             status = DftiComputeBackward(fft->inplace[0],in_data);
 935         }
 936         else
 937         {
 938             status = DftiComputeBackward(fft->ooplace[0],in_data,out_data);
 939         }
 940     }
 941
 942     if( status != 0 )
 943     {
 944         gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
 945         status = -1;
 946     }
 947
 948     return status;
 949 }
 950
 951
 952 int
 953 gmx_fft_2d_real(gmx_fft_t                  fft,
 954                 enum gmx_fft_direction     dir,
 955                 void *                     in_data,
 956                 void *                     out_data)
 957 {
 958     int inplace = (in_data == out_data);
 959     int status = 0;
 960
 961     if( (fft->real_fft != 1) || (fft->ndim != 2) ||
 962         ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 963     {
 964         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 965         return EINVAL;
 966     }
 967
 968     if(dir==GMX_FFT_REAL_TO_COMPLEX)
 969     {
 970         if(inplace)
 971         {
 972             /* real-to-complex in Y dimension, in-place */
 973             status = DftiComputeForward(fft->inplace[1],in_data);
 974
 975             /* complex-to-complex in X dimension, in-place */
 976             if ( status == 0 )
 977                 status = DftiComputeForward(fft->inplace[0],in_data);
 978         }
 979         else
 980         {
 981             /* real-to-complex in Y dimension, in_data to out_data */
 982             status = DftiComputeForward(fft->ooplace[1],in_data,out_data);
 983
 984             /* complex-to-complex in X dimension, in-place to out_data */
 985             if ( status == 0 )
 986                 status = DftiComputeForward(fft->inplace[0],out_data);
 987         }
 988     }
 989     else
 990     {
 991         if(inplace)
 992         {
 993             /* complex-to-complex in X dimension, in-place */
 994             status = DftiComputeBackward(fft->inplace[0],in_data);
 995
 996             /* complex-to-real in Y dimension, in-place */
 997             if ( status == 0 )
 998                 status = DftiComputeBackward(fft->inplace[1],in_data);
 999
1000         }
1001         else
1002         {
1003             /* complex-to-complex in X dimension, from in_data to work */
1004             status = DftiComputeBackward(fft->ooplace[0],in_data,fft->work);
1005
1006             /* complex-to-real in Y dimension, from work to out_data */
1007             if ( status == 0 )
1008                 status = DftiComputeBackward(fft->ooplace[1],fft->work,out_data);
1009
1010         }
1011     }
1012
1013     if( status != 0 )
1014     {
1015         gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
1016         status = -1;
1017     }
1018
1019     return status;
1020 }
1021
1022
1023 int
1024 gmx_fft_3d(gmx_fft_t                  fft,
1025            enum gmx_fft_direction     dir,
1026            void *                     in_data,
1027            void *                     out_data)
1028 {
1029     int inplace = (in_data == out_data);
1030     int status = 0;
1031
1032     if( (fft->real_fft == 1) || (fft->ndim != 3) ||
1033         ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
1034     {
1035         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
1036         return EINVAL;
1037     }
1038
1039     if(dir==GMX_FFT_FORWARD)
1040     {
1041         if(inplace)
1042         {
1043             status = DftiComputeForward(fft->inplace[0],in_data);
1044         }
1045         else
1046         {
1047             status = DftiComputeForward(fft->ooplace[0],in_data,out_data);
1048         }
1049     }
1050     else
1051     {
1052         if(inplace)
1053         {
1054             status = DftiComputeBackward(fft->inplace[0],in_data);
1055         }
1056         else
1057         {
1058             status = DftiComputeBackward(fft->ooplace[0],in_data,out_data);
1059         }
1060     }
1061
1062     if( status != 0 )
1063     {
1064         gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
1065         status = -1;
1066     }
1067
1068     return status;
1069 }
1070
1071
1072 int
1073 gmx_fft_3d_real(gmx_fft_t                  fft,
1074                 enum gmx_fft_direction     dir,
1075                 void *                     in_data,
1076                 void *                     out_data)
1077 {
1078     int inplace = (in_data == out_data);
1079     int status = 0;
1080     int i;
1081     int nx,ny,nzc;
1082
1083     nx  = fft->nx;
1084     ny  = fft->ny;
1085     nzc = fft->nz/2 + 1;
1086
1087     if( (fft->real_fft != 1) || (fft->ndim != 3) ||
1088         ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
1089     {
1090         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
1091         return EINVAL;
1092     }
1093
1094     if(dir==GMX_FFT_REAL_TO_COMPLEX)
1095     {
1096         if(inplace)
1097         {
1098             /* real-to-complex in Z dimension, in-place */
1099             status = DftiComputeForward(fft->inplace[2],in_data);
1100
1101             /* complex-to-complex in Y dimension, in-place */
1102             for(i=0;i<nx;i++)
1103             {
1104                 if ( status == 0 )
1105                     status = DftiComputeForward(fft->inplace[1],(t_complex *)in_data+i*ny*nzc);
1106             }
1107
1108             /* complex-to-complex in X dimension, in-place */
1109             if ( status == 0 )
1110                 status = DftiComputeForward(fft->inplace[0],in_data);
1111         }
1112         else
1113         {
1114             /* real-to-complex in Z dimension, from in_data to out_data */
1115             status = DftiComputeForward(fft->ooplace[2],in_data,out_data);
1116
1117             /* complex-to-complex in Y dimension, in-place */
1118             for(i=0;i<nx;i++)
1119             {
1120                 if ( status == 0 )
1121                     status = DftiComputeForward(fft->inplace[1],(t_complex *)out_data+i*ny*nzc);
1122             }
1123
1124             /* complex-to-complex in X dimension, in-place */
1125             if ( status == 0 )
1126                 status = DftiComputeForward(fft->inplace[0],out_data);
1127         }
1128     }
1129     else
1130     {
1131         if(inplace)
1132         {
1133             /* complex-to-complex in X dimension, in-place */
1134             status = DftiComputeBackward(fft->inplace[0],in_data);
1135
1136             /* complex-to-complex in Y dimension, in-place */
1137             for(i=0;i<nx;i++)
1138             {
1139                 if ( status == 0 )
1140                     status = DftiComputeBackward(fft->inplace[1],(t_complex *)in_data+i*ny*nzc);
1141             }
1142
1143             /* complex-to-real in Z dimension, in-place */
1144             if ( status == 0 )
1145                 status = DftiComputeBackward(fft->inplace[2],in_data);
1146         }
1147         else
1148         {
1149             /* complex-to-complex in X dimension, from in_data to work */
1150             status = DftiComputeBackward(fft->ooplace[0],in_data,fft->work);
1151
1152             /* complex-to-complex in Y dimension, in-place */
1153             for(i=0;i<nx;i++)
1154             {
1155                 if ( status == 0 )
1156                     status = DftiComputeBackward(fft->inplace[1],fft->work+i*ny*nzc);
1157             }
1158
1159             /* complex-to-real in Z dimension, work to out_data */
1160             if ( status == 0 )
1161                 status = DftiComputeBackward(fft->ooplace[2],fft->work,out_data);
1162         }
1163     }
1164
1165     if( status != 0 )
1166     {
1167         gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
1168         status = -1;
1169     }
1170
1171     return status;
1172 }
1173
1174
1175
1176 void
1177 gmx_fft_destroy(gmx_fft_t    fft)
1178 {
1179     int d;
1180
1181     if(fft != NULL)
1182     {
1183         for(d=0;d<3;d++)
1184         {
1185             if(fft->inplace[d] != NULL)
1186             {
1187                 DftiFreeDescriptor(&fft->inplace[d]);
1188             }
1189             if(fft->ooplace[d] != NULL)
1190             {
1191                 DftiFreeDescriptor(&fft->ooplace[d]);
1192             }
1193         }
1194         if(fft->ooplace[3] != NULL)
1195         {
1196             DftiFreeDescriptor(&fft->ooplace[3]);
1197         }
1198         free(fft);
1199     }
1200 }
1201
1202 #else
1203 int
1204 gmx_fft_mkl_empty;
1205 #endif /* GMX_FFT_MKL */