diff --git a/opm/simulators/linalg/GraphColoring.hpp b/opm/simulators/linalg/GraphColoring.hpp index 001662e2a..4d43c9230 100644 --- a/opm/simulators/linalg/GraphColoring.hpp +++ b/opm/simulators/linalg/GraphColoring.hpp @@ -19,38 +19,34 @@ #ifndef OPM_GRAPHCOLORING_HEADER_INCLUDED #define OPM_GRAPHCOLORING_HEADER_INCLUDED -#include -#include -#include #include +#include +#include +#include #include #include -#include -#include +#include +#include +#include -namespace Opm -{ - -namespace Detail -{ -template +namespace Opm { +namespace Detail { +template std::size_t colorGraphWelshPowell(const Graph& graph, - std::deque& orderedVertices, - std::vector& colors, - int color, int noVertices) + std::deque& orderedVertices, + std::vector& colors, + int color, + int noVertices) { std::vector forbidden(noVertices, false); std::size_t noColored = 0; - for(auto vertex = orderedVertices.begin(), - vertexEnd = orderedVertices.end(); - vertex != vertexEnd; ++vertex) - { + for (auto vertex = orderedVertices.begin(), + vertexEnd = orderedVertices.end(); vertex != vertexEnd; ++vertex) { // Skip forbidden vertices - while(vertex != vertexEnd && forbidden[*vertex]) + while (vertex != vertexEnd && forbidden[*vertex]) ++vertex; - if ( vertex == vertexEnd ) - { + if (vertex == vertexEnd) { break; } @@ -58,25 +54,24 @@ std::size_t colorGraphWelshPowell(const Graph& graph, colors[*vertex] = color; ++noColored; // Forbid neighors - for(auto edge = graph.beginEdges(*vertex), endEdge = graph.endEdges(*vertex); - edge != endEdge; ++edge) - { + for (auto edge = graph.beginEdges(*vertex), + endEdge = graph.endEdges(*vertex); edge != endEdge; ++edge) { forbidden[edge.target()] = true; } } // forbidden vertices will be colored next for coloring using Vertex = typename Graph::VertexDescriptor; - auto newEnd = std::remove_if(orderedVertices.begin(), orderedVertices.end(), - [&forbidden](const Vertex& vertex) - { - return !forbidden[vertex]; - }); - orderedVertices.resize(newEnd-orderedVertices.begin()); + auto newEnd = std::remove_if(orderedVertices.begin(), + orderedVertices.end(), + [&forbidden](const Vertex& vertex) { return !forbidden[vertex]; }); + orderedVertices.resize(newEnd - orderedVertices.begin()); return noColored; } -template -std::size_t breadthFirstSearch(const Graph& graph, typename Graph::VertexDescriptor root, - Functor functor) + +template +std::size_t breadthFirstSearch(const Graph& graph, + typename Graph::VertexDescriptor root, + Functor functor) { std::vector visited(graph.maxVertex() + 1, false); using Vertex = typename Graph::VertexDescriptor; @@ -85,15 +80,11 @@ std::size_t breadthFirstSearch(const Graph& graph, typename Graph::VertexDescrip nextVertices.push(root); visited[root] = true; // We do not visit root. - while( !nextVertices.empty() ) - { + while (!nextVertices.empty()) { auto current = nextVertices.front(); - for(auto edge = graph.beginEdges(current), - endEdge = graph.endEdges(current); - edge != endEdge; ++edge) - { - if ( ! visited[edge.target()] ) - { + for (auto edge = graph.beginEdges(current), + endEdge = graph.endEdges(current); edge != endEdge; ++edge) { + if (!visited[edge.target()]) { visited[edge.target()] = true; nextVertices.push(edge.target()); functor(edge.target()); @@ -113,44 +104,36 @@ std::size_t breadthFirstSearch(const Graph& graph, typename Graph::VertexDescrip /// \param graph The graph to color. Must adhere to the graph interface of dune-istl. /// \return A pair of a vector with the colors of the vertices and the number of colors /// assigned -template -std::tuple, int, std::vector > +template +std::tuple, int, std::vector> colorVerticesWelshPowell(const Graph& graph) { using Vertex = typename Graph::VertexDescriptor; std::deque orderedVertices; - auto noVertices = graph.maxVertex()+1; + auto noVertices = graph.maxVertex() + 1; std::vector degrees(noVertices, 0); int maxDegree = 0; std::ptrdiff_t firstDegreeChange = 0; // populate deque - for( auto vertex = graph.begin(), endVertex = graph.end(); - vertex != endVertex; ++vertex) - { + for (auto vertex = graph.begin(), + endVertex = graph.end(); vertex != endVertex; ++vertex) { auto currentVertex = *vertex; auto& degree = degrees[currentVertex]; - for(auto edge = graph.beginEdges(currentVertex), - endEdge = graph.endEdges(currentVertex); - edge != endEdge; ++edge) - { + for (auto edge = graph.beginEdges(currentVertex), + endEdge = graph.endEdges(currentVertex); edge != endEdge; ++edge) { ++degree; } - - if( degree >= maxDegree ) - { + if (degree >= maxDegree) { orderedVertices.emplace_front(currentVertex); ++firstDegreeChange; - if(degree > maxDegree) - { + if (degree > maxDegree) { firstDegreeChange = 1; maxDegree = degree; } - } - else - { + } else { orderedVertices.emplace_back(currentVertex); } } @@ -158,10 +141,8 @@ colorVerticesWelshPowell(const Graph& graph) // order deque by descending degree std::stable_sort(orderedVertices.begin() + firstDegreeChange, orderedVertices.end(), - [°rees](const Vertex& v1, const Vertex& v2) - { - return degrees[v1] > degrees[v2]; - }); + [°rees](const Vertex& v1, const Vertex& v2) + { return degrees[v1] > degrees[v2]; }); // Overwrite degree with color auto& colors = degrees; @@ -171,39 +152,38 @@ colorVerticesWelshPowell(const Graph& graph) std::vector verticesPerColor; verticesPerColor.reserve(10); - while(!orderedVertices.empty()) - { - verticesPerColor - .push_back(Detail::colorGraphWelshPowell(graph, orderedVertices, colors, - color++, noVertices)); + while (!orderedVertices.empty()) { + verticesPerColor.push_back(Detail::colorGraphWelshPowell(graph, orderedVertices, + colors, color++, noVertices)); } return std::make_tuple(colors, color, verticesPerColor); } /// \! Reorder colored graph preserving order of vertices with the same color. -template +template std::vector -reorderVerticesPreserving(const std::vector& colors, int noColors, +reorderVerticesPreserving(const std::vector& colors, + int noColors, const std::vector& verticesPerColor, const Graph& graph) { std::vector colorIndex(noColors, 0); std::vector indices(graph.maxVertex() + 1); std::partial_sum(verticesPerColor.begin(), - verticesPerColor.begin()+verticesPerColor.size() - 1, + verticesPerColor.begin() + verticesPerColor.size() - 1, colorIndex.begin() + 1); - for(const auto& vertex: graph) - { + for (const auto& vertex : graph) { indices[vertex] = colorIndex[colors[vertex]]++; } return indices; } /// \! Reorder Vetrices in spheres -template +template std::vector -reorderVerticesSpheres(const std::vector& colors, int noColors, +reorderVerticesSpheres(const std::vector& colors, + int noColors, const std::vector& verticesPerColor, const Graph& graph, typename Graph::VertexDescriptor root) @@ -213,26 +193,22 @@ reorderVerticesSpheres(const std::vector& colors, int noColors, std::vector indices(graph.maxVertex() + 1, notVisitedTag); using Vertex = typename Graph::VertexDescriptor; std::partial_sum(verticesPerColor.begin(), - verticesPerColor.begin()+verticesPerColor.size() - 1, + verticesPerColor.begin() + verticesPerColor.size() - 1, colorIndex.begin() + 1); std::size_t noVisited = 0; auto numberer = [&colorIndex, &colors, &indices](Vertex vertex) - { - indices[vertex] = colorIndex[colors[vertex]]++; - }; - - while ( noVisited < graph.maxVertex() + 1 ) { + indices[vertex] = colorIndex[colors[vertex]]++; + }; + + while (noVisited < graph.maxVertex() + 1) { numberer(root); - ++noVisited; //root node already visited and not visited in BFS + ++noVisited; // root node already visited and not visited in BFS noVisited += Detail::breadthFirstSearch(graph, root, numberer); - if ( noVisited < graph.maxVertex() + 1 ) - { + if (noVisited < graph.maxVertex() + 1) { // Graph is disconnected search for not yet visited node - for(auto vertex: graph) - { - if ( indices[vertex] == notVisitedTag ) - { + for (auto vertex : graph) { + if (indices[vertex] == notVisitedTag) { // \todo make sure that this is a peripheral node! root = vertex; break; diff --git a/tests/test_dilu.cpp b/tests/test_dilu.cpp index e3edd2a8d..9fea1bd1b 100644 --- a/tests/test_dilu.cpp +++ b/tests/test_dilu.cpp @@ -31,8 +31,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2NoZeros, T, NumericTypes) /* Tests that the dilu decomposition mathces the expected result for a 2x2 matrix with no zero blocks, with block size 2x2. - - A + + A | | 3 1| | 1 0| | | | 2 1| | 0 1| | | | @@ -52,27 +52,27 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2NoZeros, T, NumericTypes) row.insert(1); } - A[0][0][0][0]=3.0; - A[0][0][0][1]=1.0; - A[0][0][1][0]=2.0; - A[0][0][1][1]=1.0; + A[0][0][0][0] = 3.0; + A[0][0][0][1] = 1.0; + A[0][0][1][0] = 2.0; + A[0][0][1][1] = 1.0; - A[0][1][0][0]=1.0; - A[0][1][1][1]=1.0; + A[0][1][0][0] = 1.0; + A[0][1][1][1] = 1.0; - A[1][0][0][0]=2.0; - A[1][0][1][1]=2.0; + A[1][0][0][0] = 2.0; + A[1][0][1][1] = 2.0; - A[1][1][0][0]=-1.0; - A[1][1][1][1]=-1.0; + A[1][1][0][0] = -1.0; + A[1][1][1][1] = -1.0; auto D_00 = A[0][0]; auto D_00_inv = D_00; D_00_inv.invert(); // D_11 = A_11 - L_10 D_00_inv U_01 - auto D_11 = A[1][1] - A[1][0]*D_00_inv*A[0][1]; + auto D_11 = A[1][1] - A[1][0] * D_00_inv * A[0][1]; Dune::SeqDilu seqdilu(A); @@ -99,8 +99,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2, T, NumericTypes) /* Tests that the dilu decomposition mathces the expected result for a 2x2 matrix, with block size 2x2. - - A + + A | | 3 1| | 1 0| | | | 2 1| | 0 1| | | | @@ -122,16 +122,16 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2, T, NumericTypes) } } - A[0][0][0][0]=3.0; - A[0][0][0][1]=1.0; - A[0][0][1][0]=2.0; - A[0][0][1][1]=1.0; + A[0][0][0][0] = 3.0; + A[0][0][0][1] = 1.0; + A[0][0][1][0] = 2.0; + A[0][0][1][1] = 1.0; - A[0][1][0][0]=1.0; - A[0][1][1][1]=1.0; + A[0][1][0][0] = 1.0; + A[0][1][1][1] = 1.0; - A[1][1][0][0]=-1.0; - A[1][1][1][1]=-1.0; + A[1][1][0][0] = -1.0; + A[1][1][1][1] = -1.0; auto D_00 = A[0][0]; @@ -188,20 +188,20 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrectNoZeros, T, NumericTypes) row.insert(1); } - A[0][0][0][0]=3.0; - A[0][0][0][1]=1.0; - A[0][0][1][0]=2.0; - A[0][0][1][1]=1.0; + A[0][0][0][0] = 3.0; + A[0][0][0][1] = 1.0; + A[0][0][1][0] = 2.0; + A[0][0][1][1] = 1.0; - A[0][1][0][0]=1.0; - A[0][1][1][1]=1.0; + A[0][1][0][0] = 1.0; + A[0][1][1][1] = 1.0; - A[1][0][0][0]=2.0; - A[1][0][1][1]=2.0; + A[1][0][0][0] = 2.0; + A[1][0][1][1] = 2.0; - A[1][1][0][0]=-1.0; - A[1][1][1][1]=-1.0; + A[1][1][0][0] = -1.0; + A[1][1][1][1] = -1.0; Vector x(2); @@ -221,9 +221,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrectNoZeros, T, NumericTypes) auto D_00_inv = D_00; D_00_inv.invert(); // D_11= A_11 - L_10 D_00_inv U_01 - auto D_11 = A[1][1] - A[1][0]*D_00_inv*A[0][1]; + auto D_11 = A[1][1] - A[1][0] * D_00_inv * A[0][1]; auto D_11_inv = D_11; - D_11_inv.invert(); + D_11_inv.invert(); // define: z = M^-1(b - Ax) // where: M = (D + L_A) A^-1 (D + U_A) @@ -251,8 +251,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrectNoZeros, T, NumericTypes) // z_0 = y_0 - D_00_inv*A_01*z_1 z[0] = y[0]; - auto temp = D_00_inv*A[0][1]; - temp.mmv(z[1], z[0]); + auto temp = D_00_inv * A[0][1]; + temp.mmv(z[1], z[0]); // x_k+1 = x_k + z Vector new_x = x; @@ -299,16 +299,16 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect1, T, NumericTypes) } } - A[0][0][0][0]=3.0; - A[0][0][0][1]=1.0; - A[0][0][1][0]=2.0; - A[0][0][1][1]=1.0; + A[0][0][0][0] = 3.0; + A[0][0][0][1] = 1.0; + A[0][0][1][0] = 2.0; + A[0][0][1][1] = 1.0; - A[0][1][0][0]=1.0; - A[0][1][1][1]=1.0; + A[0][1][0][0] = 1.0; + A[0][1][1][1] = 1.0; - A[1][1][0][0]=-1.0; - A[1][1][1][1]=-1.0; + A[1][1][0][0] = -1.0; + A[1][1][1][1] = -1.0; Vector x(2); @@ -330,7 +330,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect1, T, NumericTypes) // D_11 = A_11 - L_10 D_0_inv U_01 = A_11 auto D_11 = A[1][1]; auto D_11_inv = D_11; - D_11_inv.invert(); + D_11_inv.invert(); // define: z = M^-1(b - Ax) // where: M = (D + L_A) A^-1 (D + U_A) @@ -356,8 +356,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect1, T, NumericTypes) // z_0 = y_0 - D_00_inv*A_01*z_1 z[0] = y[0]; - auto temp = D_00_inv*A[0][1]; - temp.mmv(z[1], z[0]); + auto temp = D_00_inv * A[0][1]; + temp.mmv(z[1], z[0]); // x_k+1 = x_k + z Vector new_x = x; @@ -402,16 +402,16 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect2, T, NumericTypes) } } - A[0][0][0][0]=3.0; - A[0][0][0][1]=1.0; - A[0][0][1][0]=2.0; - A[0][0][1][1]=1.0; + A[0][0][0][0] = 3.0; + A[0][0][0][1] = 1.0; + A[0][0][1][0] = 2.0; + A[0][0][1][1] = 1.0; - A[1][1][0][0]=2.0; - A[1][1][1][1]=2.0; + A[1][1][0][0] = 2.0; + A[1][1][1][1] = 2.0; - A[1][1][0][0]=-1.0; - A[1][1][1][1]=-1.0; + A[1][1][0][0] = -1.0; + A[1][1][1][1] = -1.0; Vector x(2); x[0][0] = 1.0; @@ -432,7 +432,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect2, T, NumericTypes) // D_11 = A_11 - L_10 D_0_inv U_01 = A_11 auto D_11 = A[1][1]; auto D_11_inv = D_11; - D_11_inv.invert(); + D_11_inv.invert(); // define: z = M^-1(b - Ax) // where: M = (D + L_A) A^-1 (D + U_A) @@ -479,13 +479,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect3x3, T, NumericTypes) Tests that the dilu decomposition mathces the expected result for a 3x3 matrix, with block size 3x3. - A + A | | 3 1 2| | 0 0 0| | 0 0 0| | | | 2 3 1| | 0 0 0| | 0 0 0| | | | 2 1 0| | 0 0 0| | 0 0 0| | | | | | 0 0 0| | 1 0 1| | 1 0 2| | - | | 0 0 0| | 4 1 0| | 0 1 1| | + | | 0 0 0| | 4 1 0| | 0 1 1| | | | 0 0 0| | 3 1 3| | 0 1 3| | | | | | 0 0 0| | 1 0 2| | 1 3 2| | @@ -504,36 +504,61 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect3x3, T, NumericTypes) for (auto row = A.createbegin(); row != A.createend(); ++row) { if (row.index() == 0) { row.insert(row.index()); - } - else if (row.index() == 1) { + } else if (row.index() == 1) { row.insert(row.index()); row.insert(row.index() + 1); - } - else if (row.index() == 2) { + } else if (row.index() == 2) { row.insert(row.index() - 1); row.insert(row.index()); } } - A[0][0][0][0]=3.0; A[1][1][0][0]=1.0; A[1][2][0][0]=1.0; - A[0][0][0][1]=1.0; A[1][1][0][1]=0.0; A[1][2][0][1]=0.0; - A[0][0][0][2]=2.0; A[1][1][0][2]=1.0; A[1][2][0][2]=2.0; - A[0][0][1][0]=2.0; A[1][1][1][0]=4.0; A[1][2][1][0]=0.0; - A[0][0][1][1]=3.0; A[1][1][1][1]=1.0; A[1][2][1][1]=1.0; - A[0][0][1][2]=1.0; A[1][1][1][2]=0.0; A[1][2][1][2]=1.0; - A[0][0][2][0]=2.0; A[1][1][2][0]=3.0; A[1][2][2][0]=0.0; - A[0][0][2][1]=1.0; A[1][1][2][1]=1.0; A[1][2][2][1]=1.0; - A[0][0][2][2]=0.0; A[1][1][2][2]=3.0; A[1][2][2][2]=3.0; + A[0][0][0][0] = 3.0; + A[1][1][0][0] = 1.0; + A[1][2][0][0] = 1.0; + A[0][0][0][1] = 1.0; + A[1][1][0][1] = 0.0; + A[1][2][0][1] = 0.0; + A[0][0][0][2] = 2.0; + A[1][1][0][2] = 1.0; + A[1][2][0][2] = 2.0; + A[0][0][1][0] = 2.0; + A[1][1][1][0] = 4.0; + A[1][2][1][0] = 0.0; + A[0][0][1][1] = 3.0; + A[1][1][1][1] = 1.0; + A[1][2][1][1] = 1.0; + A[0][0][1][2] = 1.0; + A[1][1][1][2] = 0.0; + A[1][2][1][2] = 1.0; + A[0][0][2][0] = 2.0; + A[1][1][2][0] = 3.0; + A[1][2][2][0] = 0.0; + A[0][0][2][1] = 1.0; + A[1][1][2][1] = 1.0; + A[1][2][2][1] = 1.0; + A[0][0][2][2] = 0.0; + A[1][1][2][2] = 3.0; + A[1][2][2][2] = 3.0; - A[2][1][0][0]=1.0; A[2][2][0][0]=1.0; - A[2][1][0][1]=0.0; A[2][2][0][1]=3.0; - A[2][1][0][2]=2.0; A[2][2][0][2]=2.0; - A[2][1][1][0]=0.0; A[2][2][1][0]=2.0; - A[2][1][1][1]=1.0; A[2][2][1][1]=1.0; - A[2][1][1][2]=4.0; A[2][2][1][2]=3.0; - A[2][1][2][0]=5.0; A[2][2][2][0]=3.0; - A[2][1][2][1]=1.0; A[2][2][2][1]=1.0; - A[2][1][2][2]=1.0; A[2][2][2][2]=2.0; + A[2][1][0][0] = 1.0; + A[2][2][0][0] = 1.0; + A[2][1][0][1] = 0.0; + A[2][2][0][1] = 3.0; + A[2][1][0][2] = 2.0; + A[2][2][0][2] = 2.0; + A[2][1][1][0] = 0.0; + A[2][2][1][0] = 2.0; + A[2][1][1][1] = 1.0; + A[2][2][1][1] = 1.0; + A[2][1][1][2] = 4.0; + A[2][2][1][2] = 3.0; + A[2][1][2][0] = 5.0; + A[2][2][2][0] = 3.0; + A[2][1][2][1] = 1.0; + A[2][2][2][1] = 1.0; + A[2][1][2][2] = 1.0; + A[2][2][2][2] = 2.0; auto D_00 = A[0][0]; @@ -542,11 +567,11 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect3x3, T, NumericTypes) // D_11 = A_11 - L_10 D_00_inv U_01 = A_11 auto D_11 = A[1][1]; auto D_11_inv = D_11; - D_11_inv.invert(); + D_11_inv.invert(); // D_22 = A_22 - A_20 D_00_inv A_02 - A_21 D_11_inv A_12 = A_22 - A_21 D_11_inv A_12 - auto D_22 = A[2][2] - A[2][1]*D_11_inv*A[1][2]; + auto D_22 = A[2][2] - A[2][1] * D_11_inv * A[1][2]; auto D_22_inv = D_22; - D_22_inv.invert(); + D_22_inv.invert(); Dune::SeqDilu seqdilu(A); auto Dinv = seqdilu.getDiagonal(); @@ -583,13 +608,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect3, T, NumericTypes) | | 2 1 0| | 0 0 0| | 0 0 0| | | |3| | | |2| | | | | | | | | | 0 0 0| | 1 0 1| | 1 0 2| | | |1| | | |2| | - | | 0 0 0| | 4 1 0| | 0 1 1| | x | |3| | = | |3| | + | | 0 0 0| | 4 1 0| | 0 1 1| | x | |3| | = | |3| | | | 0 0 0| | 3 1 3| | 0 1 3| | | |2| | | |2| | | | | | | | | | 0 0 0| | 1 0 2| | 1 3 2| | | |1| | | |0| | | | 0 0 0| | 0 1 4| | 2 1 3| | | |0| | | |2| | | | 0 0 0| | 5 1 1| | 3 1 2| | | |2| | | |1| | - + */ const int N = 3; @@ -602,46 +627,83 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect3, T, NumericTypes) for (auto row = A.createbegin(); row != A.createend(); ++row) { if (row.index() == 0) { row.insert(row.index()); - } - else if (row.index() == 1) { + } else if (row.index() == 1) { row.insert(row.index()); row.insert(row.index() + 1); - } - else if (row.index() == 2) { + } else if (row.index() == 2) { row.insert(row.index() - 1); row.insert(row.index()); } } - A[0][0][0][0]=3.0; A[1][1][0][0]=1.0; A[1][2][0][0]=1.0; - A[0][0][0][1]=1.0; A[1][1][0][1]=0.0; A[1][2][0][1]=0.0; - A[0][0][0][2]=2.0; A[1][1][0][2]=1.0; A[1][2][0][2]=2.0; - A[0][0][1][0]=2.0; A[1][1][1][0]=4.0; A[1][2][1][0]=0.0; - A[0][0][1][1]=3.0; A[1][1][1][1]=1.0; A[1][2][1][1]=1.0; - A[0][0][1][2]=1.0; A[1][1][1][2]=0.0; A[1][2][1][2]=1.0; - A[0][0][2][0]=2.0; A[1][1][2][0]=3.0; A[1][2][2][0]=0.0; - A[0][0][2][1]=1.0; A[1][1][2][1]=1.0; A[1][2][2][1]=1.0; - A[0][0][2][2]=0.0; A[1][1][2][2]=3.0; A[1][2][2][2]=3.0; + A[0][0][0][0] = 3.0; + A[1][1][0][0] = 1.0; + A[1][2][0][0] = 1.0; + A[0][0][0][1] = 1.0; + A[1][1][0][1] = 0.0; + A[1][2][0][1] = 0.0; + A[0][0][0][2] = 2.0; + A[1][1][0][2] = 1.0; + A[1][2][0][2] = 2.0; + A[0][0][1][0] = 2.0; + A[1][1][1][0] = 4.0; + A[1][2][1][0] = 0.0; + A[0][0][1][1] = 3.0; + A[1][1][1][1] = 1.0; + A[1][2][1][1] = 1.0; + A[0][0][1][2] = 1.0; + A[1][1][1][2] = 0.0; + A[1][2][1][2] = 1.0; + A[0][0][2][0] = 2.0; + A[1][1][2][0] = 3.0; + A[1][2][2][0] = 0.0; + A[0][0][2][1] = 1.0; + A[1][1][2][1] = 1.0; + A[1][2][2][1] = 1.0; + A[0][0][2][2] = 0.0; + A[1][1][2][2] = 3.0; + A[1][2][2][2] = 3.0; - A[2][1][0][0]=1.0; A[2][2][0][0]=1.0; - A[2][1][0][1]=0.0; A[2][2][0][1]=3.0; - A[2][1][0][2]=2.0; A[2][2][0][2]=2.0; - A[2][1][1][0]=0.0; A[2][2][1][0]=2.0; - A[2][1][1][1]=1.0; A[2][2][1][1]=1.0; - A[2][1][1][2]=4.0; A[2][2][1][2]=3.0; - A[2][1][2][0]=5.0; A[2][2][2][0]=3.0; - A[2][1][2][1]=1.0; A[2][2][2][1]=1.0; - A[2][1][2][2]=1.0; A[2][2][2][2]=2.0; + A[2][1][0][0] = 1.0; + A[2][2][0][0] = 1.0; + A[2][1][0][1] = 0.0; + A[2][2][0][1] = 3.0; + A[2][1][0][2] = 2.0; + A[2][2][0][2] = 2.0; + A[2][1][1][0] = 0.0; + A[2][2][1][0] = 2.0; + A[2][1][1][1] = 1.0; + A[2][2][1][1] = 1.0; + A[2][1][1][2] = 4.0; + A[2][2][1][2] = 3.0; + A[2][1][2][0] = 5.0; + A[2][2][2][0] = 3.0; + A[2][1][2][1] = 1.0; + A[2][2][2][1] = 1.0; + A[2][1][2][2] = 1.0; + A[2][2][2][2] = 2.0; Vector x(3); - x[0][0] = 1.0; x[1][0] = 1.0; x[2][0] = 1.0; - x[0][1] = 2.0; x[1][1] = 3.0; x[2][1] = 0.0; - x[0][2] = 3.0; x[1][2] = 2.0; x[2][2] = 2.0; + x[0][0] = 1.0; + x[1][0] = 1.0; + x[2][0] = 1.0; + x[0][1] = 2.0; + x[1][1] = 3.0; + x[2][1] = 0.0; + x[0][2] = 3.0; + x[1][2] = 2.0; + x[2][2] = 2.0; Vector b(3); - b[0][0] = 2.0; b[1][0] = 2.0; b[2][0] = 0.0; - b[0][1] = 1.0; b[1][1] = 3.0; b[2][1] = 2.0; - b[0][2] = 2.0; b[1][2] = 2.0; b[2][2] = 1.0; + b[0][0] = 2.0; + b[1][0] = 2.0; + b[2][0] = 0.0; + b[0][1] = 1.0; + b[1][1] = 3.0; + b[2][1] = 2.0; + b[0][2] = 2.0; + b[1][2] = 2.0; + b[2][2] = 1.0; // D_00 = A_00 @@ -652,17 +714,17 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect3, T, NumericTypes) // = A_11 auto D_11 = A[1][1]; auto D_11_inv = D_11; - D_11_inv.invert(); + D_11_inv.invert(); // D_22 = A_22 - A_20 D_00_inv A_02 - A_21 D_11_inv A_12 // = A_22 - A_21 D_11_inv A_12 - auto D_22 = A[2][2] - A[2][1]*D_11_inv*A[1][2]; + auto D_22 = A[2][2] - A[2][1] * D_11_inv * A[1][2]; auto D_22_inv = D_22; - D_22_inv.invert(); + D_22_inv.invert(); // define: z = M^-1(b - Ax) // where: M = (D + L_A) A^-1 (D + U_A) // lower triangular solve: (E + L) y = b - Ax - + Vector y(3); // y_0 = D_00_inv*[b_0 - (A_00*x_0 + A_01*x_1)] // = D_00_inv*[b_0 - A_00*x_0] @@ -693,8 +755,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect3, T, NumericTypes) // z_1 = y_1 - D_11_inv*A_12*z_2 z[1] = y[1]; - auto temp = D_11_inv*A[1][2]; - temp.mmv(z[2], z[1]); + auto temp = D_11_inv * A[1][2]; + temp.mmv(z[2], z[1]); // z_0 = y_0 - D_00_inv(A_01*z_1 + A_02*z_2) // z_0 = y_0 @@ -741,13 +803,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsEqualToDuneSeqILUApply, T, NumericTy row.insert(row.index()); } - A[0][0][0][0]=3.0; - A[0][0][0][1]=1.0; - A[0][0][1][0]=2.0; - A[0][0][1][1]=1.0; + A[0][0][0][0] = 3.0; + A[0][0][0][1] = 1.0; + A[0][0][1][0] = 2.0; + A[0][0][1][1] = 1.0; - A[1][1][0][0]=-1.0; - A[1][1][1][1]=-1.0; + A[1][1][0][0] = -1.0; + A[1][1][1][1] = -1.0; Dune::SeqDilu seqdilu(A);