src/gromacs/fft/tests/fft.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2016,2017 by the GROMACS development team.
   5  * Copyright (c) 2018,2019,2020, 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 /*! \internal \file
  37  * \brief
  38  * Tests utilities for fft calculations.
  39  *
  40  * Current reference data is generated in double precision using the Reference
  41  * build type, except for the compiler (Apple Clang).
  42  *
  43  * \author Roland Schulz <roland@utk.edu>
  44  * \ingroup module_fft
  45  */
  46 #include "gmxpre.h"
  47
  48 #include "gromacs/fft/fft.h"
  49
  50 #include <algorithm>
  51 #include <vector>
  52
  53 #include <gtest/gtest.h>
  54
  55 #include "gromacs/fft/parallel_3dfft.h"
  56 #include "gromacs/utility/stringutil.h"
  57
  58 #include "testutils/refdata.h"
  59 #include "testutils/testasserts.h"
  60
  61 namespace
  62 {
  63
  64 /*! \brief Input data for FFT tests.
  65  *
  66  * TODO If we require compilers that all support C++11 user literals,
  67  * then this array could be of type real, initialized with e.g. -3.5_r
  68  * that does not suffer from implicit narrowing with brace
  69  * initializers, and we would not have to do so much useless copying
  70  * during the unit tests below.
  71  */
  72 const double inputdata[] = {
  73     // print ",\n".join([",".join(["%4s"%(random.randint(-99,99)/10.,) for i in range(25)]) for j in range(20)])
  74     -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,
  75     -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,
  76     -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,
  77     -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,
  78     -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,
  79     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,
  80     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,
  81     -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,
  82     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,
  83     -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,
  84     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,
  85     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,
  86     -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,
  87     -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,
  88     -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,
  89     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,
  90     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,
  91     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,
  92     -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,
  93     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,
  94     -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,
  95     -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,
  96     -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,
  97     -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,
  98     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,
  99     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,
 100     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,
 101     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,
 102     -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,
 103     -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,
 104     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,
 105     9.5,  5.9,  5.6,  -7.8,
 106 };
 107
 108
 109 class BaseFFTTest : public ::testing::Test
 110 {
 111 public:
 112     BaseFFTTest() : checker_(data_.rootChecker()), flags_(GMX_FFT_FLAG_CONSERVATIVE)
 113     {
 114         // TODO: These tolerances are just something that has been observed
 115         // to be sufficient to pass the tests.  It would be nicer to
 116         // actually argue about why they are sufficient (or what is).
 117         checker_.setDefaultTolerance(gmx::test::relativeToleranceAsPrecisionDependentUlp(10.0, 64, 512));
 118     }
 119     ~BaseFFTTest() override { gmx_fft_cleanup(); }
 120
 121     gmx::test::TestReferenceData    data_;
 122     gmx::test::TestReferenceChecker checker_;
 123     std::vector<real>               in_, out_;
 124     int                             flags_;
 125 };
 126
 127 class FFTTest : public BaseFFTTest
 128 {
 129 public:
 130     FFTTest() : fft_(nullptr) {}
 131     ~FFTTest() override
 132     {
 133         if (fft_)
 134         {
 135             gmx_fft_destroy(fft_);
 136         }
 137     }
 138     gmx_fft_t fft_;
 139 };
 140
 141 class ManyFFTTest : public BaseFFTTest
 142 {
 143 public:
 144     ManyFFTTest() : fft_(nullptr) {}
 145     ~ManyFFTTest() override
 146     {
 147         if (fft_)
 148         {
 149             gmx_many_fft_destroy(fft_);
 150         }
 151     }
 152     gmx_fft_t fft_;
 153 };
 154
 155
 156 // TODO: Add tests for aligned/not-aligned input/output memory
 157
 158 class FFTTest1D : public FFTTest, public ::testing::WithParamInterface<int>
 159 {
 160 };
 161
 162 class FFFTest3D : public BaseFFTTest
 163 {
 164 public:
 165     FFFTest3D() : fft_(nullptr) {}
 166     ~FFFTest3D() override
 167     {
 168         if (fft_)
 169         {
 170             gmx_parallel_3dfft_destroy(fft_);
 171         }
 172     }
 173     gmx_parallel_3dfft_t fft_;
 174 };
 175
 176
 177 TEST_P(FFTTest1D, Complex)
 178 {
 179     const int nx = GetParam();
 180     ASSERT_LE(nx * 2, static_cast<int>(sizeof(inputdata) / sizeof(inputdata[0])));
 181
 182     in_ = std::vector<real>(nx * 2);
 183     std::copy(inputdata, inputdata + nx * 2, in_.begin());
 184     out_      = std::vector<real>(nx * 2);
 185     real* in  = &in_[0];
 186     real* out = &out_[0];
 187
 188     gmx_fft_init_1d(&fft_, nx, flags_);
 189
 190     gmx_fft_1d(fft_, GMX_FFT_FORWARD, in, out);
 191     checker_.checkSequenceArray(nx * 2, out, "forward");
 192     gmx_fft_1d(fft_, GMX_FFT_BACKWARD, in, out);
 193     checker_.checkSequenceArray(nx * 2, out, "backward");
 194 }
 195
 196 TEST_P(FFTTest1D, Real)
 197 {
 198     const int rx = GetParam();
 199     const int cx = (rx / 2 + 1);
 200     ASSERT_LE(cx * 2, static_cast<int>(sizeof(inputdata) / sizeof(inputdata[0])));
 201
 202     in_ = std::vector<real>(cx * 2);
 203     std::copy(inputdata, inputdata + cx * 2, in_.begin());
 204     out_      = std::vector<real>(cx * 2);
 205     real* in  = &in_[0];
 206     real* out = &out_[0];
 207
 208     gmx_fft_init_1d_real(&fft_, rx, flags_);
 209
 210     gmx_fft_1d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
 211     checker_.checkSequenceArray(cx * 2, out, "forward");
 212     gmx_fft_1d_real(fft_, GMX_FFT_COMPLEX_TO_REAL, in, out);
 213     checker_.checkSequenceArray(rx, out, "backward");
 214 }
 215
 216 INSTANTIATE_TEST_CASE_P(7_8_25_36_60, FFTTest1D, ::testing::Values(7, 8, 25, 36, 60));
 217
 218
 219 TEST_F(ManyFFTTest, Complex1DLength48Multi5Test)
 220 {
 221     const int nx = 48;
 222     const int N  = 5;
 223
 224     in_ = std::vector<real>(nx * 2 * N);
 225     std::copy(inputdata, inputdata + nx * 2 * N, in_.begin());
 226     out_      = std::vector<real>(nx * 2 * N);
 227     real* in  = &in_[0];
 228     real* out = &out_[0];
 229
 230     gmx_fft_init_many_1d(&fft_, nx, N, flags_);
 231
 232     gmx_fft_many_1d(fft_, GMX_FFT_FORWARD, in, out);
 233     checker_.checkSequenceArray(nx * 2 * N, out, "forward");
 234     gmx_fft_many_1d(fft_, GMX_FFT_BACKWARD, in, out);
 235     checker_.checkSequenceArray(nx * 2 * N, out, "backward");
 236 }
 237
 238 TEST_F(ManyFFTTest, Real1DLength48Multi5Test)
 239 {
 240     const int rx = 48;
 241     const int cx = (rx / 2 + 1);
 242     const int N  = 5;
 243
 244     in_ = std::vector<real>(cx * 2 * N);
 245     std::copy(inputdata, inputdata + cx * 2 * N, in_.begin());
 246     out_      = std::vector<real>(cx * 2 * N);
 247     real* in  = &in_[0];
 248     real* out = &out_[0];
 249
 250     gmx_fft_init_many_1d_real(&fft_, rx, N, flags_);
 251
 252     gmx_fft_many_1d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
 253     checker_.checkSequenceArray(cx * 2 * N, out, "forward");
 254     gmx_fft_many_1d_real(fft_, GMX_FFT_COMPLEX_TO_REAL, in, out);
 255     checker_.checkSequenceArray(rx * N, out, "backward");
 256 }
 257
 258 TEST_F(FFTTest, Real2DLength18_15Test)
 259 {
 260     const int rx = 18;
 261     const int cx = (rx / 2 + 1);
 262     const int ny = 15;
 263
 264     in_ = std::vector<real>(cx * 2 * ny);
 265     std::copy(inputdata, inputdata + cx * 2 * ny, in_.begin());
 266     out_      = std::vector<real>(cx * 2 * ny);
 267     real* in  = &in_[0];
 268     real* out = &out_[0];
 269
 270     gmx_fft_init_2d_real(&fft_, rx, ny, flags_);
 271
 272     gmx_fft_2d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
 273     checker_.checkSequenceArray(cx * 2 * ny, out, "forward");
 274     //    known to be wrong for gmx_fft_mkl. And not used.
 275     //    gmx_fft_2d_real(_fft,GMX_FFT_COMPLEX_TO_REAL,in,out);
 276     //    _checker.checkSequenceArray(rx*ny, out, "backward");
 277 }
 278
 279 // TODO: test with threads and more than 1 MPI ranks
 280 TEST_F(FFFTest3D, Real5_6_9)
 281 {
 282     int        ndata[] = { 5, 6, 9 };
 283     MPI_Comm   comm[]  = { MPI_COMM_NULL, MPI_COMM_NULL };
 284     real*      rdata;
 285     t_complex* cdata;
 286     ivec       local_ndata, offset, rsize, csize, complex_order;
 287
 288     gmx_parallel_3dfft_init(&fft_, ndata, &rdata, &cdata, comm, TRUE, 1);
 289
 290     gmx_parallel_3dfft_real_limits(fft_, local_ndata, offset, rsize);
 291     gmx_parallel_3dfft_complex_limits(fft_, complex_order, local_ndata, offset, csize);
 292     checker_.checkVector(rsize, "rsize");
 293     checker_.checkVector(csize, "csize");
 294     int size        = csize[0] * csize[1] * csize[2];
 295     int sizeInBytes = size * sizeof(t_complex);
 296     int sizeInReals = sizeInBytes / sizeof(real);
 297
 298     in_ = std::vector<real>(sizeInReals);
 299     // Use std::copy to convert from double to real easily
 300     std::copy(inputdata, inputdata + sizeInReals, in_.begin());
 301     // Use memcpy to convert to t_complex easily
 302     memcpy(rdata, in_.data(), sizeInBytes);
 303     gmx_parallel_3dfft_execute(fft_, GMX_FFT_REAL_TO_COMPLEX, 0, nullptr);
 304     // TODO use std::complex and add checkComplex for it
 305     checker_.checkSequenceArray(size * 2, reinterpret_cast<real*>(cdata), "forward");
 306
 307     // Use std::copy to convert from double to real easily
 308     std::copy(inputdata, inputdata + sizeInReals, in_.begin());
 309     // Use memcpy to convert to t_complex easily
 310     memcpy(cdata, in_.data(), sizeInBytes);
 311     gmx_parallel_3dfft_execute(fft_, GMX_FFT_COMPLEX_TO_REAL, 0, nullptr);
 312     for (int i = 0; i < ndata[0] * ndata[1]; i++) // check sequence but skip unused data
 313     {
 314         checker_.checkSequenceArray(ndata[2], rdata + i * rsize[2],
 315                                     gmx::formatString("backward %d", i).c_str());
 316     }
 317 }
 318
 319 } // namespace