diff --git a/opm/simulators/linalg/GraphColoring.hpp b/opm/simulators/linalg/GraphColoring.hpp index 4d43c9230..047168039 100644 --- a/opm/simulators/linalg/GraphColoring.hpp +++ b/opm/simulators/linalg/GraphColoring.hpp @@ -19,6 +19,10 @@ #ifndef OPM_GRAPHCOLORING_HEADER_INCLUDED #define OPM_GRAPHCOLORING_HEADER_INCLUDED +#include + +#include + #include #include #include @@ -218,5 +222,94 @@ reorderVerticesSpheres(const std::vector& colors, } return indices; } -} // end namespace Opm + +/// \brief Specify coloring type. +/// \details The coloring types have been implemented initially to parallelize DILU +/// preconditioner and parallel sparse triangular solves. +/// Symmetric coloring will create a dependency from row i to j if +/// both element A_ij and A_ji exists. +/// Lower coloring creates a dependency from row i to j (where i < j) if +/// A_ij is nonzero. +/// Upper coloring creates a dependecy from row i to j (where i > j) if A_ij is nonzero. +enum class ColoringType { SYMMETRIC, LOWER, UPPER }; + +/// This coloring algorithm interprets the sparsity structure of a matrix as a graph. Each row is given a color +/// or level where all the rows in the same level only have dependencies from lower levels. The level computation +/// is done with dynamic programming, and to improve caching the rows in the same level stay in matrix order. +/// \brief Given a matrix and dependecy type, returns a SparseTable grouping the rows by which +/// can be executed in parallel without breaking dependencies +/// \param matrix A dune sparse matrix +/// \param coloringType The coloringtype determines what constitutes a dependency, +/// see ColoringType definition above +/// \return SparseTable with rows of the matrix grouped into least number of groups +/// while dependencies only come from groups with lower index +template +Opm::SparseTable +getMatrixRowColoring(const M& matrix, ColoringType coloringType) +{ + OPM_TIMEBLOCK(createMatrix); + + std::vector color(matrix.N(), 0); + std::vector rowIndices(matrix.N(), 0); + +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (std::size_t i = 0; i < matrix.N(); i++) { + rowIndices[i] = i; + } + + std::vector colorCnt; + // These dynamic programming computations only rely on the following observation: + // level[row_i] = 1 + max{level[row_j]} for all j which i depend on + // This minimizes the level of each row, and every rows dependencies belong to a lower level set. + if (coloringType == ColoringType::SYMMETRIC) { + for (auto i = matrix.begin(); i != matrix.end(); ++i) { + for (auto a_ij = i->begin(); a_ij.index() != i.index(); ++a_ij) { + auto a_ji = matrix[a_ij.index()].find(i.index()); + if (a_ji != matrix[a_ij.index()].end()) { + color[i.index()] = std::max({color[i.index()], color[a_ij.index()] + 1}); + } + } + if (color[i.index()] >= colorCnt.size()) { + colorCnt.push_back(1); + } else { + ++colorCnt[color[i.index()]]; + } + } + } else if (coloringType == ColoringType::UPPER) { + for (auto i = matrix.beforeEnd(); i != matrix.beforeBegin(); --i) { + for (auto a_ij = ++(i->find(i.index())); a_ij != i->end(); ++a_ij) { + color[i.index()] = std::max({color[i.index()], color[a_ij.index()] + 1}); + } + if (color[i.index()] >= colorCnt.size()) { + colorCnt.push_back(1); + } else { + ++colorCnt[color[i.index()]]; + } + } + } else if (coloringType == ColoringType::LOWER) { + for (auto i = matrix.begin(); i != matrix.end(); ++i) { + for (auto a_ij = i->begin(); a_ij.index() != i.index(); ++a_ij) { + color[i.index()] = std::max({color[i.index()], color[a_ij.index()] + 1}); + } + if (color[i.index()] >= colorCnt.size()) { + colorCnt.push_back(1); + } else { + ++colorCnt[color[i.index()]]; + } + } + } + + std::stable_sort(rowIndices.begin(), + rowIndices.end(), + [&color](const std::size_t a, const std::size_t b) + { return color[a] < color[b]; }); + + return {rowIndices.data(), rowIndices.data() + rowIndices.size(), + colorCnt.data(), colorCnt.data() + colorCnt.size()}; +} + +} // end namespace Opm + #endif diff --git a/tests/test_graphcoloring.cpp b/tests/test_graphcoloring.cpp index c05a60741..e56ca8202 100644 --- a/tests/test_graphcoloring.cpp +++ b/tests/test_graphcoloring.cpp @@ -99,3 +99,323 @@ BOOST_AUTO_TEST_CASE(TestWelschPowell) graph, 0); checkAllIndices(newOrder); } + +// The following tests verify the graph coloring in the context of revealing which rows +// can be operated on at the same time in the DILU preconditioner +BOOST_AUTO_TEST_CASE(TestColoredDiluParallelisms3x3Matrix) +{ + /* + Matrix on this form: + |x | + | xx| + | xx| + We only expect a DILU dependency from the second to the third row, + hence row 1 and 2 should have color 0, row 3 should have color 1 + */ + const int N = 3; + const int bz = 3; + const int nonZeroes = 5; + + // creating some shorthand typenames + using Matrix = Dune::BCRSMatrix>; + + Matrix testMatrix(N, N, nonZeroes, Matrix::row_wise); + for (auto row = testMatrix.createbegin(); row != testMatrix.createend(); ++row) { + if (row.index() == 0) { + row.insert(row.index()); + } + else if (row.index() == 1) { + row.insert(row.index()); + row.insert(row.index() + 1); + } + else if (row.index() == 2) { + row.insert(row.index() - 1); + row.insert(row.index()); + } + } + + testMatrix[0][0][0][0] = 3.0; + testMatrix[1][1][0][0] = 1.0; + testMatrix[1][2][0][0] = 1.0; + testMatrix[2][1][0][0] = 1.0; + testMatrix[2][2][0][0] = 1.0; + + Opm::SparseTable coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::SYMMETRIC); + + std::vector> correctColor = {{0, 1}, {2}}; + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } + + coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::UPPER); + correctColor = {{0, 2}, {1}}; + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } + coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::LOWER); + correctColor = {{0, 1}, {2}}; + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } +} + +BOOST_AUTO_TEST_CASE(TestColoredDiluParallelisms5x5Simple) +{ + /* + Test matrix: + |xxx | + |xx | + |x xx | + | x | + | xx| + */ + const int N = 5; + const int bz = 3; + const int nonZeroes = 11; + + // creating some shorthand typenames + using Matrix = Dune::BCRSMatrix>; + + Matrix testMatrix(N, N, nonZeroes, Matrix::row_wise); + for (auto row = testMatrix.createbegin(); row != testMatrix.createend(); ++row) { + if (row.index() == 0) { + row.insert(row.index()); + row.insert(row.index()+1); + row.insert(row.index()+2); + } + else if (row.index() == 1) { + row.insert(row.index()); + row.insert(row.index() - 1); + } + else if (row.index() == 2) { + row.insert(row.index() - 2); + row.insert(row.index()); + row.insert(row.index() + 1); + } + else if (row.index() == 3) { + row.insert(row.index()); + } + else if (row.index() == 4) { + row.insert(row.index() - 1); + row.insert(row.index()); + } + } + + testMatrix[0][0][0][0] = 1.0; + testMatrix[0][1][0][0] = 1.0; + testMatrix[0][2][0][0] = 1.0; + testMatrix[1][0][0][0] = 1.0; + testMatrix[1][1][0][0] = 1.0; + testMatrix[2][0][0][0] = 1.0; + testMatrix[2][2][0][0] = 1.0; + testMatrix[2][3][0][0] = 1.0; + testMatrix[3][3][0][0] = 1.0; + testMatrix[4][3][0][0] = 1.0; + testMatrix[4][4][0][0] = 1.0; + + Opm::SparseTable coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::SYMMETRIC); + + std::vector> correctColor = {{0, 3, 4}, {1, 2}}; + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } + + coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::UPPER); + correctColor = {{1, 3, 4}, {2}, {0}}; + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } + + coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::LOWER); + correctColor = {{0, 3}, {1, 2, 4}}; + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } +} + +BOOST_AUTO_TEST_CASE(TestColoredDiluParallelisms5x5Tridiagonal) +{ + /* + Test matrix: + |xx | + |xxx | + | xxx | + | xxx| + | xx| + The tridiagonal structure will force a strictly serial computation stage + */ + const int N = 5; + const int bz = 3; + const int nonZeroes = 13; + + // creating some shorthand typenames + using Matrix = Dune::BCRSMatrix>; + + Matrix testMatrix(N, N, nonZeroes, Matrix::row_wise); + for (auto row = testMatrix.createbegin(); row != testMatrix.createend(); ++row) { + if (row.index() == 0) { + row.insert(row.index()); + row.insert(row.index()+1); + } + else if (row.index() > 0 && row.index() < 4) { + row.insert(row.index() - 1); + row.insert(row.index()); + row.insert(row.index() + 1); + } + else if (row.index() == 4) { + row.insert(row.index() - 1); + row.insert(row.index()); + } + } + + testMatrix[0][0][0][0] = 1.0; + testMatrix[0][1][0][0] = 1.0; + testMatrix[1][0][0][0] = 1.0; + testMatrix[1][1][0][0] = 1.0; + testMatrix[1][2][0][0] = 1.0; + testMatrix[2][1][0][0] = 1.0; + testMatrix[2][2][0][0] = 1.0; + testMatrix[2][3][0][0] = 1.0; + testMatrix[3][2][0][0] = 1.0; + testMatrix[3][3][0][0] = 1.0; + testMatrix[3][4][0][0] = 1.0; + testMatrix[4][3][0][0] = 1.0; + testMatrix[4][4][0][0] = 1.0; + + Opm::SparseTable coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::SYMMETRIC); + + std::vector> correctColor = {{0}, {1}, {2}, {3}, {4}}; + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } + + coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::LOWER); + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } + + + coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::UPPER); + correctColor = {{4}, {3}, {2}, {1}, {0}}; + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } +} + +BOOST_AUTO_TEST_CASE(TestColoredDiluParallelisms5x5Complex) +{ + /* + Test matrix: + |xxx x| + |xx x | + |x x x| + | x x | + |x x x| + */ + const int N = 5; + const int bz = 3; + const int nonZeroes = 15; + + // creating some shorthand typenames + using Matrix = Dune::BCRSMatrix>; + + Matrix testMatrix(N, N, nonZeroes, Matrix::row_wise); + for (auto row = testMatrix.createbegin(); row != testMatrix.createend(); ++row) { + if (row.index() == 0) { + row.insert(row.index()); + row.insert(row.index()+1); + row.insert(row.index()+2); + row.insert(row.index()+4); + } + else if (row.index() == 1) { + row.insert(row.index() - 1); + row.insert(row.index()); + row.insert(row.index() + 2); + } + else if (row.index() == 2) { + row.insert(row.index() - 2); + row.insert(row.index()); + row.insert(row.index() + 2); + } + else if (row.index() == 3) { + row.insert(row.index() - 2); + row.insert(row.index()); + } + else if (row.index() == 4) { + row.insert(row.index() - 4); + row.insert(row.index() - 2); + row.insert(row.index()); + } + } + + testMatrix[0][0][0][0] = 1.0; + testMatrix[0][1][0][0] = 1.0; + testMatrix[0][2][0][0] = 1.0; + testMatrix[0][4][0][0] = 1.0; + testMatrix[1][0][0][0] = 1.0; + testMatrix[1][1][0][0] = 1.0; + testMatrix[1][3][0][0] = 1.0; + testMatrix[2][0][0][0] = 1.0; + testMatrix[2][2][0][0] = 1.0; + testMatrix[2][4][0][0] = 1.0; + testMatrix[3][1][0][0] = 1.0; + testMatrix[3][3][0][0] = 1.0; + testMatrix[4][0][0][0] = 1.0; + testMatrix[4][2][0][0] = 1.0; + testMatrix[4][4][0][0] = 1.0; + + Opm::SparseTable coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::SYMMETRIC); + + std::vector> correctColor = {{0}, {1, 2}, {3, 4}}; + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } + + coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::LOWER); + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } + + coloring = Opm::getMatrixRowColoring(testMatrix, Opm::ColoringType::UPPER); + correctColor = {{3, 4}, {1, 2}, {0}}; + coloring.print(std::cout); + + for (std::size_t i = 0; i < correctColor.size(); ++i){ + for (std::size_t j = 0; j < correctColor[i].size(); ++j){ + BOOST_CHECK(coloring[i][j] == correctColor[i][j]); + } + } +} \ No newline at end of file