2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2016,2017,2018,2019, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Tests utilities for fft calculations.
39 * Current reference data is generated in double precision using the Reference
40 * build type, except for the compiler (Apple Clang).
42 * \author Roland Schulz <roland@utk.edu>
47 #include "gromacs/fft/fft.h"
52 #include <gtest/gtest.h>
54 #include "gromacs/fft/parallel_3dfft.h"
55 #include "gromacs/utility/stringutil.h"
57 #include "testutils/refdata.h"
58 #include "testutils/testasserts.h"
63 /*! \brief Input data for FFT tests.
65 * TODO If we require compilers that all support C++11 user literals,
66 * then this array could be of type real, initialized with e.g. -3.5_r
67 * that does not suffer from implicit narrowing with brace
68 * initializers, and we would not have to do so much useless copying
69 * during the unit tests below.
71 const double inputdata[] = {
72 // print ",\n".join([",".join(["%4s"%(random.randint(-99,99)/10.,) for i in range(25)]) for j in range(20)])
73 -3.5, 6.3, 1.2, 0.3, 1.1, -5.7, 5.8, -1.9, -6.3, -1.4, 7.4, 2.4, -9.9, -7.2, 5.4, 6.1,
74 -1.9, -7.6, 1.4, -3.5, 0.7, 5.6, -4.2, -1.1, -4.4, -6.3, -7.2, 4.6, -3.0, -0.9, 7.2, 2.5,
75 -3.6, 6.1, -3.2, -2.1, 6.5, -0.4, -9.0, 2.3, 8.4, 4.0, -5.2, -9.0, 4.7, -3.7, -2.0, -9.5,
76 -3.9, -3.6, 7.1, 0.8, -0.6, 5.2, -9.3, -4.5, 5.9, 2.2, -5.8, 5.0, 1.2, -0.1, 2.2, 0.2,
77 -7.7, 1.9, -8.4, 4.4, 2.3, -2.9, 6.7, 2.7, 5.8, -3.6, 8.9, 8.9, 4.3, 9.1, 9.3, -8.7,
78 4.1, 9.6, -6.2, 6.6, -9.3, 8.2, 4.5, 6.2, 9.4, -8.0, -6.8, -3.3, 7.2, 1.7, 0.6, -4.9,
79 9.8, 1.3, 3.2, -0.2, 9.9, 4.4, -9.9, -7.2, 4.4, 4.7, 7.2, -0.3, 0.3, -2.1, 8.4, -2.1,
80 -6.1, 4.1, -5.9, -2.2, -3.8, 5.2, -8.2, -7.8, -8.8, 6.7, -9.5, -4.2, 0.8, 8.3, 5.2, -9.0,
81 8.7, 9.8, -9.9, -7.8, -8.3, 9.0, -2.8, -9.2, -9.6, 8.4, 2.5, 6.0, -0.4, 1.3, -0.5, 9.1,
82 -9.5, -0.8, 1.9, -6.2, 4.3, -3.8, 8.6, -1.9, -2.1, -0.4, -7.1, -3.7, 9.1, -6.4, -0.6, 2.5,
83 8.0, -5.2, -9.8, -4.3, 4.5, 1.7, 9.3, 9.2, 1.0, 5.3, -4.5, 6.4, -6.6, 3.1, -6.8, 2.1,
84 2.0, 7.3, 8.6, 5.0, 5.2, 0.4, -7.1, 4.5, -9.2, -9.1, 0.2, -6.3, -1.1, -9.6, 7.4, -3.7,
85 -5.5, 2.6, -3.5, -0.7, 9.0, 9.8, -8.0, 3.6, 3.0, -2.2, -2.8, 0.8, 9.0, 2.8, 7.7, -0.7,
86 -5.0, -1.8, -2.3, -0.4, -6.2, -9.1, -9.2, 0.5, 5.7, -3.9, 2.1, 0.6, 0.4, 9.1, 7.4, 7.1,
87 -2.5, 7.3, 7.8, -4.3, 6.3, -0.8, -3.8, -1.5, 6.6, 2.3, 3.9, -4.6, 5.8, -7.4, 5.9, 2.8,
88 4.7, 3.9, -5.4, 9.1, -1.6, -1.9, -4.2, -2.6, 0.6, -5.1, 1.8, 5.2, 4.0, -6.2, 6.5, -9.1,
89 0.5, 2.1, 7.1, -8.6, 7.6, -9.7, -4.6, -5.7, 6.1, -1.8, -7.3, 9.4, 8.0, -2.6, -1.8, 5.7,
90 9.3, -7.9, 7.4, 6.3, 2.0, 9.6, -4.5, -6.2, 6.1, 2.3, 0.8, 5.9, -2.8, -3.5, -1.5, 6.0,
91 -4.9, 3.5, 7.7, -4.2, -9.7, 2.4, 8.1, 5.9, 3.4, -7.5, 7.5, 2.6, 4.7, 2.7, 2.2, 2.6,
92 6.2, 7.5, 0.2, -6.4, -2.8, -0.5, -0.3, 0.4, 1.2, 3.5, -4.0, -0.5, 9.3, -7.2, 8.5, -5.5,
93 -1.7, -5.3, 0.3, 3.9, -3.6, -3.6, 4.7, -8.1, 1.4, 4.0, 1.3, -4.3, -8.8, -7.3, 6.3, -7.5,
94 -9.0, 9.1, 4.5, -1.9, 1.9, 9.9, -1.7, -9.1, -5.1, 8.5, -9.3, 2.1, -5.8, -3.6, -0.8, -0.9,
95 -3.3, -2.7, 7.0, -7.2, -5.0, 7.4, -1.4, 0.0, -4.5, -9.7, 0.7, -1.0, -9.1, -5.3, 4.3, 3.4,
96 -6.6, 9.8, -1.1, 8.9, 5.0, 2.9, 0.2, -2.9, 0.8, 6.7, -0.6, 0.6, 4.1, 5.3, -1.7, -0.3,
97 4.2, 3.7, -8.3, 4.0, 1.3, 6.3, 0.2, 1.3, -1.1, -3.5, 2.8, -7.7, 6.2, -4.9, -9.9, 9.6,
98 3.0, -9.2, -8.0, -3.9, 7.9, -6.1, 6.0, 5.9, 9.6, 1.2, 6.2, 3.6, 2.1, 5.8, 9.2, -8.8,
99 8.8, -3.3, -9.2, 4.6, 1.8, 4.6, 2.9, -2.7, 4.2, 7.3, -0.4, 7.7, -7.0, 2.1, 0.3, 3.7,
100 3.3, -8.6, 9.8, 3.6, 3.1, 6.5, -2.4, 7.8, 7.5, 8.4, -2.8, -6.3, -5.1, -2.7, 9.3, -0.8,
101 -9.2, 7.9, 8.9, 3.4, 0.1, -5.3, -6.8, 4.9, 4.3, -0.7, -2.2, -3.2, -7.5, -2.3, 0.0, 8.1,
102 -9.2, -2.3, -5.7, 2.1, 2.6, 2.0, 0.3, -8.0, -2.0, -7.9, 6.6, 8.4, 4.0, -6.2, -6.9, -7.2,
103 7.7, -5.0, 5.3, 1.9, -5.3, -7.5, 8.8, 8.3, 9.0, 8.1, 3.2, 1.2, -5.4, -0.2, 2.1, -5.2,
108 class BaseFFTTest : public ::testing::Test
111 BaseFFTTest() : checker_(data_.rootChecker()), flags_(GMX_FFT_FLAG_CONSERVATIVE)
113 // TODO: These tolerances are just something that has been observed
114 // to be sufficient to pass the tests. It would be nicer to
115 // actually argue about why they are sufficient (or what is).
116 checker_.setDefaultTolerance(gmx::test::relativeToleranceAsPrecisionDependentUlp(10.0, 64, 512));
118 ~BaseFFTTest() override { gmx_fft_cleanup(); }
120 gmx::test::TestReferenceData data_;
121 gmx::test::TestReferenceChecker checker_;
122 std::vector<real> in_, out_;
126 class FFTTest : public BaseFFTTest
129 FFTTest() : fft_(nullptr) {}
134 gmx_fft_destroy(fft_);
140 class ManyFFTTest : public BaseFFTTest
143 ManyFFTTest() : fft_(nullptr) {}
144 ~ManyFFTTest() override
148 gmx_many_fft_destroy(fft_);
155 // TODO: Add tests for aligned/not-aligned input/output memory
157 class FFTTest1D : public FFTTest, public ::testing::WithParamInterface<int>
161 class FFFTest3D : public BaseFFTTest
164 FFFTest3D() : fft_(nullptr) {}
165 ~FFFTest3D() override
169 gmx_parallel_3dfft_destroy(fft_);
172 gmx_parallel_3dfft_t fft_;
176 TEST_P(FFTTest1D, Complex)
178 const int nx = GetParam();
179 ASSERT_LE(nx * 2, static_cast<int>(sizeof(inputdata) / sizeof(inputdata[0])));
181 in_ = std::vector<real>(nx * 2);
182 std::copy(inputdata, inputdata + nx * 2, in_.begin());
183 out_ = std::vector<real>(nx * 2);
185 real* out = &out_[0];
187 gmx_fft_init_1d(&fft_, nx, flags_);
189 gmx_fft_1d(fft_, GMX_FFT_FORWARD, in, out);
190 checker_.checkSequenceArray(nx * 2, out, "forward");
191 gmx_fft_1d(fft_, GMX_FFT_BACKWARD, in, out);
192 checker_.checkSequenceArray(nx * 2, out, "backward");
195 TEST_P(FFTTest1D, Real)
197 const int rx = GetParam();
198 const int cx = (rx / 2 + 1);
199 ASSERT_LE(cx * 2, static_cast<int>(sizeof(inputdata) / sizeof(inputdata[0])));
201 in_ = std::vector<real>(cx * 2);
202 std::copy(inputdata, inputdata + cx * 2, in_.begin());
203 out_ = std::vector<real>(cx * 2);
205 real* out = &out_[0];
207 gmx_fft_init_1d_real(&fft_, rx, flags_);
209 gmx_fft_1d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
210 checker_.checkSequenceArray(cx * 2, out, "forward");
211 gmx_fft_1d_real(fft_, GMX_FFT_COMPLEX_TO_REAL, in, out);
212 checker_.checkSequenceArray(rx, out, "backward");
215 INSTANTIATE_TEST_CASE_P(7_8_25_36_60, FFTTest1D, ::testing::Values(7, 8, 25, 36, 60));
218 TEST_F(ManyFFTTest, Complex1DLength48Multi5Test)
223 in_ = std::vector<real>(nx * 2 * N);
224 std::copy(inputdata, inputdata + nx * 2 * N, in_.begin());
225 out_ = std::vector<real>(nx * 2 * N);
227 real* out = &out_[0];
229 gmx_fft_init_many_1d(&fft_, nx, N, flags_);
231 gmx_fft_many_1d(fft_, GMX_FFT_FORWARD, in, out);
232 checker_.checkSequenceArray(nx * 2 * N, out, "forward");
233 gmx_fft_many_1d(fft_, GMX_FFT_BACKWARD, in, out);
234 checker_.checkSequenceArray(nx * 2 * N, out, "backward");
237 TEST_F(ManyFFTTest, Real1DLength48Multi5Test)
240 const int cx = (rx / 2 + 1);
243 in_ = std::vector<real>(cx * 2 * N);
244 std::copy(inputdata, inputdata + cx * 2 * N, in_.begin());
245 out_ = std::vector<real>(cx * 2 * N);
247 real* out = &out_[0];
249 gmx_fft_init_many_1d_real(&fft_, rx, N, flags_);
251 gmx_fft_many_1d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
252 checker_.checkSequenceArray(cx * 2 * N, out, "forward");
253 gmx_fft_many_1d_real(fft_, GMX_FFT_COMPLEX_TO_REAL, in, out);
254 checker_.checkSequenceArray(rx * N, out, "backward");
257 TEST_F(FFTTest, Real2DLength18_15Test)
260 const int cx = (rx / 2 + 1);
263 in_ = std::vector<real>(cx * 2 * ny);
264 std::copy(inputdata, inputdata + cx * 2 * ny, in_.begin());
265 out_ = std::vector<real>(cx * 2 * ny);
267 real* out = &out_[0];
269 gmx_fft_init_2d_real(&fft_, rx, ny, flags_);
271 gmx_fft_2d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
272 checker_.checkSequenceArray(cx * 2 * ny, out, "forward");
273 // known to be wrong for gmx_fft_mkl. And not used.
274 // gmx_fft_2d_real(_fft,GMX_FFT_COMPLEX_TO_REAL,in,out);
275 // _checker.checkSequenceArray(rx*ny, out, "backward");
278 // TODO: test with threads and more than 1 MPI ranks
279 TEST_F(FFFTest3D, Real5_6_9)
281 int ndata[] = { 5, 6, 9 };
282 MPI_Comm comm[] = { MPI_COMM_NULL, MPI_COMM_NULL };
285 ivec local_ndata, offset, rsize, csize, complex_order;
287 gmx_parallel_3dfft_init(&fft_, ndata, &rdata, &cdata, comm, TRUE, 1);
289 gmx_parallel_3dfft_real_limits(fft_, local_ndata, offset, rsize);
290 gmx_parallel_3dfft_complex_limits(fft_, complex_order, local_ndata, offset, csize);
291 checker_.checkVector(rsize, "rsize");
292 checker_.checkVector(csize, "csize");
293 int size = csize[0] * csize[1] * csize[2];
294 int sizeInBytes = size * sizeof(t_complex);
295 int sizeInReals = sizeInBytes / sizeof(real);
297 in_ = std::vector<real>(sizeInReals);
298 // Use std::copy to convert from double to real easily
299 std::copy(inputdata, inputdata + sizeInReals, in_.begin());
300 // Use memcpy to convert to t_complex easily
301 memcpy(rdata, in_.data(), sizeInBytes);
302 gmx_parallel_3dfft_execute(fft_, GMX_FFT_REAL_TO_COMPLEX, 0, nullptr);
303 // TODO use std::complex and add checkComplex for it
304 checker_.checkSequenceArray(size * 2, reinterpret_cast<real*>(cdata), "forward");
306 // Use std::copy to convert from double to real easily
307 std::copy(inputdata, inputdata + sizeInReals, in_.begin());
308 // Use memcpy to convert to t_complex easily
309 memcpy(cdata, in_.data(), sizeInBytes);
310 gmx_parallel_3dfft_execute(fft_, GMX_FFT_COMPLEX_TO_REAL, 0, nullptr);
311 for (int i = 0; i < ndata[0] * ndata[1]; i++) // check sequence but skip unused data
313 checker_.checkSequenceArray(ndata[2], rdata + i * rsize[2],
314 gmx::formatString("backward %d", i).c_str());