src/gromacs/fft/fft_mkl.cpp

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