Merge branch 'mpi_reduce_improvement' into 'master'

Adding the possibility to either use MPI Allreduce or Reduce, depending on...

See merge request !4

alpine3d/MPIControl.cc

@@ -204,31 +204,43 @@ void MPIControl::gather(const int& send_value, std::vector<int>& receive_vector,
 	MPI_Gather(&value, 1, MPI_INT, recv_data, 1, MPI_INT, static_cast<int>(root), MPI_COMM_WORLD);
 }
 
-void MPIControl::allreduce_max(double& value)
+void MPIControl::reduce_max(double& value, const bool all)
 {
 	double buffer;
-	MPI_Allreduce(&value, &buffer, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+	if (all)
+		MPI_Allreduce(&value, &buffer, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+	else
+		MPI_Reduce(&value, &buffer, 1, MPI_DOUBLE, MPI_MAX, master_rank(), MPI_COMM_WORLD);
 	value = buffer;
 }
 
-void MPIControl::allreduce_min(double& value)
+void MPIControl::reduce_min(double& value, const bool all)
 {
 	double buffer;
-	MPI_Allreduce(&value, &buffer, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+	if (all)
+		MPI_Allreduce(&value, &buffer, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+	else
+		MPI_Reduce(&value, &buffer, 1, MPI_DOUBLE, MPI_MIN, master_rank(), MPI_COMM_WORLD);
 	value = buffer;
 }
 
-void MPIControl::allreduce_sum(double& value)
+void MPIControl::reduce_sum(double& value, const bool all)
 {
 	double buffer;
-	MPI_Allreduce(&value, &buffer, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+	if (all)
+		MPI_Allreduce(&value, &buffer, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+	else
+		MPI_Reduce(&value, &buffer, 1, MPI_DOUBLE, MPI_SUM, master_rank(), MPI_COMM_WORLD);
 	value = buffer;
 }
 
-void MPIControl::allreduce_sum(int& value)
+void MPIControl::reduce_sum(int& value, const bool all)
 {
 	int buffer;
-	MPI_Allreduce(&value, &buffer, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+	if (all)
+		MPI_Allreduce(&value, &buffer, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+	else
+		MPI_Reduce(&value, &buffer, 1, MPI_INT, MPI_SUM, master_rank(), MPI_COMM_WORLD);
 	value = buffer;
 }
 
@@ -268,10 +280,10 @@ MPIControl::MPIControl() : rank_(0), size_(1), name_( getHostName() )
 
 MPIControl::~MPIControl() {}
 void MPIControl::barrier() const {}
-void MPIControl::allreduce_max(double&) {}
-void MPIControl::allreduce_min(double&) {}
-void MPIControl::allreduce_sum(double&) {}
-void MPIControl::allreduce_sum(int&) {}
+void MPIControl::reduce_max(double&, bool) {}
+void MPIControl::reduce_min(double&, bool) {}
+void MPIControl::reduce_sum(double&, bool) {}
+void MPIControl::reduce_sum(int&, bool) {}
 void MPIControl::gather(const int& val, std::vector<int>& vec, const size_t&) { vec.resize(1, val); }
 #endif
 

alpine3d/MPIControl.h

@@ -117,16 +117,18 @@ class MPIControl
 
 		//@{
 		/**
-		 * Combines values from all processes and distributes the result back to all processes.
-		 * - allreduce_max distributes the maximum of all processes
-		 * - allreduce_min distributes the minimum of all processes
-		 * - allreduce_sum distributes the sum of all processes
+		 * Combines values from all processes and, if all==true, distributes the result back to all processes.
+		 * - reduce_max distributes the maximum of all processes
+		 * - reduce_min distributes the minimum of all processes
+		 * - reduce_sum distributes the sum of all processes
 		 * @param[in,out] value The value that is used to perform the reduction and to hold the result
+		 * @param[in] all True (default): value will be distributed back to all processes (i.e., MPI_Allreduce).
+		 *                False: value will only be available at the master process (i.e., MPI_reduce).
 		 */
-		void allreduce_max(double& value);
-		void allreduce_min(double& value);
-		void allreduce_sum(double& value);
-		void allreduce_sum(int& value);
+		void reduce_max(double& value, const bool all=true);
+		void reduce_min(double& value, const bool all=true);
+		void reduce_sum(double& value, const bool all=true);
+		void reduce_sum(int& value, const bool all=true);
 		//@}
 
 		/**
@@ -191,12 +193,14 @@ class MPIControl
 		#endif
 
 		#ifdef ENABLE_MPI
-		/** @brief Adds up the values of type T from all processes and distributes the sum back to all processes.
-		 * @param obj The obj that is a summand of the global sum and will hold the result
+		/** @brief Adds up the values of type T from all processes and, if all==true, distributes the sum back to all processes.
+		 * @param obj The obj that is summed and will hold the result
+		 * @param all True (default): obj will be distributed back to all processes (i.e., MPI_Allreduce).
+		 *            False: obj will only be available at the master process (i.e., MPI_reduce).
 		 * @note class T must support the extraction and insertion operators, \>\> and \<\<
 		 *       and furthermore support the operator+
 		 */
-		template <class T> void allreduce_sum(T& obj)
+		template <class T> void reduce_sum(T& obj, const bool all=true)
 		{
 			if (size_ <= 1) return;
 
@@ -213,16 +217,19 @@ class MPIControl
 			MPI_Type_contiguous(static_cast<int>(len), MPI_CHAR, &chars_datatype);
 			MPI_Type_commit(&chars_datatype);
 
-			MPI_Allreduce(in_obj_char, out_obj_char, 1, chars_datatype, op, MPI_COMM_WORLD);
+			if (all)
+				MPI_Allreduce(in_obj_char, out_obj_char, 1, chars_datatype, op, MPI_COMM_WORLD);
+			else
+				MPI_Reduce(in_obj_char, out_obj_char, 1, chars_datatype, op, master_rank(), MPI_COMM_WORLD);
 			MPI_Op_free(&op);
 
-			deserialize(out_obj_char, len, obj);
+			if (all || rank() == master_rank()) deserialize(out_obj_char, len, obj);
 
 			delete[] (char*)out_obj_char;
 			delete[] (char*)in_obj_char;
 		}
 		#else
-		template <class T> void allreduce_sum(T& /*obj*/) {}
+		template <class T> void reduce_sum(T& /*obj*/, const bool /*all*/=true) {}
 		#endif
 
 		#ifdef ENABLE_MPI
@@ -405,7 +412,7 @@ class MPIControl
 			const size_t sum = vec_local.size();
 
 			gather(vec_local.size(), vec_sizes); //global offsets
-			allreduce_sum(sum); //global amount of T*
+			reduce_sum(sum, false); //global amount of T*
 
 			if (rank_ == root) {
 				std::string obj_string;

alpine3d/SnowpackInterface.cc

@@ -428,7 +428,7 @@ void SnowpackInterface::writeOutput(const mio::Date& date)
 	const bool isMaster = mpicontrol.master();
 
 	if (do_grid_output(date)) {
-		//no OpenMP pragma here, otherwise multiple threads might call an MPI allreduce_sum()
+		//no OpenMP pragma here, otherwise multiple threads might call an MPI reduce_sum()
 		for (size_t ii=0; ii<output_grids.size(); ii++) {
 			const size_t SnGrids_idx = SnGrids::getParameterIndex( output_grids[ii] );
 			mio::Grid2DObject grid( getGrid( static_cast<SnGrids::Parameters>(SnGrids_idx)) );
@@ -773,7 +773,7 @@ mio::Grid2DObject SnowpackInterface::getGrid(const SnGrids::Parameters& param) c
 		default: ; //so compilers do not complain about missing conditions
 	}
 
-	mio::Grid2DObject o_grid2D(dem, 0.); //so the allreduce_sum works
+	mio::Grid2DObject o_grid2D(dem, 0.); //so the reduce_sum works
 	size_t errCount = 0;
 	mio::Grid2DObject tmp_grid2D_(dem,mpi_offset, 0, mpi_nx, dimy);
 	mio::Grid2DObject tmp_grid2D(tmp_grid2D_,mio::IOUtils::nodata);
@@ -783,10 +783,7 @@ mio::Grid2DObject SnowpackInterface::getGrid(const SnGrids::Parameters& param) c
 		if (!tmp.empty()) {
 			for (size_t i=0; i<tmp_grid2D.getNx(); ++i){
 				for (size_t j=0; j<tmp_grid2D.getNy(); ++j){
-					if (tmp(i,j)!=mio::IOUtils::nodata)
-					{
-							tmp_grid2D(i,j)=tmp(i,j);
-					}
+					if (tmp(i,j)!=mio::IOUtils::nodata) tmp_grid2D(i,j)=tmp(i,j);
 				}
 			}
 		} else {
@@ -795,7 +792,7 @@ mio::Grid2DObject SnowpackInterface::getGrid(const SnGrids::Parameters& param) c
 	}
 	o_grid2D.grid2D.fill(tmp_grid2D.grid2D, mpi_offset, 0, mpi_nx, dimy);
 
-	MPIControl::instance().allreduce_sum(o_grid2D);
+	MPIControl::instance().reduce_sum(o_grid2D, false);
 	//with some MPI implementations, when transfering large amounts of data, the buffers might get full and lead to a crash
 
 	if (errCount>0) {

alpine3d/ebalance/EnergyBalance.cc

@@ -151,9 +151,9 @@ void EnergyBalance::setMeteo(const mio::Grid2DObject& in_ilwr,
 		diffuse.fill(band_diffuse, startx, 0, nx, dimy);
 		direct_unshaded_horizontal.fill(band_direct_unshaded_horizontal, startx, 0, nx, dimy);
 	}
-	MPIControl::instance().allreduce_sum(direct);
-	MPIControl::instance().allreduce_sum(diffuse);
-	MPIControl::instance().allreduce_sum(direct_unshaded_horizontal);
+	MPIControl::instance().reduce_sum(direct);
+	MPIControl::instance().reduce_sum(diffuse);
+	MPIControl::instance().reduce_sum(direct_unshaded_horizontal);
 	double solarAzimuth, solarElevation;
 	radfields[0].getPositionSun(solarAzimuth, solarElevation);
 

alpine3d/ebalance/TerrainRadiationComplex.cc

@@ -748,7 +748,7 @@ void TerrainRadiationComplex::getRadiation(mio::Array2D<double> &direct, mio::Ar
 		}
 	}
 
-	mpicontrol.allreduce_sum(TList_sky_aniso);
+	mpicontrol.reduce_sum(TList_sky_aniso);
 
 	///////////////////////////////
 	///// LAND MAIN LOOP [MT eq. 2.97]
@@ -815,8 +815,8 @@ void TerrainRadiationComplex::getRadiation(mio::Array2D<double> &direct, mio::Ar
 			}
 		}
 
-		mpicontrol.allreduce_sum(terrain_flux_new);
-		mpicontrol.allreduce_sum(TList_ms_new);
+		mpicontrol.reduce_sum(terrain_flux_new);
+		mpicontrol.reduce_sum(TList_ms_new);
 
 		// Test if convergence is below treshold [MT eq. 2.100]
 		if (number_rounds != 0 && TerrainBiggestDifference(terrain_flux_new, terrain_flux_old) < delta_F_max)
@@ -886,9 +886,9 @@ void TerrainRadiationComplex::getRadiation(mio::Array2D<double> &direct, mio::Ar
 	// If SolarPanel-module is used, send all needed data
 	if (_hasSP)
 	{
-		mpicontrol.allreduce_sum(TList_sky_iso);
-		mpicontrol.allreduce_sum(TList_direct);
-		mpicontrol.allreduce_sum(TList_ms_new);
+		mpicontrol.reduce_sum(TList_sky_iso);
+		mpicontrol.reduce_sum(TList_direct);
+		mpicontrol.reduce_sum(TList_ms_new);
 		SP.setTLists(TList_direct, TList_sky_iso, TList_sky_aniso, TList_ms_new + TList_sky_aniso);
 	}
 
@@ -901,8 +901,8 @@ void TerrainRadiationComplex::getRadiation(mio::Array2D<double> &direct, mio::Ar
 			terrain_temp(ii, jj) = (terrain_flux_new(ii, jj, 1) + terrain_flux_new(ii, jj, 0)) / 2.;
 		}
 	}
-	mpicontrol.allreduce_sum(diffuse_temp);
-	mpicontrol.allreduce_sum(terrain_temp);
+	mpicontrol.reduce_sum(diffuse_temp);
+	mpicontrol.reduce_sum(terrain_temp);
 	diffuse = diffuse_temp;
 	terrain = terrain_temp;
 
@@ -1141,7 +1141,7 @@ double TerrainRadiationComplex::computeSkyViewFactor(size_t ii_dem, size_t jj_de
 void TerrainRadiationComplex::getSkyViewFactor(mio::Array2D<double> &o_sky_vf)
 {
 	o_sky_vf = sky_vf_mean;
-	MPIControl::instance().allreduce_sum(o_sky_vf);
+	MPIControl::instance().reduce_sum(o_sky_vf);
 }
 
 double TerrainRadiationComplex::getSkyViewFactor(size_t ii_dem, size_t jj_dem, size_t which_triangle)

alpine3d/ebalance/TerrainRadiationSimple.cc

@@ -44,8 +44,8 @@ void TerrainRadiationSimple::getRadiation(mio::Array2D<double>& /*direct*/,
                                           double /*solarAzimuth*/, double /*solarElevation*/)
 {
 	MPIControl& mpicontrol = MPIControl::instance();
-	terrain.resize(dimx, dimy, 0.);  //so allreduce_sum works properly when it sums full grids
-	Array2D<double> diff_corr(dimx, dimy, 0.); //so allreduce_sum works properly when it sums full grids
+	terrain.resize(dimx, dimy, 0.);  //so reduce_sum works properly when it sums full grids
+	Array2D<double> diff_corr(dimx, dimy, 0.); //so reduce_sum works properly when it sums full grids
 
 	if (mpicontrol.master()) {
 		std::cout << "[i] Calculating terrain radiation with simple method, using " << mpicontrol.size();
@@ -68,8 +68,8 @@ void TerrainRadiationSimple::getRadiation(mio::Array2D<double>& /*direct*/,
 		}
 	}
 
-	mpicontrol.allreduce_sum(terrain);
-	mpicontrol.allreduce_sum(diff_corr);
+	mpicontrol.reduce_sum(terrain);
+	mpicontrol.reduce_sum(diff_corr);
 	diffuse = diff_corr; //return the corrected diffuse radiation
 }
 
@@ -108,7 +108,7 @@ void TerrainRadiationSimple::setMeteo(const mio::Array2D<double>& albedo, const
 
 void TerrainRadiationSimple::getSkyViewFactor(mio::Array2D<double> &o_sky_vf) {
 	o_sky_vf = sky_vf;
-	MPIControl::instance().allreduce_sum(o_sky_vf);
+	MPIControl::instance().reduce_sum(o_sky_vf);
 }
 
 void TerrainRadiationSimple::initSkyViewFactor(const mio::DEMObject &dem)