492 hpcg on steroids

include/graphblas/algorithms/hpcg/greedy_coloring.hpp

@@ -0,0 +1,191 @@
+
+/*
+ *   Copyright 2022 Huawei Technologies Co., Ltd.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/**
+ * @file greedy_coloring.hpp
+ * @author Alberto Scolari ([email protected])
+ * Utilities to partition the elements of a mesh via a simple, greedy coloring algorithm.
+ */
+
+#ifndef _H_GRB_ALGORITHMS_HPCG_GREEDY_COLORING
+#define _H_GRB_ALGORITHMS_HPCG_GREEDY_COLORING
+
+#include <cstddef>
+#include <vector>
+
+#include <graphblas/utils/multigrid/linearized_halo_ndim_system.hpp>
+
+namespace grb {
+	namespace algorithms {
+
+		/**
+		 * Coloring algorithm for matrix generated by a \p DIMS - dimensional system.
+		 *
+		 * This function implements a < b>greedy heuristics< /b> to color the rows of a matrix generated by
+		 * a \p DIMS - dimensional generator \p system, so that no two connected elements \a i,j
+		 * in the system (corresponding to a nonzero \a (i,j) entry in the matrix) have the same color.
+		 * If \p reorder_rows_per_color is false (as per default), the coloring information is stored into
+		 * \p row_colors, while \p color_counters stores the number of rows for each color.
+		 *
+		 * If \p reorder_rows_per_color is true, the function performs the additional step of \b re-ordering
+		 * the rows depending on their color: rows of color \a 0 are moved first, then rows of color \a 1
+		 * are moved to the following positions and so on. In this case, \p row_colors stores the new row number
+		 * while \p color_counters stores at each position \a i the new position of the first row of color \a i.
+		 *
+		 * In both cases, \a color_counters.size() gives the number of found colors.
+		 *
+		 * This algorithm performs a \a global coloring of the input system, i.e. it must run on the entire system
+		 * \a before any partitioning occurs. Although this is not scalable, it should not be a problem for
+		 * most sizes, as the constants in front of this algorithms are very small. Implementing a distributed
+		 * coloring algorithm is anyway out of the scope of this prototype.
+		 *
+		 * Colors are by default assigned in a greedy way from the lowest one up, making this coloring scheme very
+		 * regular: close elements tend to have similar colors. This can be changed with \p lower_color_first
+		 * \c = \c false , which assigns colors from the highest one. This may avoid "destructive interference"
+		 * with following coarsening schemes.
+		 *
+		 * @tparam DIMS dimensions of the system
+		 * @tparam CoordType type of the coordinates
+		 * @tparam lower_color_first start greedy assignment of colors from lowest first
+		 *
+		 * @param[in] system generator for an \p DIMS - dimesional system with halo
+		 * @param[out] row_colors if \p reorder_rows_per_color is false, stores the color of each row;
+		 * 	if \p reorder_rows_per_color is true, stores the new position of each row, so that rows
+		 * 	of the same color are grouped together; the initial content of the vector is destroyed
+		 * @param[out] color_counters if \p reorder_rows_per_color is false, stores the number of rows per color;
+		 * 	if \p reorder_rows_per_color is true, stores at each position \a i the offset in \p color_counters
+		 * 	where the (clustered) rows of color \a i start from; the initial content of the vector is destroyed
+		 * @param[in] reorder_rows_per_color whether to do the clustering after the coloring
+		 */
+		template<
+			size_t DIMS,
+			typename CoordType,
+			bool lowest_color_first = true
+		> void hpcg_greedy_color_ndim_system(
+			const grb::utils::multigrid::LinearizedHaloNDimSystem< DIMS, CoordType > & system,
+			std::vector< CoordType > & row_colors,
+			std::vector< CoordType > & color_counters,
+			bool reorder_rows_per_color = false
+		) {
+			CoordType nrows = system.system_size();
+			// value `nrows' means `uninitialized'; initialized colors go from 0 to nrow-1
+			row_colors.insert( row_colors.begin(), nrows, nrows );
+			CoordType totalColors = 1;
+			row_colors[ 0 ] = 0; // first point gets color 0
+
+			// Finds colors in a greedy (a likely non-optimal) fashion.
+			typename grb::utils::multigrid::LinearizedHaloNDimSystem< DIMS, CoordType >::Iterator begin = system.begin();
+			begin.next_element(); // skip first row
+
+			std::vector< bool > assigned( totalColors );
+			while( begin.has_more_elements() ) {
+				CoordType curRow = begin->get_element_linear();
+
+				if( row_colors[ curRow ] != nrows ) {
+					// if color already assigned to curRow
+					continue;
+				}
+				assigned.assign( totalColors, false );
+				CoordType currentlyAssigned = 0;
+
+				while( begin.has_more_neighbours() ) {
+					CoordType curCol = begin->get_neighbor_linear();
+					if( curCol < curRow ) {
+						assert( row_colors[ curCol ] < nrows ); // if curCol < curRow, curCol has already a color assigned
+						std::vector< bool >::reference color_is_assigned = assigned[ row_colors[ curCol ] ];
+						if( ! color_is_assigned ) {
+							// count how many colors are already assigned
+							(void)currentlyAssigned++;
+						}
+						// track which colors are assigned
+						color_is_assigned = true;
+					} // else // could take advantage of indices being sorted
+					begin.next_neighbour();
+				}
+
+				if( currentlyAssigned < totalColors ) {
+					// if there is at least one color left to use, look for it
+					// smallest possible
+					if( lowest_color_first ) {
+						// here, assign colors greedily starting from the lowest available one
+						for( CoordType j = 0; j < totalColors; ++j ) {
+							if( ! assigned[ j ] ) {
+								// if no neighbor with this color, use it for this row
+								row_colors[ curRow ] = j;
+								break;
+							}
+						}
+					} else {
+						// here, assign colors greedily starting from the highest available one
+						for( CoordType j = totalColors; j > 0; --j ) {
+							CoordType color = j - 1;
+							if( ! assigned[ color ] ) {
+								// if no neighbor with this color, use it for this row
+								row_colors[ curRow ] = color;
+								break;
+							}
+						}
+					}
+				} else {
+					assert( row_colors[ curRow ] == nrows );
+					if( row_colors[ curRow ] == nrows ) {
+						row_colors[ curRow ] = totalColors;
+						(void)totalColors++;
+					} else {
+						assert( 0 ); // should never get here
+					}
+				}
+				begin.next_element();
+			}
+
+#ifdef _DEBUG
+			std::cout << "assigned colors: " << totalColors << " [ <row> -> <color>]\n";
+			for( size_t i = 0; i < row_colors.size(); i++ ) {
+				std::cout << i << " -> " << row_colors[ i ] << ", ";
+			}
+			std::cout << std::endl;
+#endif
+
+			// count number of vertices per color
+			color_counters.insert( color_counters.begin(), totalColors, 0 );
+			for( CoordType i = 0; i < nrows; ++i ) {
+				(void)color_counters[ row_colors[ i ] ]++;
+			}
+
+			if( ! reorder_rows_per_color ) {
+				return;
+			}
+
+			// form in-place prefix scan
+			CoordType old = 0, old0;
+			for( CoordType i = 1; i < totalColors; ++i ) {
+				old0 = color_counters[ i ];
+				color_counters[ i ] = color_counters[ i - 1 ] + old;
+				old = old0;
+			}
+			color_counters[ 0 ] = 0;
+
+			// translate `colors' into a permutation
+			for( CoordType i = 0; i < nrows; ++i ) {
+				row_colors[ i ] = color_counters[ row_colors[ i ] ]++;
+			}
+		}
+
+	} // namespace algorithms
+} // namespace grb
+
+#endif // _H_GRB_ALGORITHMS_HPCG_GREEDY_COLORING