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