src/gromacs/fft/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 #include <errno.h>
  23 #include <stdlib.h>
  24
  25 #include <mkl_dfti.h>
  26 #include <mkl_service.h>
  27
  28 #include "gromacs/fft/fft.h"
  29 #include "gmx_fatal.h"
  30
  31
  32 /* For MKL version (<10.0), we should define MKL_LONG. */
  33 #ifndef MKL_LONG
  34 #define MKL_LONG long int
  35 #endif
  36
  37
  38 #ifdef GMX_DOUBLE
  39 #define GMX_DFTI_PREC  DFTI_DOUBLE
  40 #else
  41 #define GMX_DFTI_PREC  DFTI_SINGLE
  42 #endif
  43
  44 /*! \internal \brief
  45  * Contents of the Intel MKL FFT fft datatype.
  46  *
  47  * Note that this is one of several possible implementations of gmx_fft_t.
  48  *
  49  *  The MKL _API_ supports 1D,2D, and 3D transforms, including real-to-complex.
  50  *  Unfortunately the actual library implementation does not support 3D real
  51  *  transforms as of version 7.2, and versions before 7.0 don't support 2D real
  52  *  either. In addition, the multi-dimensional storage format for real data
  53  *  is not compatible with our padding.
  54  *
  55  *  To work around this we roll our own 2D and 3D real-to-complex transforms,
  56  *  using separate X/Y/Z handles defined to perform (ny*nz), (nx*nz), and
  57  *  (nx*ny) transforms at once when necessary. To perform strided multiple
  58  *  transforms out-of-place (i.e., without padding in the last dimension)
  59  *  on the fly we also need to separate the forward and backward
  60  *  handles for real-to-complex/complex-to-real data permutation.
  61  *
  62  *  This makes it necessary to define 3 handles for in-place FFTs, and 4 for
  63  *  the out-of-place transforms. Still, whenever possible we try to use
  64  *  a single 3D-transform handle instead.
  65  *
  66  *  So, the handles are enumerated as follows:
  67  *
  68  *  1D FFT (real too):    Index 0 is the handle for the entire FFT
  69  *  2D complex FFT:       Index 0 is the handle for the entire FFT
  70  *  3D complex FFT:       Index 0 is the handle for the entire FFT
  71  *  2D, inplace real FFT: 0=FFTx, 1=FFTy handle
  72  *  2D, ooplace real FFT: 0=FFTx, 1=real-to-complex FFTy, 2=complex-to-real FFTy
  73  *  3D, inplace real FFT: 0=FFTx, 1=FFTy, 2=FFTz handle
  74  *  3D, ooplace real FFT: 0=FFTx, 1=FFTy, 2=r2c FFTz, 3=c2r FFTz
  75  *
  76  *  Intel people reading this: Learn from FFTW what a good interface looks like :-)
  77  */
  78 #ifdef DOXYGEN
  79 struct gmx_fft_mkl
  80 #else
  81 struct gmx_fft
  82 #endif
  83 {
  84     int                ndim;          /**< Number of dimensions in FFT  */
  85     int                nx;            /**< Length of X transform        */
  86     int                ny;            /**< Length of Y transform        */
  87     int                nz;            /**< Length of Z transform        */
  88     int                real_fft;      /**< 1 if real FFT, otherwise 0   */
  89     DFTI_DESCRIPTOR *  inplace[3];    /**< in-place FFT                 */
  90     DFTI_DESCRIPTOR *  ooplace[4];    /**< out-of-place FFT             */
  91     t_complex     *    work;          /**< Enable out-of-place c2r FFT  */
  92 };
  93
  94
  95
  96 int
  97 gmx_fft_init_1d(gmx_fft_t *        pfft,
  98                 int                nx,
  99                 gmx_fft_flag       flags)
 100 {
 101     gmx_fft_t      fft;
 102     int            d;
 103     int            status;
 104
 105     if (pfft == NULL)
 106     {
 107         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 108         return EINVAL;
 109     }
 110     *pfft = NULL;
 111
 112     if ( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 113     {
 114         return ENOMEM;
 115     }
 116
 117     /* Mark all handles invalid */
 118     for (d = 0; d < 3; d++)
 119     {
 120         fft->inplace[d] = fft->ooplace[d] = NULL;
 121     }
 122     fft->ooplace[3] = NULL;
 123
 124
 125     status = DftiCreateDescriptor(&fft->inplace[0], GMX_DFTI_PREC, DFTI_COMPLEX, 1, (MKL_LONG)nx);
 126
 127     if (status == 0)
 128     {
 129         status = DftiSetValue(fft->inplace[0], DFTI_PLACEMENT, DFTI_INPLACE);
 130     }
 131
 132     if (status == 0)
 133     {
 134         status = DftiCommitDescriptor(fft->inplace[0]);
 135     }
 136
 137
 138     if (status == 0)
 139     {
 140         status = DftiCreateDescriptor(&fft->ooplace[0], GMX_DFTI_PREC, DFTI_COMPLEX, 1, (MKL_LONG)nx);
 141     }
 142
 143     if (status == 0)
 144     {
 145         DftiSetValue(fft->ooplace[0], DFTI_PLACEMENT, DFTI_NOT_INPLACE);
 146     }
 147
 148     if (status == 0)
 149     {
 150         DftiCommitDescriptor(fft->ooplace[0]);
 151     }
 152
 153
 154     if (status != 0)
 155     {
 156         gmx_fatal(FARGS, "Error initializing Intel MKL FFT; status=%d", status);
 157         gmx_fft_destroy(fft);
 158         return status;
 159     }
 160
 161     fft->ndim     = 1;
 162     fft->nx       = nx;
 163     fft->real_fft = 0;
 164     fft->work     = NULL;
 165
 166     *pfft = fft;
 167     return 0;
 168 }
 169
 170
 171
 172 int
 173 gmx_fft_init_1d_real(gmx_fft_t *        pfft,
 174                      int                nx,
 175                      gmx_fft_flag       flags)
 176 {
 177     gmx_fft_t      fft;
 178     int            d;
 179     int            status;
 180
 181     if (pfft == NULL)
 182     {
 183         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 184         return EINVAL;
 185     }
 186     *pfft = NULL;
 187
 188     if ( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 189     {
 190         return ENOMEM;
 191     }
 192
 193     /* Mark all handles invalid */
 194     for (d = 0; d < 3; d++)
 195     {
 196         fft->inplace[d] = fft->ooplace[d] = NULL;
 197     }
 198     fft->ooplace[3] = NULL;
 199
 200     status = DftiCreateDescriptor(&fft->inplace[0], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)nx);
 201
 202     if (status == 0)
 203     {
 204         status = DftiSetValue(fft->inplace[0], DFTI_PLACEMENT, DFTI_INPLACE);
 205     }
 206
 207     if (status == 0)
 208     {
 209         status = DftiCommitDescriptor(fft->inplace[0]);
 210     }
 211
 212
 213     if (status == 0)
 214     {
 215         status = DftiCreateDescriptor(&fft->ooplace[0], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)nx);
 216     }
 217
 218     if (status == 0)
 219     {
 220         status = DftiSetValue(fft->ooplace[0], DFTI_PLACEMENT, DFTI_NOT_INPLACE);
 221     }
 222
 223     if (status == 0)
 224     {
 225         status = DftiCommitDescriptor(fft->ooplace[0]);
 226     }
 227
 228
 229     if (status == DFTI_UNIMPLEMENTED)
 230     {
 231         gmx_fatal(FARGS,
 232                   "The linked Intel MKL version (<6.0?) cannot do real FFTs.");
 233         gmx_fft_destroy(fft);
 234         return status;
 235     }
 236
 237
 238     if (status != 0)
 239     {
 240         gmx_fatal(FARGS, "Error initializing Intel MKL FFT; status=%d", status);
 241         gmx_fft_destroy(fft);
 242         return status;
 243     }
 244
 245     fft->ndim     = 1;
 246     fft->nx       = nx;
 247     fft->real_fft = 1;
 248     fft->work     = NULL;
 249
 250     *pfft = fft;
 251     return 0;
 252 }
 253
 254
 255
 256 int
 257 gmx_fft_init_2d_real(gmx_fft_t *        pfft,
 258                      int                nx,
 259                      int                ny,
 260                      gmx_fft_flag       flags)
 261 {
 262     gmx_fft_t      fft;
 263     int            d;
 264     int            status;
 265     MKL_LONG       stride[2];
 266     MKL_LONG       nyc;
 267
 268     if (pfft == NULL)
 269     {
 270         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 271         return EINVAL;
 272     }
 273     *pfft = NULL;
 274
 275     if ( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
 276     {
 277         return ENOMEM;
 278     }
 279
 280     nyc = (ny/2 + 1);
 281
 282     /* Mark all handles invalid */
 283     for (d = 0; d < 3; d++)
 284     {
 285         fft->inplace[d] = fft->ooplace[d] = NULL;
 286     }
 287     fft->ooplace[3] = NULL;
 288
 289     /* Roll our own 2D real transform using multiple transforms in MKL,
 290      * since the current MKL versions does not support our storage format,
 291      * and all but the most recent don't even have 2D real FFTs.
 292      */
 293
 294     /* In-place X FFT */
 295     status = DftiCreateDescriptor(&fft->inplace[0], GMX_DFTI_PREC, DFTI_COMPLEX, 1, (MKL_LONG)nx);
 296
 297     if (status == 0)
 298     {
 299         stride[0]  = 0;
 300         stride[1]  = nyc;
 301
 302         status =
 303             (DftiSetValue(fft->inplace[0], DFTI_PLACEMENT, DFTI_INPLACE)    ||
 304              DftiSetValue(fft->inplace[0], DFTI_NUMBER_OF_TRANSFORMS, nyc)  ||
 305              DftiSetValue(fft->inplace[0], DFTI_INPUT_DISTANCE, 1)          ||
 306              DftiSetValue(fft->inplace[0], DFTI_INPUT_STRIDES, stride)      ||
 307              DftiSetValue(fft->inplace[0], DFTI_OUTPUT_DISTANCE, 1)         ||
 308              DftiSetValue(fft->inplace[0], DFTI_OUTPUT_STRIDES, stride));
 309     }
 310
 311     if (status == 0)
 312     {
 313         status = DftiCommitDescriptor(fft->inplace[0]);
 314     }
 315
 316     /* Out-of-place X FFT */
 317     if (status == 0)
 318     {
 319         status = DftiCreateDescriptor(&(fft->ooplace[0]), GMX_DFTI_PREC, DFTI_COMPLEX, 1, (MKL_LONG)nx);
 320     }
 321
 322     if (status == 0)
 323     {
 324         stride[0] = 0;
 325         stride[1] = nyc;
 326
 327         status =
 328             (DftiSetValue(fft->ooplace[0], DFTI_PLACEMENT, DFTI_NOT_INPLACE) ||
 329              DftiSetValue(fft->ooplace[0], DFTI_NUMBER_OF_TRANSFORMS, nyc)   ||
 330              DftiSetValue(fft->ooplace[0], DFTI_INPUT_DISTANCE, 1)           ||
 331              DftiSetValue(fft->ooplace[0], DFTI_INPUT_STRIDES, stride)       ||
 332              DftiSetValue(fft->ooplace[0], DFTI_OUTPUT_DISTANCE, 1)          ||
 333              DftiSetValue(fft->ooplace[0], DFTI_OUTPUT_STRIDES, stride));
 334     }
 335
 336     if (status == 0)
 337     {
 338         status = DftiCommitDescriptor(fft->ooplace[0]);
 339     }
 340
 341
 342     /* In-place Y FFT  */
 343     if (status == 0)
 344     {
 345         status = DftiCreateDescriptor(&fft->inplace[1], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)ny);
 346     }
 347
 348     if (status == 0)
 349     {
 350         stride[0] = 0;
 351         stride[1] = 1;
 352
 353         status =
 354             (DftiSetValue(fft->inplace[1], DFTI_PLACEMENT, DFTI_INPLACE)             ||
 355              DftiSetValue(fft->inplace[1], DFTI_NUMBER_OF_TRANSFORMS, (MKL_LONG)nx)  ||
 356              DftiSetValue(fft->inplace[1], DFTI_INPUT_DISTANCE, 2*nyc)               ||
 357              DftiSetValue(fft->inplace[1], DFTI_INPUT_STRIDES, stride)               ||
 358              DftiSetValue(fft->inplace[1], DFTI_OUTPUT_DISTANCE, 2*nyc)              ||
 359              DftiSetValue(fft->inplace[1], DFTI_OUTPUT_STRIDES, stride)              ||
 360              DftiCommitDescriptor(fft->inplace[1]));
 361     }
 362
 363
 364     /* Out-of-place real-to-complex (affects output distance) Y FFT */
 365     if (status == 0)
 366     {
 367         status = DftiCreateDescriptor(&fft->ooplace[1], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)ny);
 368     }
 369
 370     if (status == 0)
 371     {
 372         stride[0] = 0;
 373         stride[1] = 1;
 374
 375         status =
 376             (DftiSetValue(fft->ooplace[1], DFTI_PLACEMENT, DFTI_NOT_INPLACE)           ||
 377              DftiSetValue(fft->ooplace[1], DFTI_NUMBER_OF_TRANSFORMS, (MKL_LONG)nx)    ||
 378              DftiSetValue(fft->ooplace[1], DFTI_INPUT_DISTANCE, (MKL_LONG)ny)          ||
 379              DftiSetValue(fft->ooplace[1], DFTI_INPUT_STRIDES, stride)                 ||
 380              DftiSetValue(fft->ooplace[1], DFTI_OUTPUT_DISTANCE, 2*nyc)                ||
 381              DftiSetValue(fft->ooplace[1], DFTI_OUTPUT_STRIDES, stride)                ||
 382              DftiCommitDescriptor(fft->ooplace[1]));
 383     }
 384
 385
 386     /* Out-of-place complex-to-real (affects output distance) Y FFT */
 387     if (status == 0)
 388     {
 389         status = DftiCreateDescriptor(&fft->ooplace[2], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)ny);
 390     }
 391
 392     if (status == 0)
 393     {
 394         stride[0] = 0;
 395         stride[1] = 1;
 396
 397         status =
 398             (DftiSetValue(fft->ooplace[2], DFTI_PLACEMENT, DFTI_NOT_INPLACE)           ||
 399              DftiSetValue(fft->ooplace[2], DFTI_NUMBER_OF_TRANSFORMS, (MKL_LONG)nx)    ||
 400              DftiSetValue(fft->ooplace[2], DFTI_INPUT_DISTANCE, 2*nyc)                 ||
 401              DftiSetValue(fft->ooplace[2], DFTI_INPUT_STRIDES, stride)                 ||
 402              DftiSetValue(fft->ooplace[2], DFTI_OUTPUT_DISTANCE, (MKL_LONG)ny)         ||
 403              DftiSetValue(fft->ooplace[2], DFTI_OUTPUT_STRIDES, stride)                ||
 404              DftiCommitDescriptor(fft->ooplace[2]));
 405     }
 406
 407
 408     if (status == 0)
 409     {
 410         if ((fft->work = (t_complex *)malloc(sizeof(t_complex)*(nx*(ny/2+1)))) == NULL)
 411         {
 412             status = ENOMEM;
 413         }
 414     }
 415
 416     if (status != 0)
 417     {
 418         gmx_fatal(FARGS, "Error initializing Intel MKL FFT; status=%d", status);
 419         gmx_fft_destroy(fft);
 420         return status;
 421     }
 422
 423     fft->ndim     = 2;
 424     fft->nx       = nx;
 425     fft->ny       = ny;
 426     fft->real_fft = 1;
 427
 428     *pfft = fft;
 429     return 0;
 430 }
 431
 432 int
 433 gmx_fft_1d(gmx_fft_t                  fft,
 434            enum gmx_fft_direction     dir,
 435            void *                     in_data,
 436            void *                     out_data)
 437 {
 438     int inplace = (in_data == out_data);
 439     int status  = 0;
 440
 441     if ( (fft->real_fft == 1) || (fft->ndim != 1) ||
 442          ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
 443     {
 444         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 445         return EINVAL;
 446     }
 447
 448     if (dir == GMX_FFT_FORWARD)
 449     {
 450         if (inplace)
 451         {
 452             status = DftiComputeForward(fft->inplace[0], in_data);
 453         }
 454         else
 455         {
 456             status = DftiComputeForward(fft->ooplace[0], in_data, out_data);
 457         }
 458     }
 459     else
 460     {
 461         if (inplace)
 462         {
 463             status = DftiComputeBackward(fft->inplace[0], in_data);
 464         }
 465         else
 466         {
 467             status = DftiComputeBackward(fft->ooplace[0], in_data, out_data);
 468         }
 469     }
 470
 471     if (status != 0)
 472     {
 473         gmx_fatal(FARGS, "Error executing Intel MKL FFT.");
 474         status = -1;
 475     }
 476
 477     return status;
 478 }
 479
 480
 481
 482 int
 483 gmx_fft_1d_real(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 status  = 0;
 490
 491     if ( (fft->real_fft != 1) || (fft->ndim != 1) ||
 492          ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 493     {
 494         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 495         return EINVAL;
 496     }
 497
 498     if (dir == GMX_FFT_REAL_TO_COMPLEX)
 499     {
 500         if (inplace)
 501         {
 502             status = DftiComputeForward(fft->inplace[0], in_data);
 503         }
 504         else
 505         {
 506             status = DftiComputeForward(fft->ooplace[0], in_data, out_data);
 507         }
 508     }
 509     else
 510     {
 511         if (inplace)
 512         {
 513             status = DftiComputeBackward(fft->inplace[0], in_data);
 514         }
 515         else
 516         {
 517             status = DftiComputeBackward(fft->ooplace[0], in_data, out_data);
 518         }
 519     }
 520
 521     if (status != 0)
 522     {
 523         gmx_fatal(FARGS, "Error executing Intel MKL FFT.");
 524         status = -1;
 525     }
 526
 527     return status;
 528 }
 529
 530
 531 int
 532 gmx_fft_2d_real(gmx_fft_t                  fft,
 533                 enum gmx_fft_direction     dir,
 534                 void *                     in_data,
 535                 void *                     out_data)
 536 {
 537     int inplace = (in_data == out_data);
 538     int status  = 0;
 539
 540     if ( (fft->real_fft != 1) || (fft->ndim != 2) ||
 541          ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 542     {
 543         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 544         return EINVAL;
 545     }
 546
 547     if (dir == GMX_FFT_REAL_TO_COMPLEX)
 548     {
 549         if (inplace)
 550         {
 551             /* real-to-complex in Y dimension, in-place */
 552             status = DftiComputeForward(fft->inplace[1], in_data);
 553
 554             /* complex-to-complex in X dimension, in-place */
 555             if (status == 0)
 556             {
 557                 status = DftiComputeForward(fft->inplace[0], in_data);
 558             }
 559         }
 560         else
 561         {
 562             /* real-to-complex in Y dimension, in_data to out_data */
 563             status = DftiComputeForward(fft->ooplace[1], in_data, out_data);
 564
 565             /* complex-to-complex in X dimension, in-place to out_data */
 566             if (status == 0)
 567             {
 568                 status = DftiComputeForward(fft->inplace[0], out_data);
 569             }
 570         }
 571     }
 572     else
 573     {
 574         /* prior implementation was incorrect. See fft.cpp unit test */
 575         gmx_incons("Complex -> Real is not supported by MKL.");
 576     }
 577
 578     if (status != 0)
 579     {
 580         gmx_fatal(FARGS, "Error executing Intel MKL FFT.");
 581         status = -1;
 582     }
 583
 584     return status;
 585 }
 586
 587 void
 588 gmx_fft_destroy(gmx_fft_t    fft)
 589 {
 590     int d;
 591
 592     if (fft != NULL)
 593     {
 594         for (d = 0; d < 3; d++)
 595         {
 596             if (fft->inplace[d] != NULL)
 597             {
 598                 DftiFreeDescriptor(&fft->inplace[d]);
 599             }
 600             if (fft->ooplace[d] != NULL)
 601             {
 602                 DftiFreeDescriptor(&fft->ooplace[d]);
 603             }
 604         }
 605         if (fft->ooplace[3] != NULL)
 606         {
 607             DftiFreeDescriptor(&fft->ooplace[3]);
 608         }
 609         if (fft->work != NULL)
 610         {
 611             free(fft->work);
 612         }
 613         free(fft);
 614     }
 615 }
 616
 617 void gmx_fft_cleanup()
 618 {
 619     mkl_free_buffers();
 620 }
 621
 622 const char *gmx_fft_get_version_info()
 623 {
 624     return "Intel MKL";
 625 }