2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
42 #include <mkl_service.h>
44 #include "gromacs/fft/fft.h"
45 #include "gromacs/utility/fatalerror.h"
48 /* For MKL version (<10.0), we should define MKL_LONG. */
50 # define MKL_LONG long int
55 # define GMX_DFTI_PREC DFTI_DOUBLE
57 # define GMX_DFTI_PREC DFTI_SINGLE
62 * Contents of the Intel MKL FFT fft datatype.
64 * Note that this is one of several possible implementations of gmx_fft_t.
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.
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.
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.
83 * So, the handles are enumerated as follows:
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
93 * Intel people reading this: Learn from FFTW what a good interface looks like :-)
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 */
112 int gmx_fft_init_1d(gmx_fft_t* pfft, int nx, gmx_fft_flag gmx_unused flags)
120 gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
125 if ((fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
130 /* Mark all handles invalid */
131 for (d = 0; d < 3; d++)
133 fft->inplace[d] = fft->ooplace[d] = NULL;
135 fft->ooplace[3] = NULL;
138 status = DftiCreateDescriptor(&fft->inplace[0], GMX_DFTI_PREC, DFTI_COMPLEX, 1, (MKL_LONG)nx);
142 status = DftiSetValue(fft->inplace[0], DFTI_PLACEMENT, DFTI_INPLACE);
147 status = DftiCommitDescriptor(fft->inplace[0]);
153 status = DftiCreateDescriptor(&fft->ooplace[0], GMX_DFTI_PREC, DFTI_COMPLEX, 1, (MKL_LONG)nx);
158 DftiSetValue(fft->ooplace[0], DFTI_PLACEMENT, DFTI_NOT_INPLACE);
163 DftiCommitDescriptor(fft->ooplace[0]);
169 gmx_fatal(FARGS, "Error initializing Intel MKL FFT; status=%d", status);
170 gmx_fft_destroy(fft);
184 int gmx_fft_init_1d_real(gmx_fft_t* pfft, int nx, gmx_fft_flag gmx_unused flags)
192 gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
197 if ((fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
202 /* Mark all handles invalid */
203 for (d = 0; d < 3; d++)
205 fft->inplace[d] = fft->ooplace[d] = NULL;
207 fft->ooplace[3] = NULL;
209 status = DftiCreateDescriptor(&fft->inplace[0], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)nx);
213 status = DftiSetValue(fft->inplace[0], DFTI_PLACEMENT, DFTI_INPLACE);
218 status = DftiCommitDescriptor(fft->inplace[0]);
224 status = DftiCreateDescriptor(&fft->ooplace[0], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)nx);
229 status = DftiSetValue(fft->ooplace[0], DFTI_PLACEMENT, DFTI_NOT_INPLACE);
234 status = DftiCommitDescriptor(fft->ooplace[0]);
238 if (status == DFTI_UNIMPLEMENTED)
240 gmx_fatal(FARGS, "The linked Intel MKL version (<6.0?) cannot do real FFTs.");
241 gmx_fft_destroy(fft);
248 gmx_fatal(FARGS, "Error initializing Intel MKL FFT; status=%d", status);
249 gmx_fft_destroy(fft);
263 int gmx_fft_init_2d_real(gmx_fft_t* pfft, int nx, int ny, gmx_fft_flag gmx_unused flags)
273 gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
278 if ((fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
285 /* Mark all handles invalid */
286 for (d = 0; d < 3; d++)
288 fft->inplace[d] = fft->ooplace[d] = NULL;
290 fft->ooplace[3] = NULL;
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.
298 status = DftiCreateDescriptor(&fft->inplace[0], GMX_DFTI_PREC, DFTI_COMPLEX, 1, (MKL_LONG)nx);
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));
315 status = DftiCommitDescriptor(fft->inplace[0]);
318 /* Out-of-place X FFT */
321 status = DftiCreateDescriptor(&(fft->ooplace[0]), GMX_DFTI_PREC, DFTI_COMPLEX, 1, (MKL_LONG)nx);
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));
339 status = DftiCommitDescriptor(fft->ooplace[0]);
346 status = DftiCreateDescriptor(&fft->inplace[1], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)ny);
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]));
364 /* Out-of-place real-to-complex (affects output distance) Y FFT */
367 status = DftiCreateDescriptor(&fft->ooplace[1], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)ny);
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]));
385 /* Out-of-place complex-to-real (affects output distance) Y FFT */
388 status = DftiCreateDescriptor(&fft->ooplace[2], GMX_DFTI_PREC, DFTI_REAL, 1, (MKL_LONG)ny);
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]));
408 void* memory = malloc(sizeof(t_complex) * (nx * (ny / 2 + 1)));
409 if (nullptr == memory)
413 fft->work = static_cast<t_complex*>(memory);
418 gmx_fatal(FARGS, "Error initializing Intel MKL FFT; status=%d", status);
419 gmx_fft_destroy(fft);
432 int gmx_fft_1d(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
434 int inplace = (in_data == out_data);
437 if ((fft->real_fft == 1) || (fft->ndim != 1) || ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)))
439 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
443 if (dir == GMX_FFT_FORWARD)
447 status = DftiComputeForward(fft->inplace[0], in_data);
451 status = DftiComputeForward(fft->ooplace[0], in_data, out_data);
458 status = DftiComputeBackward(fft->inplace[0], in_data);
462 status = DftiComputeBackward(fft->ooplace[0], in_data, out_data);
468 gmx_fatal(FARGS, "Error executing Intel MKL FFT.");
476 int gmx_fft_1d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
478 int inplace = (in_data == out_data);
481 if ((fft->real_fft != 1) || (fft->ndim != 1)
482 || ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)))
484 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
488 if (dir == GMX_FFT_REAL_TO_COMPLEX)
492 status = DftiComputeForward(fft->inplace[0], in_data);
496 status = DftiComputeForward(fft->ooplace[0], in_data, out_data);
503 status = DftiComputeBackward(fft->inplace[0], in_data);
507 status = DftiComputeBackward(fft->ooplace[0], in_data, out_data);
513 gmx_fatal(FARGS, "Error executing Intel MKL FFT.");
521 int gmx_fft_2d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
523 int inplace = (in_data == out_data);
526 if ((fft->real_fft != 1) || (fft->ndim != 2)
527 || ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)))
529 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
533 if (dir == GMX_FFT_REAL_TO_COMPLEX)
537 /* real-to-complex in Y dimension, in-place */
538 status = DftiComputeForward(fft->inplace[1], in_data);
540 /* complex-to-complex in X dimension, in-place */
543 status = DftiComputeForward(fft->inplace[0], in_data);
548 /* real-to-complex in Y dimension, in_data to out_data */
549 status = DftiComputeForward(fft->ooplace[1], in_data, out_data);
551 /* complex-to-complex in X dimension, in-place to out_data */
554 status = DftiComputeForward(fft->inplace[0], out_data);
560 /* prior implementation was incorrect. See fft.cpp unit test */
561 gmx_incons("Complex -> Real is not supported by MKL.");
566 gmx_fatal(FARGS, "Error executing Intel MKL FFT.");
573 void gmx_fft_destroy(gmx_fft_t fft)
579 for (d = 0; d < 3; d++)
581 if (fft->inplace[d] != NULL)
583 DftiFreeDescriptor(&fft->inplace[d]);
585 if (fft->ooplace[d] != NULL)
587 DftiFreeDescriptor(&fft->ooplace[d]);
590 if (fft->ooplace[3] != NULL)
592 DftiFreeDescriptor(&fft->ooplace[3]);
594 if (fft->work != NULL)
602 void gmx_fft_cleanup()