src/gromacs/fft/tests/fft.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Tests utilities for fft calculations.
  34  *
  35  * \author Roland Schulz <roland@utk.edu>
  36  * \ingroup module_fft
  37  */
  38 #include <vector>
  39
  40 #include <gtest/gtest.h>
  41
  42 #include "gromacs/fft/fft.h"
  43 #include "gromacs/fft/parallel_3dfft.h"
  44 #include "gromacs/utility/stringutil.h"
  45 #include "testutils/refdata.h"
  46
  47 namespace
  48 {
  49
  50 //! Input data for FFT tests.
  51 const real inputdata[] = { //print ",\n".join([",".join(["%4s"%(random.randint(-99,99)/10.,) for i in range(25)]) for j in range(20)])
  52     -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,
  53     -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,
  54     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,
  55     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,
  56     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,
  57     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,
  58     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,
  59     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,
  60     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,
  61     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,
  62     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,
  63     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,
  64     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,
  65     -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,
  66     -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,
  67     -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,
  68     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,
  69     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,
  70     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,
  71     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,
  72 };
  73
  74
  75 class BaseFFTTest : public ::testing::Test
  76 {
  77     public:
  78         BaseFFTTest() : checker_(data_.rootChecker())
  79         {
  80         }
  81         ~BaseFFTTest()
  82         {
  83             gmx_fft_cleanup();
  84         }
  85         gmx::test::TestReferenceData    data_;
  86         gmx::test::TestReferenceChecker checker_;
  87         std::vector<real>               in_, out_;
  88 };
  89
  90 class FFTTest : public BaseFFTTest
  91 {
  92     public:
  93         FFTTest() : fft_(NULL)
  94         {
  95         }
  96         ~FFTTest()
  97         {
  98             if (fft_)
  99             {
 100                 gmx_fft_destroy(fft_);
 101             }
 102         }
 103         gmx_fft_t fft_;
 104 };
 105
 106 class ManyFFTTest : public BaseFFTTest
 107 {
 108     public:
 109         ManyFFTTest() : fft_(NULL)
 110         {
 111         }
 112         ~ManyFFTTest()
 113         {
 114             if (fft_)
 115             {
 116                 gmx_many_fft_destroy(fft_);
 117             }
 118         }
 119         gmx_fft_t fft_;
 120 };
 121
 122
 123 //TODO: Add tests for aligned/not-aligned input/output memory
 124
 125 class FFTTest1D : public FFTTest, public ::testing::WithParamInterface<int>
 126 {
 127
 128 };
 129
 130 class FFFTest3D : public BaseFFTTest
 131 {
 132     public:
 133         FFFTest3D() : fft_(NULL)
 134         {
 135         }
 136         ~FFFTest3D()
 137         {
 138             if (fft_)
 139             {
 140                 gmx_parallel_3dfft_destroy(fft_);
 141             }
 142         }
 143         gmx_parallel_3dfft_t fft_;
 144 };
 145
 146
 147 TEST_P(FFTTest1D, Complex)
 148 {
 149     const int nx = GetParam();
 150     ASSERT_LE(nx*2, static_cast<int>(sizeof(inputdata)/sizeof(real)));
 151
 152     in_  = std::vector<real>(inputdata, inputdata+nx*2);
 153     out_ = std::vector<real>(nx*2);
 154     real* in  = &in_[0];
 155     real* out = &out_[0];
 156
 157     gmx_fft_init_1d(&fft_, nx, 0);
 158
 159     gmx_fft_1d(fft_, GMX_FFT_FORWARD, in, out);
 160     checker_.checkSequenceArray(nx*2, out, "forward");
 161     gmx_fft_1d(fft_, GMX_FFT_BACKWARD, in, out);
 162     checker_.checkSequenceArray(nx*2, out, "backward");
 163 }
 164
 165 TEST_P(FFTTest1D, Real)
 166 {
 167     const int rx = GetParam();
 168     const int cx = (rx/2+1);
 169     ASSERT_LE(cx*2, static_cast<int>(sizeof(inputdata)/sizeof(real)));
 170
 171     in_  = std::vector<real>(inputdata, inputdata+cx*2);
 172     out_ = std::vector<real>(cx*2);
 173     real* in  = &in_[0];
 174     real* out = &out_[0];
 175
 176     gmx_fft_init_1d_real(&fft_, rx, 0);
 177
 178     gmx_fft_1d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
 179     checker_.checkSequenceArray(cx*2, out, "forward");
 180     gmx_fft_1d_real(fft_, GMX_FFT_COMPLEX_TO_REAL, in, out);
 181     checker_.checkSequenceArray(rx, out, "backward");
 182 }
 183
 184 INSTANTIATE_TEST_CASE_P(7_8_25_36_60,
 185                         FFTTest1D, ::testing::Values(7, 8, 25, 36, 60));
 186
 187
 188 TEST_F(ManyFFTTest, Complex1DLength48Multi5Test)
 189 {
 190     const int nx = 48;
 191     const int N  = 5;
 192
 193     in_  = std::vector<real>(inputdata, inputdata+nx*2*N);
 194     out_ = std::vector<real>(nx*2*N);
 195     real* in  = &in_[0];
 196     real* out = &out_[0];
 197
 198     gmx_fft_init_many_1d(&fft_, nx, N, 0);
 199
 200     gmx_fft_many_1d(fft_, GMX_FFT_FORWARD, in, out);
 201     checker_.checkSequenceArray(nx*2*N, out, "forward");
 202     gmx_fft_many_1d(fft_, GMX_FFT_BACKWARD, in, out);
 203     checker_.checkSequenceArray(nx*2*N, out, "backward");
 204 }
 205
 206 TEST_F(ManyFFTTest, Real1DLength48Multi5Test)
 207 {
 208     const int rx = 48;
 209     const int cx = (rx/2+1);
 210     const int N  = 5;
 211
 212     in_  = std::vector<real>(inputdata, inputdata+cx*2*N);
 213     out_ = std::vector<real>(cx*2*N);
 214     real* in  = &in_[0];
 215     real* out = &out_[0];
 216
 217     gmx_fft_init_many_1d_real(&fft_, rx, N, 0);
 218
 219     gmx_fft_many_1d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
 220     checker_.checkSequenceArray(cx*2*N, out, "forward");
 221     gmx_fft_many_1d_real(fft_, GMX_FFT_COMPLEX_TO_REAL, in, out);
 222     checker_.checkSequenceArray(rx*N, out, "backward");
 223 }
 224
 225 TEST_F(FFTTest, Real2DLength18_15Test)
 226 {
 227     const int rx = 18;
 228     const int cx = (rx/2+1);
 229     const int ny = 15;
 230
 231     in_  = std::vector<real>(inputdata, inputdata+cx*2*ny);
 232     out_ = std::vector<real>(cx*2*ny);
 233     real* in  = &in_[0];
 234     real* out = &out_[0];
 235
 236     gmx_fft_init_2d_real(&fft_, rx, ny, 0);
 237
 238     gmx_fft_2d_real(fft_, GMX_FFT_REAL_TO_COMPLEX, in, out);
 239     checker_.checkSequenceArray(cx*2*ny, out, "forward");
 240 //    known to be wrong for gmx_fft_mkl. And not used.
 241 //    gmx_fft_2d_real(_fft,GMX_FFT_COMPLEX_TO_REAL,in,out);
 242 //    _checker.checkSequenceArray(rx*ny, out, "backward");
 243 }
 244
 245 //TODO: test with threads and more than 1 MPI ranks
 246 TEST_F(FFFTest3D, Real5_6_9)
 247 {
 248     int        ndata[] = {5, 6, 9};
 249     MPI_Comm   comm[]  = {MPI_COMM_NULL, MPI_COMM_NULL};
 250     real     * rdata;
 251     t_complex* cdata;
 252     ivec       local_ndata, offset, rsize, csize, complex_order;
 253
 254     gmx_parallel_3dfft_init(&fft_, ndata, &rdata, &cdata,
 255                             comm, NULL, NULL, TRUE, 1);
 256
 257     gmx_parallel_3dfft_real_limits(fft_, local_ndata, offset, rsize);
 258     gmx_parallel_3dfft_complex_limits(fft_, complex_order,
 259                                       local_ndata, offset, csize);
 260     checker_.checkVector(rsize, "rsize");
 261     checker_.checkVector(csize, "csize");
 262     int size = csize[0]*csize[1]*csize[2];
 263
 264     memcpy(rdata, inputdata, size*sizeof(t_complex));
 265     gmx_parallel_3dfft_execute(fft_, GMX_FFT_REAL_TO_COMPLEX, rdata, cdata,
 266                                0, NULL);
 267     //TODO use std::complex and add checkComplex for it
 268     checker_.checkSequenceArray(size*2,
 269                                 reinterpret_cast<real*>(cdata), "forward");
 270
 271     memcpy(cdata, inputdata, size*sizeof(t_complex));
 272     gmx_parallel_3dfft_execute(fft_, GMX_FFT_COMPLEX_TO_REAL, rdata, cdata,
 273                                0, NULL);
 274     for (int i = 0; i < ndata[0]*ndata[1]; i++) //check sequence but skip unused data
 275     {
 276         checker_.checkSequenceArray(ndata[2], rdata+i*rsize[2],
 277                                     gmx::formatString("backward %d", i).c_str());
 278     }
 279 }
 280
 281 } // namespace