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