diff --git a/fsgrid.hpp b/fsgrid.hpp index 7746966..4390140 100644 --- a/fsgrid.hpp +++ b/fsgrid.hpp @@ -141,8 +141,19 @@ template class FsGrid : public FsGridTools{ * \param MPI_Comm The MPI communicator this grid should use. * \param isPeriodic An array specifying, for each dimension, whether it is to be treated as periodic. */ - FsGrid(std::array globalSize, MPI_Comm parent_comm, std::array isPeriodic, FsGridCouplingInformation& coupling) - : globalSize(globalSize),coupling(coupling) { + FsGrid(std::array globalSize, MPI_Comm parent_comm, std::array isPeriodic, FsGridCouplingInformation &coupling) + : parentCom(parent_comm),globalSize(globalSize), coupling(coupling){ + this->init(parent_comm,isPeriodic); + } + + //Copy constructor + FsGrid(const FsGrid &other) + : parentCom(other.parentCom),periodic(other.periodic), globalSize(other.globalSize), coupling(other.coupling){ + init(parentCom,periodic); + this->data = other.data; + } + + void init(MPI_Comm &parent_comm,std::array isPeriodic){ int status; int size; @@ -474,7 +485,7 @@ template class FsGrid : public FsGridTools{ retval[1] = y - localStart[1]; retval[2] = z - localStart[2]; - if(retval[0] > localSize[0] || retval[1] > localSize[1] || retval[2] > localSize[2] + if(retval[0] >= localSize[0] || retval[1] >= localSize[1] || retval[2] >= localSize[2] || retval[0] < 0 || retval[1] < 0 || retval[2] < 0) { return {-1,-1,-1}; } @@ -965,27 +976,68 @@ template class FsGrid : public FsGridTools{ //! Copy the entire data from another FsGrid of the same signature over. FsGrid& operator=(const FsGrid& other) { - // Don't copy if sizes mismatch. - // (Should this instead crash the program?) - if(other.localSize[0] != localSize[0] || - other.localSize[1] != localSize[1] || - other.localSize[2] != localSize[2]) { - return *this; - } + // Fail if sizes mismatch. + assert(other.localSize[0] == this->localSize[0]); + assert(other.localSize[1] == this->localSize[1]); + assert(other.localSize[2] == this->localSize[2]); data = other.data; + return *this; + } + FsGrid& operator+=(const FsGrid& rhs) { + assert(this->data.size()==rhs.data.size()); + for (size_t i=0; i< this->data.size(); i++){ + this->data[i] = this->data[i]+rhs.data[i]; + } + return *this; + } + + FsGrid& operator-=(const FsGrid& rhs) { + assert(this->data.size()==rhs.data.size()); + for (size_t i=0; i< this->data.size(); i++){ + this->data[i] = this->data[i]-rhs.data[i]; + } + return *this; + } + + FsGrid& operator*=(const FsGrid& rhs) { + assert(this->data.size()==rhs.data.size()); + for (size_t i=0; i< this->data.size(); i++){ + this->data[i] = this->data[i]*rhs.data[i]; + } + return *this; + } + + FsGrid& operator/=(const FsGrid& rhs) { + assert(this->data.size()==rhs.data.size()); + for (size_t i=0; i< this->data.size(); i++){ + this->data[i] = this->data[i]/rhs.data[i]; + } + return *this; + } + + template ::value>::type * = nullptr> + FsGrid &operator*=(S factor){ + for (size_t i=0; i< this->data.size(); i++){ + this->data[i] = this->data[i]*factor; + } + return *this; + } + + template ::value>::type * = nullptr> + FsGrid &operator/=(S factor){ + for (size_t i=0; i< this->data.size(); i++){ + this->data[i] = this->data[i]/factor; + } return *this; } - - - private: //! MPI Cartesian communicator used in this grid - MPI_Comm comm3d; - int rank; //!< This task's rank in the communicator - std::vector requests; - uint numRequests; + MPI_Comm comm3d, parentCom; + int rank; //!< This task's rank in the communicator + std::vector requests; + uint numRequests; std::array neighbour; //!< Tasks of the 26 neighbours (plus ourselves) std::vector neighbour_index; //!< Lookup table from rank to index in the neighbour array @@ -1040,3 +1092,110 @@ template class FsGrid : public FsGridTools{ array[2] = a; } }; + + +//Array Operator Overloads +template +static inline std::array operator+(const std::array& arr1, const std::array& arr2){ + std::array res; + for (size_t i = 0; i < N; ++i){ + res[i] = arr1[i] + arr2[i]; + } + return res; +} + +template +static inline std::array operator-(const std::array& arr1, const std::array& arr2){ + std::array res; + for (size_t i = 0; i < N; ++i){ + res[i] = arr1[i] - arr2[i]; + } + return res; +} + +template +static inline std::array operator*(const std::array &arr1, const std::array &arr2){ + std::array res; + for (size_t i = 0; i < N; ++i){ + res[i] = arr1[i] * arr2[i]; + } + return res; +} + +template +static inline std::array operator/(const std::array& arr1, const std::array& arr2){ + std::array res; + for (size_t i = 0; i < N; ++i){ + res[i] = arr1[i] / arr2[i]; + } + return res; +} + +template ::value>::type* = nullptr> +std::array operator*(std::array& arr, S factor){ + std::array res; + for (size_t i = 0; i < N; ++i){ + res[i]=arr[i]*factor; + } + return res; +} + +template ::value>::type* = nullptr> +std::array operator/(std::array& arr, S factor){ + std::array res; + for (size_t i = 0; i < N; ++i){ + res[i]=arr[i]/factor; + } + return res; +} + + +// FsGrid Operator Overloads +template +static inline FsGrid operator+(const FsGrid &lhs, const FsGrid &rhs){ + FsGrid a(lhs); + return a+=rhs; +} + +template +static inline FsGrid operator-(const FsGrid &lhs, const FsGrid &rhs){ + FsGrid a(lhs); + return a-=rhs; +} + +template +static inline FsGrid operator*(const FsGrid &lhs, const FsGrid &rhs){ + FsGrid a(lhs); + return a *= rhs; +} + +template ::value>::type * = nullptr> +static inline FsGrid operator*(const FsGrid& lhs, S factor){ + FsGrid a(lhs); + return a*=factor; +} + +template ::value>::type * = nullptr> +static inline FsGrid operator/(const FsGrid& lhs, S factor){ + FsGrid a(lhs); + return a/=factor; +} + +template +static inline FsGrid operator/(const FsGrid &lhs, const FsGrid &rhs){ + FsGrid a(lhs); + return a /= rhs; +} + + +/* Here tNorm is the interpolation time but normalized. + * So t0=0 , t1=1 and then for any t>1 we have a + * linear extrapolation */ +template ::value>::type* = nullptr> +static inline FsGrid lerp_t(const FsGrid &lhs, const FsGrid &rhs,S tNorm){ + //If sizes mismatch then this will make it crash in the '=' operator later on + //Do interpolate-targetGrid=lhs(1-t)+rhs*t1 + return lhs*(1. - tNorm) + rhs*tNorm; +} + + diff --git a/tests/Makefile b/tests/Makefile index 1c2427e..0618e02 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -1,14 +1,16 @@ -CXX=CC +CXX=mpic++ CXXFLAGS= -O3 -std=c++11 -march=native -g -Wall -all: benchmark test ddtest +all: benchmark test ddtest operators benchmark: benchmark.cpp ../fsgrid.hpp $(CXX) $(CXXFLAGS) -o $@ $< +operators: operators.cpp ../fsgrid.hpp + $(CXX) $(CXXFLAGS) -o $@ $< test: test.cpp ../fsgrid.hpp $(CXX) $(CXXFLAGS) -o $@ $< ddtest: ddtest.cpp $(CXX) $(CXXFLAGS) -o $@ $< clean: - -rm test + -rm test benchmark ddtest operators diff --git a/tests/benchmark.cpp b/tests/benchmark.cpp index 25808ea..14e25fd 100644 --- a/tests/benchmark.cpp +++ b/tests/benchmark.cpp @@ -22,8 +22,9 @@ */ template void timeit(std::array globalSize, std::array isPeriodic, int iterations){ - double t1,t2; - FsGrid testGrid(globalSize, MPI_COMM_WORLD, isPeriodic); + double t1,t2; + FsGridCouplingInformation gridCoupling; + FsGrid testGrid(globalSize, MPI_COMM_WORLD, isPeriodic,gridCoupling); int rank,size; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &size); diff --git a/tests/operators.cpp b/tests/operators.cpp new file mode 100644 index 0000000..4b71804 --- /dev/null +++ b/tests/operators.cpp @@ -0,0 +1,81 @@ +#include +#include "../fsgrid.hpp" +/* + Copyright (C) 2016 Finnish Meteorological Institute + Copyright (C) 2016 CSC -IT Center for Science + + This file is part of fsgrid + + fsgrid 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 3 of the License, or + (at your option) any later version. + + fsgrid 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 fsgrid. If not, see . +*/ + + + +int main(int argc, char** argv) { + typedef float Real; + MPI_Init(&argc,&argv); + + int rank,size; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + + // Create a 8×8 Testgrid + std::array globalSize{2000,1000,1}; + std::array globalSize2{3000,1000,1}; + std::array isPeriodic{false,false,true}; + + + + FsGridCouplingInformation gridCoupling; + FsGrid, 2> A(globalSize, MPI_COMM_WORLD, isPeriodic, gridCoupling); + FsGrid, 2> B(globalSize, MPI_COMM_WORLD, isPeriodic, gridCoupling); + FsGrid, 2> C(globalSize, MPI_COMM_WORLD, isPeriodic, gridCoupling); + FsGrid, 2> D(globalSize2, MPI_COMM_WORLD, isPeriodic, gridCoupling); + + // FsGrid, 2> newGuy(A + B); + // + + A=B+C; + // += + A+=B; + // - + A=A-C; + // -= + A-=B; + D=D*3.0; + C=lerp_t(A,B,5.,10.,7.); + //supposed to fail; + // A=B+D; + + + + + + + + + + + + + + + + + + + + MPI_Finalize(); + return 0; +} diff --git a/tests/test.cpp b/tests/test.cpp index 86ed4f1..df75760 100644 --- a/tests/test.cpp +++ b/tests/test.cpp @@ -33,8 +33,9 @@ int main(int argc, char** argv) { // Create a 8×8 Testgrid std::array globalSize{20,20,1}; std::array isPeriodic{false,false,true}; + FsGridCouplingInformation gridCoupling; { - FsGrid testGrid(globalSize, MPI_COMM_WORLD, isPeriodic); + FsGrid testGrid(globalSize, MPI_COMM_WORLD, isPeriodic,gridCoupling); /* if(rank == 0) { std::cerr << " --- Test task mapping functions ---" << std::endl; @@ -190,9 +191,7 @@ int main(int argc, char** argv) { } } } - - MPI_Finalize(); return 0;