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