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,2018, 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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Tests utilities for fft calculations.
  38  *
  39  * Current reference data is generated in double precision using the Reference
  40  * build type, except for the compiler (Apple Clang).
  41  *
  42  * \author Roland Schulz <roland@utk.edu>
  43  * \ingroup module_fft
  44  */
  45 #include "gmxpre.h"
  46
  47 #include "gromacs/fft/fft.h"
  48
  49 #include <algorithm>
  50 #include <vector>
  51
  52 #include <gtest/gtest.h>
  53
  54 #include "gromacs/fft/parallel_3dfft.h"
  55 #include "gromacs/utility/stringutil.h"
  56
  57 #include "testutils/refdata.h"
  58 #include "testutils/testasserts.h"
  59
  60 namespace
  61 {
  62
  63 /*! \brief Input data for FFT tests.
  64  *
  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.
  70  */
  71 const double inputdata[] = { //print ",\n".join([",".join(["%4s"%(random.randint(-99,99)/10.,) for i in range(25)]) for j in range(20)])
  72     -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, -1.9, -7.6, 1.4, -3.5, 0.7, 5.6, -4.2, -1.1, -4.4,
  73     -6.3, -7.2, 4.6, -3.0, -0.9, 7.2, 2.5, -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, -3.9, -3.6,
  74     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, -7.7, 1.9, -8.4, 4.4, 2.3, -2.9, 6.7, 2.7, 5.8, -3.6, 8.9,
  75     8.9, 4.3, 9.1, 9.3, -8.7, 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, 9.8, 1.3, 3.2, -0.2,
  76     9.9, 4.4, -9.9, -7.2, 4.4, 4.7, 7.2, -0.3, 0.3, -2.1, 8.4, -2.1, -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,
  77     8.3, 5.2, -9.0, 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, -9.5, -0.8, 1.9, -6.2, 4.3, -3.8,
  78     8.6, -1.9, -2.1, -0.4, -7.1, -3.7, 9.1, -6.4, -0.6, 2.5, 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,
  79     2.1, 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, -5.5, 2.6, -3.5, -0.7, 9.0, 9.8, -8.0, 3.6,
  80     3.0, -2.2, -2.8, 0.8, 9.0, 2.8, 7.7, -0.7, -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, -2.5,
  81     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, 4.7, 3.9, -5.4, 9.1, -1.6, -1.9, -4.2, -2.6, 0.6, -5.1,
  82     1.8, 5.2, 4.0, -6.2, 6.5, -9.1, 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, 9.3, -7.9, 7.4,
  83     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, -4.9, 3.5, 7.7, -4.2, -9.7, 2.4, 8.1, 5.9, 3.4, -7.5, 7.5, 2.6,
  84     4.7, 2.7, 2.2, 2.6, 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, -1.7, -5.3, 0.3, 3.9, -3.6,
  85     -3.6, 4.7, -8.1, 1.4, 4.0, 1.3, -4.3, -8.8, -7.3, 6.3, -7.5, -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,
  86     -0.8, -0.9, -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, -6.6, 9.8, -1.1, 8.9, 5.0, 2.9, 0.2,
  87     -2.9, 0.8, 6.7, -0.6, 0.6, 4.1, 5.3, -1.7, -0.3, 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,
  88     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, 8.8, -3.3, -9.2, 4.6, 1.8, 4.6, 2.9, -2.7, 4.2,
  89     7.3, -0.4, 7.7, -7.0, 2.1, 0.3, 3.7, 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, -9.2, 7.9,
  90     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, -9.2, -2.3, -5.7, 2.1, 2.6, 2.0, 0.3, -8.0, -2.0, -7.9, 6.6,
  91     8.4, 4.0, -6.2, -6.9, -7.2, 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, 9.5, 5.9, 5.6, -7.8,
  92 };
  93
  94
  95 class BaseFFTTest : public ::testing::Test
  96 {
  97     public:
  98         BaseFFTTest()
  99             : checker_(data_.rootChecker()), flags_(GMX_FFT_FLAG_CONSERVATIVE)
 100         {
 101             // TODO: These tolerances are just something that has been observed
 102             // to be sufficient to pass the tests.  It would be nicer to
 103             // actually argue about why they are sufficient (or what is).
 104             checker_.setDefaultTolerance(
 105                     gmx::test::relativeToleranceAsPrecisionDependentUlp(10.0, 64, 512));
 106         }
 107         ~BaseFFTTest() override
 108         {
 109             gmx_fft_cleanup();
 110         }
 111
 112         gmx::test::TestReferenceData    data_;
 113         gmx::test::TestReferenceChecker checker_;
 114         std::vector<real>               in_, out_;
 115         int                             flags_;
 116 };
 117
 118 class FFTTest : public BaseFFTTest
 119 {
 120     public:
 121         FFTTest() : fft_(nullptr)
 122         {
 123         }
 124         ~FFTTest() override
 125         {
 126             if (fft_)
 127             {
 128                 gmx_fft_destroy(fft_);
 129             }
 130         }
 131         gmx_fft_t fft_;
 132 };
 133
 134 class ManyFFTTest : public BaseFFTTest
 135 {
 136     public:
 137         ManyFFTTest() : fft_(nullptr)
 138         {
 139         }
 140         ~ManyFFTTest() override
 141         {
 142             if (fft_)
 143             {
 144                 gmx_many_fft_destroy(fft_);
 145             }
 146         }
 147         gmx_fft_t fft_;
 148 };
 149
 150
 151 //TODO: Add tests for aligned/not-aligned input/output memory
 152
 153 class FFTTest1D : public FFTTest, public ::testing::WithParamInterface<int>
 154 {
 155
 156 };
 157
 158 class FFFTest3D : public BaseFFTTest
 159 {
 160     public:
 161         FFFTest3D() : fft_(nullptr)
 162         {
 163         }
 164         ~FFFTest3D() override
 165         {
 166             if (fft_)
 167             {
 168                 gmx_parallel_3dfft_destroy(fft_);
 169             }
 170         }
 171         gmx_parallel_3dfft_t fft_;
 172 };
 173
 174
 175 TEST_P(FFTTest1D, Complex)
 176 {
 177     const int nx = GetParam();
 178     ASSERT_LE(nx*2, static_cast<int>(sizeof(inputdata)/sizeof(inputdata[0])));
 179
 180     in_  = std::vector<real>(nx*2);
 181     std::copy(inputdata, inputdata+nx*2, in_.begin());
 182     out_ = std::vector<real>(nx*2);
 183     real* in  = &in_[0];
 184     real* out = &out_[0];
 185
 186     gmx_fft_init_1d(&fft_, nx, flags_);
 187
 188     gmx_fft_1d(fft_, GMX_FFT_FORWARD, in, out);
 189     checker_.checkSequenceArray(nx*2, out, "forward");
 190     gmx_fft_1d(fft_, GMX_FFT_BACKWARD, in, out);
 191     checker_.checkSequenceArray(nx*2, out, "backward");
 192 }
 193
 194 TEST_P(FFTTest1D, Real)
 195 {
 196     const int rx = GetParam();
 197     const int cx = (rx/2+1);
 198     ASSERT_LE(cx*2, static_cast<int>(sizeof(inputdata)/sizeof(inputdata[0])));
 199
 200     in_  = std::vector<real>(cx*2);
 201     std::copy(inputdata, inputdata+cx*2, in_.begin());
 202     out_ = std::vector<real>(cx*2);
 203     real* in  = &in_[0];
 204     real* out = &out_[0];
 205
 206     gmx_fft_init_1d_real(&fft_, rx, flags_);
 207
 208     gmx_fft_1d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
 209     checker_.checkSequenceArray(cx*2, out, "forward");
 210     gmx_fft_1d_real(fft_, GMX_FFT_COMPLEX_TO_REAL, in, out);
 211     checker_.checkSequenceArray(rx, out, "backward");
 212 }
 213
 214 INSTANTIATE_TEST_CASE_P(7_8_25_36_60,
 215                         FFTTest1D, ::testing::Values(7, 8, 25, 36, 60));
 216
 217
 218 TEST_F(ManyFFTTest, Complex1DLength48Multi5Test)
 219 {
 220     const int nx = 48;
 221     const int N  = 5;
 222
 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);
 226     real* in  = &in_[0];
 227     real* out = &out_[0];
 228
 229     gmx_fft_init_many_1d(&fft_, nx, N, flags_);
 230
 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");
 235 }
 236
 237 TEST_F(ManyFFTTest, Real1DLength48Multi5Test)
 238 {
 239     const int rx = 48;
 240     const int cx = (rx/2+1);
 241     const int N  = 5;
 242
 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);
 246     real* in  = &in_[0];
 247     real* out = &out_[0];
 248
 249     gmx_fft_init_many_1d_real(&fft_, rx, N, flags_);
 250
 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");
 255 }
 256
 257 TEST_F(FFTTest, Real2DLength18_15Test)
 258 {
 259     const int rx = 18;
 260     const int cx = (rx/2+1);
 261     const int ny = 15;
 262
 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);
 266     real* in  = &in_[0];
 267     real* out = &out_[0];
 268
 269     gmx_fft_init_2d_real(&fft_, rx, ny, flags_);
 270
 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");
 276 }
 277
 278 //TODO: test with threads and more than 1 MPI ranks
 279 TEST_F(FFFTest3D, Real5_6_9)
 280 {
 281     int        ndata[] = {5, 6, 9};
 282     MPI_Comm   comm[]  = {MPI_COMM_NULL, MPI_COMM_NULL};
 283     real     * rdata;
 284     t_complex* cdata;
 285     ivec       local_ndata, offset, rsize, csize, complex_order;
 286
 287     gmx_parallel_3dfft_init(&fft_, ndata, &rdata, &cdata,
 288                             comm, TRUE, 1);
 289
 290     gmx_parallel_3dfft_real_limits(fft_, local_ndata, offset, rsize);
 291     gmx_parallel_3dfft_complex_limits(fft_, complex_order,
 292                                       local_ndata, offset, csize);
 293     checker_.checkVector(rsize, "rsize");
 294     checker_.checkVector(csize, "csize");
 295     int size        = csize[0]*csize[1]*csize[2];
 296     int sizeInBytes = size*sizeof(t_complex);
 297     int sizeInReals = sizeInBytes/sizeof(real);
 298
 299     in_  = std::vector<real>(sizeInReals);
 300     // Use std::copy to convert from double to real easily
 301     std::copy(inputdata, inputdata+sizeInReals, in_.begin());
 302     // Use memcpy to convert to t_complex easily
 303     memcpy(rdata, in_.data(), sizeInBytes);
 304     gmx_parallel_3dfft_execute(fft_, GMX_FFT_REAL_TO_COMPLEX, 0, nullptr);
 305     //TODO use std::complex and add checkComplex for it
 306     checker_.checkSequenceArray(size*2,
 307                                 reinterpret_cast<real*>(cdata), "forward");
 308
 309     // Use std::copy to convert from double to real easily
 310     std::copy(inputdata, inputdata+sizeInReals, in_.begin());
 311     // Use memcpy to convert to t_complex easily
 312     memcpy(cdata, in_.data(), sizeInBytes);
 313     gmx_parallel_3dfft_execute(fft_, GMX_FFT_COMPLEX_TO_REAL, 0, nullptr);
 314     for (int i = 0; i < ndata[0]*ndata[1]; i++) //check sequence but skip unused data
 315     {
 316         checker_.checkSequenceArray(ndata[2], rdata+i*rsize[2],
 317                                     gmx::formatString("backward %d", i).c_str());
 318     }
 319 }
 320
 321 } // namespace