From 17097a93ece0df40fbc09877245ac6298d198f92 Mon Sep 17 00:00:00 2001 From: paboyle Date: Wed, 17 Aug 2016 01:33:55 +0100 Subject: [PATCH] FFTW test ran over 4 mpi processes. --- lib/FFT.h | 198 +++++++++++++++++++++++++++++++++ lib/Grid.h | 4 +- lib/lattice/Lattice_transfer.h | 75 ++++++++++++- lib/simd/Grid_vector_types.h | 6 + lib/simd/Grid_vector_unops.h | 10 +- lib/tensors/Tensor_unary.h | 1 + tests/core/Test_fft.cc | 87 +++++++++++++++ 7 files changed, 372 insertions(+), 9 deletions(-) create mode 100644 lib/FFT.h create mode 100644 tests/core/Test_fft.cc diff --git a/lib/FFT.h b/lib/FFT.h new file mode 100644 index 00000000..262e525b --- /dev/null +++ b/lib/FFT.h @@ -0,0 +1,198 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/Cshift.h + + Copyright (C) 2015 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#ifndef _GRID_FFT_H_ +#define _GRID_FFT_H_ + +#include + +namespace Grid { + + + class FFT { + private: + + GridCartesian *vgrid; + GridCartesian *sgrid; + + int Nd; + std::vector dimensions; + std::vector processors; + std::vector processor_coor; + + public: + + static const int forward=FFTW_FORWARD; + static const int backward=FFTW_BACKWARD; + + FFT ( GridCartesian * grid ) : + vgrid(grid), + Nd(grid->_ndimension), + dimensions(grid->_fdimensions), + processors(grid->_processors), + processor_coor(grid->_processor_coor) + { + std::vector layout(Nd,1); + sgrid = new GridCartesian(dimensions,layout,processors); + }; + + ~FFT ( void) { + delete sgrid; + } + + template + void FFT_dim(Lattice &result,const Lattice &source,int dim, int inverse){ + + conformable(result._grid,vgrid); + conformable(source._grid,vgrid); + + int L = vgrid->_ldimensions[dim]; + int G = vgrid->_fdimensions[dim]; + + std::vector layout(Nd,1); + std::vector pencil_gd(vgrid->_fdimensions); + std::vector pencil_ld(processors); + + pencil_gd[dim] = G*processors[dim]; + pencil_ld[dim] = G*processors[dim]; + + // Pencil global vol LxLxGxLxL per node + GridCartesian pencil_g(pencil_gd,layout,processors); + GridCartesian pencil_l(pencil_ld,layout,processors); + + // Construct pencils + typedef typename vobj::scalar_object sobj; + Lattice ssource(vgrid); ssource =source; + Lattice pgsource(&pencil_g); + Lattice pgresult(&pencil_g); + Lattice plsource(&pencil_l); + Lattice plresult(&pencil_l); + + { + + assert(sizeof(typename sobj::scalar_type)==sizeof(ComplexD)); + assert(sizeof(fftw_complex)==sizeof(ComplexD)); + assert(sizeof(fftw_complex)==sizeof(ComplexD)); + + int Ncomp = sizeof(sobj)/sizeof(fftw_complex); + + std::cout << "Ncomp = "<lSites();idx++) { + + std::vector lcoor(Nd); + sgrid->LocalIndexToLocalCoor(idx,lcoor); + + sobj s; + + peekLocalSite(s,ssource,lcoor); + + lcoor[dim]+=p*L; + + pokeLocalSite(s,pgsource,lcoor); + } + + ssource = Cshift(ssource,dim,L); + } + + std::cout << " pgsource pencil " << pgsource<lSites();idx++) { + + std::vector pcoor(Nd,0); + std::vector lcoor(Nd); + sgrid->LocalIndexToLocalCoor(idx,lcoor); + + if ( lcoor[dim] == 0 ) { // restricts loop to plane at lcoor[dim]==0 + + // Project to local pencil array + for(int l=0;l #include #include #include +#include #include #include #include @@ -78,7 +79,8 @@ Author: paboyle #include #include #include -#include + +#include #include #include diff --git a/lib/lattice/Lattice_transfer.h b/lib/lattice/Lattice_transfer.h index 2fa72014..cc4617de 100644 --- a/lib/lattice/Lattice_transfer.h +++ b/lib/lattice/Lattice_transfer.h @@ -349,7 +349,7 @@ void localConvert(const Lattice &in,Lattice &out) assert(ig->_ldimensions[d] == og->_ldimensions[d]); } -PARALLEL_FOR_LOOP + //PARALLEL_FOR_LOOP for(int idx=0;idxlSites();idx++){ std::vector lcoor(ni); ig->LocalIndexToLocalCoor(idx,lcoor); @@ -446,6 +446,79 @@ void ExtractSlice(Lattice &lowDim, Lattice & higherDim,int slice, in } + +template +void InsertSliceLocal(Lattice &lowDim, Lattice & higherDim,int slice_lo,int slice_hi, int orthog) +{ + typedef typename vobj::scalar_object sobj; + sobj s; + + GridBase *lg = lowDim._grid; + GridBase *hg = higherDim._grid; + int nl = lg->_ndimension; + int nh = hg->_ndimension; + + assert(nl == nh); + assert(orthog=0); + + for(int d=0;d_processors[d] == hg->_processors[d]); + assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + } + + // the above should guarantee that the operations are local + //PARALLEL_FOR_LOOP + for(int idx=0;idxlSites();idx++){ + std::vector lcoor(nl); + std::vector hcoor(nh); + lg->LocalIndexToLocalCoor(idx,lcoor); + if( lcoor[orthog] == slice_lo ) { + hcoor=lcoor; + hcoor[orthog] = slice_hi; + peekLocalSite(s,lowDim,lcoor); + pokeLocalSite(s,higherDim,hcoor); + } + } +} + + +template +void ExtractSliceLocal(Lattice &lowDim, Lattice & higherDim,int slice_lo,int slice_hi, int orthog) +{ + typedef typename vobj::scalar_object sobj; + sobj s; + + GridBase *lg = lowDim._grid; + GridBase *hg = higherDim._grid; + int nl = lg->_ndimension; + int nh = hg->_ndimension; + + assert(nl == nh); + assert(orthog=0); + + for(int d=0;d_processors[d] == hg->_processors[d]); + assert(lg->_ldimensions[d] == hg->_ldimensions[d]); + } + + // the above should guarantee that the operations are local + //PARALLEL_FOR_LOOP + for(int idx=0;idxlSites();idx++){ + std::vector lcoor(nl); + std::vector hcoor(nh); + lg->LocalIndexToLocalCoor(idx,lcoor); + if( lcoor[orthog] == slice_lo ) { + hcoor=lcoor; + hcoor[orthog] = slice_hi; + peekLocalSite(s,higherDim,hcoor); + pokeLocalSite(s,lowDim,lcoor); + } + } +} + + template void Replicate(Lattice &coarse,Lattice & fine) { diff --git a/lib/simd/Grid_vector_types.h b/lib/simd/Grid_vector_types.h index 2f2b70c8..30576fb6 100644 --- a/lib/simd/Grid_vector_types.h +++ b/lib/simd/Grid_vector_types.h @@ -388,6 +388,12 @@ class Grid_simd { }; // end of Grid_simd class definition + +inline void permute(ComplexD &y,ComplexD b, int perm) { y=b; } +inline void permute(ComplexF &y,ComplexF b, int perm) { y=b; } +inline void permute(RealD &y,RealD b, int perm) { y=b; } +inline void permute(RealF &y,RealF b, int perm) { y=b; } + //////////////////////////////////////////////////////////////////// // General rotate //////////////////////////////////////////////////////////////////// diff --git a/lib/simd/Grid_vector_unops.h b/lib/simd/Grid_vector_unops.h index a67f4e7d..2afac190 100644 --- a/lib/simd/Grid_vector_unops.h +++ b/lib/simd/Grid_vector_unops.h @@ -67,15 +67,13 @@ template struct AsinRealFunctor { scalar operator()(const scalar &a) const { return asin(real(a)); } }; - template struct LogRealFunctor { scalar operator()(const scalar &a) const { return log(real(a)); } }; - template -struct ExpRealFunctor { - scalar operator()(const scalar &a) const { return exp(real(a)); } +struct ExpFunctor { + scalar operator()(const scalar &a) const { return exp(a); } }; template struct NotFunctor { @@ -85,7 +83,6 @@ template struct AbsRealFunctor { scalar operator()(const scalar &a) const { return std::abs(real(a)); } }; - template struct PowRealFunctor { double y; @@ -135,7 +132,6 @@ template inline Scalar rsqrt(const Scalar &r) { return (RSqrtRealFunctor(), r); } - template inline Grid_simd cos(const Grid_simd &r) { return SimdApply(CosRealFunctor(), r); @@ -162,7 +158,7 @@ inline Grid_simd abs(const Grid_simd &r) { } template inline Grid_simd exp(const Grid_simd &r) { - return SimdApply(ExpRealFunctor(), r); + return SimdApply(ExpFunctor(), r); } template inline Grid_simd Not(const Grid_simd &r) { diff --git a/lib/tensors/Tensor_unary.h b/lib/tensors/Tensor_unary.h index dd05a4a7..92a968df 100644 --- a/lib/tensors/Tensor_unary.h +++ b/lib/tensors/Tensor_unary.h @@ -36,6 +36,7 @@ template inline auto func(const iScalar &z) -> iScalar\ {\ iScalar ret;\ ret._internal = func( (z._internal));\ + std::cout << "Unary "<<#func<<" " << z._internal <<" -> "<< ret._internal <<" "<< typeid(obj).name() < inline auto func(const iVector &z) -> iVector\ diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc new file mode 100644 index 00000000..deed6f0a --- /dev/null +++ b/tests/core/Test_fft.cc @@ -0,0 +1,87 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_cshift.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#include + +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout( { vComplexD::Nsimd(),1,1,1}); + std::vector mpi_layout = GridDefaultMpi(); + + GridCartesian Fine(latt_size,simd_layout,mpi_layout); + + LatticeComplexD one(&Fine); + LatticeComplexD zz(&Fine); + LatticeComplexD C(&Fine); + LatticeComplexD Ctilde(&Fine); + LatticeComplexD coor(&Fine); + + std::vector p({1.0,2.0,3.0,2.0}); + + one = ComplexD(1.0,0.0); + zz = ComplexD(0.0,0.0); + + ComplexD ci(0.0,1.0); + + C=zero; + for(int mu=0;mu<4;mu++){ + RealD TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + C = C - TwoPiL * p[mu] * coor; + } + + std::cout << GridLogMessage<< " C " << C<