src/gromacs/fft/fft_mkl.c

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