OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5006 from akva2/janitoring

Browse Source 
Some janitoring
This commit is contained in:
Arne Morten Kvarving 2023-11-20 14:29:45 +01:00 committed by GitHub
commit 4ccc5cddc7
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 252 additions and 214 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

152
opm/simulators/linalg/GraphColoring.hpp
View File

@ -19,38 +19,34 @@
#ifndef OPM_GRAPHCOLORING_HEADER_INCLUDED #ifndef OPM_GRAPHCOLORING_HEADER_INCLUDED
#define OPM_GRAPHCOLORING_HEADER_INCLUDED #define OPM_GRAPHCOLORING_HEADER_INCLUDED
#include <vector>
#include <deque>
#include <tuple>
#include <algorithm> #include <algorithm>
#include <cstddef>
#include <deque>
#include <limits>
#include <numeric> #include <numeric>
#include <queue> #include <queue>
#include <cstddef> #include <string>
#include <limits> #include <tuple>
#include <vector>
namespace Opm namespace Opm {
{ namespace Detail {
template <class Graph>
namespace Detail
{
template<class Graph>
std::size_t colorGraphWelshPowell(const Graph& graph, std::size_t colorGraphWelshPowell(const Graph& graph,
std::deque<typename Graph::VertexDescriptor>& orderedVertices, std::deque<typename Graph::VertexDescriptor>& orderedVertices,
std::vector<int>& colors, std::vector<int>& colors,
int color, int noVertices) int color,
int noVertices)
{ {
std::vector<int> forbidden(noVertices, false); std::vector<int> forbidden(noVertices, false);
std::size_t noColored = 0; std::size_t noColored = 0;
for(auto vertex = orderedVertices.begin(), for (auto vertex = orderedVertices.begin(),
vertexEnd = orderedVertices.end(); vertexEnd = orderedVertices.end(); vertex != vertexEnd; ++vertex) {
vertex != vertexEnd; ++vertex)
{
// Skip forbidden vertices // Skip forbidden vertices
while(vertex != vertexEnd && forbidden[*vertex]) while (vertex != vertexEnd && forbidden[*vertex])
++vertex; ++vertex;
if ( vertex == vertexEnd ) if (vertex == vertexEnd) {
{
break; break;
} }
@ -58,25 +54,24 @@ std::size_t colorGraphWelshPowell(const Graph& graph,
colors[*vertex] = color; colors[*vertex] = color;
++noColored; ++noColored;
// Forbid neighors // Forbid neighors
for(auto edge = graph.beginEdges(*vertex), endEdge = graph.endEdges(*vertex); for (auto edge = graph.beginEdges(*vertex),
edge != endEdge; ++edge) endEdge = graph.endEdges(*vertex); edge != endEdge; ++edge) {
{
forbidden[edge.target()] = true; forbidden[edge.target()] = true;
} }
} }
// forbidden vertices will be colored next for coloring // forbidden vertices will be colored next for coloring
using Vertex = typename Graph::VertexDescriptor; using Vertex = typename Graph::VertexDescriptor;
auto newEnd = std::remove_if(orderedVertices.begin(), orderedVertices.end(), auto newEnd = std::remove_if(orderedVertices.begin(),
[&forbidden](const Vertex& vertex) orderedVertices.end(),
{ [&forbidden](const Vertex& vertex) { return !forbidden[vertex]; });
return !forbidden[vertex]; orderedVertices.resize(newEnd - orderedVertices.begin());
});
orderedVertices.resize(newEnd-orderedVertices.begin());
return noColored; return noColored;
} }
template<class Graph, class Functor>
std::size_t breadthFirstSearch(const Graph& graph, typename Graph::VertexDescriptor root, template <class Graph, class Functor>
Functor functor) std::size_t breadthFirstSearch(const Graph& graph,
typename Graph::VertexDescriptor root,
Functor functor)
{ {
std::vector<int> visited(graph.maxVertex() + 1, false); std::vector<int> visited(graph.maxVertex() + 1, false);
using Vertex = typename Graph::VertexDescriptor; using Vertex = typename Graph::VertexDescriptor;
@ -85,15 +80,11 @@ std::size_t breadthFirstSearch(const Graph& graph, typename Graph::VertexDescrip
nextVertices.push(root); nextVertices.push(root);
visited[root] = true; // We do not visit root. visited[root] = true; // We do not visit root.
while( !nextVertices.empty() ) while (!nextVertices.empty()) {
{
auto current = nextVertices.front(); auto current = nextVertices.front();
for(auto edge = graph.beginEdges(current), for (auto edge = graph.beginEdges(current),
endEdge = graph.endEdges(current); endEdge = graph.endEdges(current); edge != endEdge; ++edge) {
edge != endEdge; ++edge) if (!visited[edge.target()]) {
{
if ( ! visited[edge.target()] )
{
visited[edge.target()] = true; visited[edge.target()] = true;
nextVertices.push(edge.target()); nextVertices.push(edge.target());
functor(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. /// \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 /// \return A pair of a vector with the colors of the vertices and the number of colors
/// assigned /// assigned
template<class Graph> template <class Graph>
std::tuple<std::vector<int>, int, std::vector<std::size_t> > std::tuple<std::vector<int>, int, std::vector<std::size_t>>
colorVerticesWelshPowell(const Graph& graph) colorVerticesWelshPowell(const Graph& graph)
{ {
using Vertex = typename Graph::VertexDescriptor; using Vertex = typename Graph::VertexDescriptor;
std::deque<Vertex> orderedVertices; std::deque<Vertex> orderedVertices;
auto noVertices = graph.maxVertex()+1; auto noVertices = graph.maxVertex() + 1;
std::vector<int> degrees(noVertices, 0); std::vector<int> degrees(noVertices, 0);
int maxDegree = 0; int maxDegree = 0;
std::ptrdiff_t firstDegreeChange = 0; std::ptrdiff_t firstDegreeChange = 0;
// populate deque // populate deque
for( auto vertex = graph.begin(), endVertex = graph.end(); for (auto vertex = graph.begin(),
vertex != endVertex; ++vertex) endVertex = graph.end(); vertex != endVertex; ++vertex) {
{
auto currentVertex = *vertex; auto currentVertex = *vertex;
auto& degree = degrees[currentVertex]; auto& degree = degrees[currentVertex];
for(auto edge = graph.beginEdges(currentVertex), for (auto edge = graph.beginEdges(currentVertex),
endEdge = graph.endEdges(currentVertex); endEdge = graph.endEdges(currentVertex); edge != endEdge; ++edge) {
edge != endEdge; ++edge)
{
++degree; ++degree;
} }
if (degree >= maxDegree) {
if( degree >= maxDegree )
{
orderedVertices.emplace_front(currentVertex); orderedVertices.emplace_front(currentVertex);
++firstDegreeChange; ++firstDegreeChange;
if(degree > maxDegree) if (degree > maxDegree) {
{
firstDegreeChange = 1; firstDegreeChange = 1;
maxDegree = degree; maxDegree = degree;
} }
} } else {
else
{
orderedVertices.emplace_back(currentVertex); orderedVertices.emplace_back(currentVertex);
} }
} }
@ -158,10 +141,8 @@ colorVerticesWelshPowell(const Graph& graph)
// order deque by descending degree // order deque by descending degree
std::stable_sort(orderedVertices.begin() + firstDegreeChange, std::stable_sort(orderedVertices.begin() + firstDegreeChange,
orderedVertices.end(), orderedVertices.end(),
[&degrees](const Vertex& v1, const Vertex& v2) [&degrees](const Vertex& v1, const Vertex& v2)
{ { return degrees[v1] > degrees[v2]; });
return degrees[v1] > degrees[v2];
});
// Overwrite degree with color // Overwrite degree with color
auto& colors = degrees; auto& colors = degrees;
@ -171,39 +152,38 @@ colorVerticesWelshPowell(const Graph& graph)
std::vector<std::size_t> verticesPerColor; std::vector<std::size_t> verticesPerColor;
verticesPerColor.reserve(10); verticesPerColor.reserve(10);
while(!orderedVertices.empty()) while (!orderedVertices.empty()) {
{ verticesPerColor.push_back(Detail::colorGraphWelshPowell(graph, orderedVertices,
verticesPerColor colors, color++, noVertices));
.push_back(Detail::colorGraphWelshPowell(graph, orderedVertices, colors,
color++, noVertices));
} }
return std::make_tuple(colors, color, verticesPerColor); return std::make_tuple(colors, color, verticesPerColor);
} }
/// \! Reorder colored graph preserving order of vertices with the same color. /// \! Reorder colored graph preserving order of vertices with the same color.
template<class Graph> template <class Graph>
std::vector<std::size_t> std::vector<std::size_t>
reorderVerticesPreserving(const std::vector<int>& colors, int noColors, reorderVerticesPreserving(const std::vector<int>& colors,
int noColors,
const std::vector<std::size_t>& verticesPerColor, const std::vector<std::size_t>& verticesPerColor,
const Graph& graph) const Graph& graph)
{ {
std::vector<std::size_t> colorIndex(noColors, 0); std::vector<std::size_t> colorIndex(noColors, 0);
std::vector<std::size_t> indices(graph.maxVertex() + 1); std::vector<std::size_t> indices(graph.maxVertex() + 1);
std::partial_sum(verticesPerColor.begin(), std::partial_sum(verticesPerColor.begin(),
verticesPerColor.begin()+verticesPerColor.size() - 1, verticesPerColor.begin() + verticesPerColor.size() - 1,
colorIndex.begin() + 1); colorIndex.begin() + 1);
for(const auto& vertex: graph) for (const auto& vertex : graph) {
{
indices[vertex] = colorIndex[colors[vertex]]++; indices[vertex] = colorIndex[colors[vertex]]++;
} }
return indices; return indices;
} }
/// \! Reorder Vetrices in spheres /// \! Reorder Vetrices in spheres
template<class Graph> template <class Graph>
std::vector<std::size_t> std::vector<std::size_t>
reorderVerticesSpheres(const std::vector<int>& colors, int noColors, reorderVerticesSpheres(const std::vector<int>& colors,
int noColors,
const std::vector<std::size_t>& verticesPerColor, const std::vector<std::size_t>& verticesPerColor,
const Graph& graph, const Graph& graph,
typename Graph::VertexDescriptor root) typename Graph::VertexDescriptor root)
@ -213,26 +193,22 @@ reorderVerticesSpheres(const std::vector<int>& colors, int noColors,
std::vector<std::size_t> indices(graph.maxVertex() + 1, notVisitedTag); std::vector<std::size_t> indices(graph.maxVertex() + 1, notVisitedTag);
using Vertex = typename Graph::VertexDescriptor; using Vertex = typename Graph::VertexDescriptor;
std::partial_sum(verticesPerColor.begin(), std::partial_sum(verticesPerColor.begin(),
verticesPerColor.begin()+verticesPerColor.size() - 1, verticesPerColor.begin() + verticesPerColor.size() - 1,
colorIndex.begin() + 1); colorIndex.begin() + 1);
std::size_t noVisited = 0; std::size_t noVisited = 0;
auto numberer = [&colorIndex, &colors, &indices](Vertex vertex) 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); 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); 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 // Graph is disconnected search for not yet visited node
for(auto vertex: graph) for (auto vertex : graph) {
{ if (indices[vertex] == notVisitedTag) {
if ( indices[vertex] == notVisitedTag )
{
// \todo make sure that this is a peripheral node! // \todo make sure that this is a peripheral node!
root = vertex; root = vertex;
break; break;

314
tests/test_dilu.cpp
View File

@ -31,8 +31,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2NoZeros, T, NumericTypes)
/* /*
Tests that the dilu decomposition mathces the expected result Tests that the dilu decomposition mathces the expected result
for a 2x2 matrix with no zero blocks, with block size 2x2. for a 2x2 matrix with no zero blocks, with block size 2x2.
A A
| | 3 1| | 1 0| | | | 3 1| | 1 0| |
| | 2 1| | 0 1| | | | 2 1| | 0 1| |
| | | |
@ -52,27 +52,27 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2NoZeros, T, NumericTypes)
row.insert(1); row.insert(1);
} }
A[0][0][0][0]=3.0; A[0][0][0][0] = 3.0;
A[0][0][0][1]=1.0; A[0][0][0][1] = 1.0;
A[0][0][1][0]=2.0; A[0][0][1][0] = 2.0;
A[0][0][1][1]=1.0; A[0][0][1][1] = 1.0;
A[0][1][0][0]=1.0; A[0][1][0][0] = 1.0;
A[0][1][1][1]=1.0; A[0][1][1][1] = 1.0;
A[1][0][0][0]=2.0; A[1][0][0][0] = 2.0;
A[1][0][1][1]=2.0; A[1][0][1][1] = 2.0;
A[1][1][0][0]=-1.0; A[1][1][0][0] = -1.0;
A[1][1][1][1]=-1.0; A[1][1][1][1] = -1.0;
auto D_00 = A[0][0]; auto D_00 = A[0][0];
auto D_00_inv = D_00; auto D_00_inv = D_00;
D_00_inv.invert(); D_00_inv.invert();
// D_11 = A_11 - L_10 D_00_inv U_01 // 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<Matrix, Vector, Vector> seqdilu(A); Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
@ -99,8 +99,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2, T, NumericTypes)
/* /*
Tests that the dilu decomposition mathces the expected result Tests that the dilu decomposition mathces the expected result
for a 2x2 matrix, with block size 2x2. for a 2x2 matrix, with block size 2x2.
A A
| | 3 1| | 1 0| | | | 3 1| | 1 0| |
| | 2 1| | 0 1| | | | 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][0] = 3.0;
A[0][0][0][1]=1.0; A[0][0][0][1] = 1.0;
A[0][0][1][0]=2.0; A[0][0][1][0] = 2.0;
A[0][0][1][1]=1.0; A[0][0][1][1] = 1.0;
A[0][1][0][0]=1.0; A[0][1][0][0] = 1.0;
A[0][1][1][1]=1.0; A[0][1][1][1] = 1.0;
A[1][1][0][0]=-1.0; A[1][1][0][0] = -1.0;
A[1][1][1][1]=-1.0; A[1][1][1][1] = -1.0;
auto D_00 = A[0][0]; auto D_00 = A[0][0];
@ -188,20 +188,20 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrectNoZeros, T, NumericTypes)
row.insert(1); row.insert(1);
} }
A[0][0][0][0]=3.0; A[0][0][0][0] = 3.0;
A[0][0][0][1]=1.0; A[0][0][0][1] = 1.0;
A[0][0][1][0]=2.0; A[0][0][1][0] = 2.0;
A[0][0][1][1]=1.0; A[0][0][1][1] = 1.0;
A[0][1][0][0]=1.0; A[0][1][0][0] = 1.0;
A[0][1][1][1]=1.0; A[0][1][1][1] = 1.0;
A[1][0][0][0]=2.0; A[1][0][0][0] = 2.0;
A[1][0][1][1]=2.0; A[1][0][1][1] = 2.0;
A[1][1][0][0]=-1.0; A[1][1][0][0] = -1.0;
A[1][1][1][1]=-1.0; A[1][1][1][1] = -1.0;
Vector x(2); Vector x(2);
@ -221,9 +221,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrectNoZeros, T, NumericTypes)
auto D_00_inv = D_00; auto D_00_inv = D_00;
D_00_inv.invert(); D_00_inv.invert();
// D_11= A_11 - L_10 D_00_inv U_01 // 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; auto D_11_inv = D_11;
D_11_inv.invert(); D_11_inv.invert();
// define: z = M^-1(b - Ax) // define: z = M^-1(b - Ax)
// where: M = (D + L_A) A^-1 (D + U_A) // 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 - D_00_inv*A_01*z_1
z[0] = y[0]; z[0] = y[0];
auto temp = D_00_inv*A[0][1]; auto temp = D_00_inv * A[0][1];
temp.mmv(z[1], z[0]); temp.mmv(z[1], z[0]);
// x_k+1 = x_k + z // x_k+1 = x_k + z
Vector new_x = x; 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][0] = 3.0;
A[0][0][0][1]=1.0; A[0][0][0][1] = 1.0;
A[0][0][1][0]=2.0; A[0][0][1][0] = 2.0;
A[0][0][1][1]=1.0; A[0][0][1][1] = 1.0;
A[0][1][0][0]=1.0; A[0][1][0][0] = 1.0;
A[0][1][1][1]=1.0; A[0][1][1][1] = 1.0;
A[1][1][0][0]=-1.0; A[1][1][0][0] = -1.0;
A[1][1][1][1]=-1.0; A[1][1][1][1] = -1.0;
Vector x(2); 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 // D_11 = A_11 - L_10 D_0_inv U_01 = A_11
auto D_11 = A[1][1]; auto D_11 = A[1][1];
auto D_11_inv = D_11; auto D_11_inv = D_11;
D_11_inv.invert(); D_11_inv.invert();
// define: z = M^-1(b - Ax) // define: z = M^-1(b - Ax)
// where: M = (D + L_A) A^-1 (D + U_A) // 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 - D_00_inv*A_01*z_1
z[0] = y[0]; z[0] = y[0];
auto temp = D_00_inv*A[0][1]; auto temp = D_00_inv * A[0][1];
temp.mmv(z[1], z[0]); temp.mmv(z[1], z[0]);
// x_k+1 = x_k + z // x_k+1 = x_k + z
Vector new_x = x; 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][0] = 3.0;
A[0][0][0][1]=1.0; A[0][0][0][1] = 1.0;
A[0][0][1][0]=2.0; A[0][0][1][0] = 2.0;
A[0][0][1][1]=1.0; A[0][0][1][1] = 1.0;
A[1][1][0][0]=2.0; A[1][1][0][0] = 2.0;
A[1][1][1][1]=2.0; A[1][1][1][1] = 2.0;
A[1][1][0][0]=-1.0; A[1][1][0][0] = -1.0;
A[1][1][1][1]=-1.0; A[1][1][1][1] = -1.0;
Vector x(2); Vector x(2);
x[0][0] = 1.0; 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 // D_11 = A_11 - L_10 D_0_inv U_01 = A_11
auto D_11 = A[1][1]; auto D_11 = A[1][1];
auto D_11_inv = D_11; auto D_11_inv = D_11;
D_11_inv.invert(); D_11_inv.invert();
// define: z = M^-1(b - Ax) // define: z = M^-1(b - Ax)
// where: M = (D + L_A) A^-1 (D + U_A) // 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 Tests that the dilu decomposition mathces the expected result
for a 3x3 matrix, with block size 3x3. for a 3x3 matrix, with block size 3x3.
A A
| | 3 1 2| | 0 0 0| | 0 0 0| | | | 3 1 2| | 0 0 0| | 0 0 0| |
| | 2 3 1| | 0 0 0| | 0 0 0| | | | 2 3 1| | 0 0 0| | 0 0 0| |
| | 2 1 0| | 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| | 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| | 3 1 3| | 0 1 3| |
| | | |
| | 0 0 0| | 1 0 2| | 1 3 2| | | | 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) { for (auto row = A.createbegin(); row != A.createend(); ++row) {
if (row.index() == 0) { if (row.index() == 0) {
row.insert(row.index()); row.insert(row.index());
} } else if (row.index() == 1) {
else if (row.index() == 1) {
row.insert(row.index()); row.insert(row.index());
row.insert(row.index() + 1); row.insert(row.index() + 1);
} } else if (row.index() == 2) {
else if (row.index() == 2) {
row.insert(row.index() - 1); row.insert(row.index() - 1);
row.insert(row.index()); 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][0] = 3.0;
A[0][0][0][1]=1.0; A[1][1][0][1]=0.0; A[1][2][0][1]=0.0; A[1][1][0][0] = 1.0;
A[0][0][0][2]=2.0; A[1][1][0][2]=1.0; A[1][2][0][2]=2.0; A[1][2][0][0] = 1.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][0][1] = 1.0;
A[0][0][1][1]=3.0; A[1][1][1][1]=1.0; A[1][2][1][1]=1.0; A[1][1][0][1] = 0.0;
A[0][0][1][2]=1.0; A[1][1][1][2]=0.0; A[1][2][1][2]=1.0; A[1][2][0][1] = 0.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][0][2] = 2.0;
A[0][0][2][1]=1.0; A[1][1][2][1]=1.0; A[1][2][2][1]=1.0; A[1][1][0][2] = 1.0;
A[0][0][2][2]=0.0; A[1][1][2][2]=3.0; A[1][2][2][2]=3.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][0] = 1.0;
A[2][1][0][1]=0.0; A[2][2][0][1]=3.0; A[2][2][0][0] = 1.0;
A[2][1][0][2]=2.0; A[2][2][0][2]=2.0; A[2][1][0][1] = 0.0;
A[2][1][1][0]=0.0; A[2][2][1][0]=2.0; A[2][2][0][1] = 3.0;
A[2][1][1][1]=1.0; A[2][2][1][1]=1.0; A[2][1][0][2] = 2.0;
A[2][1][1][2]=4.0; A[2][2][1][2]=3.0; A[2][2][0][2] = 2.0;
A[2][1][2][0]=5.0; A[2][2][2][0]=3.0; A[2][1][1][0] = 0.0;
A[2][1][2][1]=1.0; A[2][2][2][1]=1.0; A[2][2][1][0] = 2.0;
A[2][1][2][2]=1.0; A[2][2][2][2]=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]; 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 // D_11 = A_11 - L_10 D_00_inv U_01 = A_11
auto D_11 = A[1][1]; auto D_11 = A[1][1];
auto D_11_inv = D_11; 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 // 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; auto D_22_inv = D_22;
D_22_inv.invert(); D_22_inv.invert();
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A); Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
auto Dinv = seqdilu.getDiagonal(); 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| | | | 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| | 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| | 3 1 3| | 0 1 3| | | |2| | | |2| |
| | | | | | | | | | | |
| | 0 0 0| | 1 0 2| | 1 3 2| | | |1| | | |0| | | | 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| | 0 1 4| | 2 1 3| | | |0| | | |2| |
| | 0 0 0| | 5 1 1| | 3 1 2| | | |2| | | |1| | | | 0 0 0| | 5 1 1| | 3 1 2| | | |2| | | |1| |
*/ */
const int N = 3; 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) { for (auto row = A.createbegin(); row != A.createend(); ++row) {
if (row.index() == 0) { if (row.index() == 0) {
row.insert(row.index()); row.insert(row.index());
} } else if (row.index() == 1) {
else if (row.index() == 1) {
row.insert(row.index()); row.insert(row.index());
row.insert(row.index() + 1); row.insert(row.index() + 1);
} } else if (row.index() == 2) {
else if (row.index() == 2) {
row.insert(row.index() - 1); row.insert(row.index() - 1);
row.insert(row.index()); 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][0] = 3.0;
A[0][0][0][1]=1.0; A[1][1][0][1]=0.0; A[1][2][0][1]=0.0; A[1][1][0][0] = 1.0;
A[0][0][0][2]=2.0; A[1][1][0][2]=1.0; A[1][2][0][2]=2.0; A[1][2][0][0] = 1.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][0][1] = 1.0;
A[0][0][1][1]=3.0; A[1][1][1][1]=1.0; A[1][2][1][1]=1.0; A[1][1][0][1] = 0.0;
A[0][0][1][2]=1.0; A[1][1][1][2]=0.0; A[1][2][1][2]=1.0; A[1][2][0][1] = 0.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][0][2] = 2.0;
A[0][0][2][1]=1.0; A[1][1][2][1]=1.0; A[1][2][2][1]=1.0; A[1][1][0][2] = 1.0;
A[0][0][2][2]=0.0; A[1][1][2][2]=3.0; A[1][2][2][2]=3.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][0] = 1.0;
A[2][1][0][1]=0.0; A[2][2][0][1]=3.0; A[2][2][0][0] = 1.0;
A[2][1][0][2]=2.0; A[2][2][0][2]=2.0; A[2][1][0][1] = 0.0;
A[2][1][1][0]=0.0; A[2][2][1][0]=2.0; A[2][2][0][1] = 3.0;
A[2][1][1][1]=1.0; A[2][2][1][1]=1.0; A[2][1][0][2] = 2.0;
A[2][1][1][2]=4.0; A[2][2][1][2]=3.0; A[2][2][0][2] = 2.0;
A[2][1][2][0]=5.0; A[2][2][2][0]=3.0; A[2][1][1][0] = 0.0;
A[2][1][2][1]=1.0; A[2][2][2][1]=1.0; A[2][2][1][0] = 2.0;
A[2][1][2][2]=1.0; A[2][2][2][2]=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); Vector x(3);
x[0][0] = 1.0; x[1][0] = 1.0; x[2][0] = 1.0; x[0][0] = 1.0;
x[0][1] = 2.0; x[1][1] = 3.0; x[2][1] = 0.0; x[1][0] = 1.0;
x[0][2] = 3.0; x[1][2] = 2.0; x[2][2] = 2.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); Vector b(3);
b[0][0] = 2.0; b[1][0] = 2.0; b[2][0] = 0.0; b[0][0] = 2.0;
b[0][1] = 1.0; b[1][1] = 3.0; b[2][1] = 2.0; b[1][0] = 2.0;
b[0][2] = 2.0; b[1][2] = 2.0; b[2][2] = 1.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 // D_00 = A_00
@ -652,17 +714,17 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect3, T, NumericTypes)
// = A_11 // = A_11
auto D_11 = A[1][1]; auto D_11 = A[1][1];
auto D_11_inv = D_11; 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 // 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 // = 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; auto D_22_inv = D_22;
D_22_inv.invert(); D_22_inv.invert();
// define: z = M^-1(b - Ax) // define: z = M^-1(b - Ax)
// where: M = (D + L_A) A^-1 (D + U_A) // where: M = (D + L_A) A^-1 (D + U_A)
// lower triangular solve: (E + L) y = b - Ax // lower triangular solve: (E + L) y = b - Ax
Vector y(3); Vector y(3);
// y_0 = D_00_inv*[b_0 - (A_00*x_0 + A_01*x_1)] // y_0 = D_00_inv*[b_0 - (A_00*x_0 + A_01*x_1)]
// = D_00_inv*[b_0 - A_00*x_0] // = 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 - D_11_inv*A_12*z_2
z[1] = y[1]; z[1] = y[1];
auto temp = D_11_inv*A[1][2]; auto temp = D_11_inv * A[1][2];
temp.mmv(z[2], z[1]); 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 - D_00_inv(A_01*z_1 + A_02*z_2)
// z_0 = y_0 // z_0 = y_0
@ -741,13 +803,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsEqualToDuneSeqILUApply, T, NumericTy
row.insert(row.index()); row.insert(row.index());
} }
A[0][0][0][0]=3.0; A[0][0][0][0] = 3.0;
A[0][0][0][1]=1.0; A[0][0][0][1] = 1.0;
A[0][0][1][0]=2.0; A[0][0][1][0] = 2.0;
A[0][0][1][1]=1.0; A[0][0][1][1] = 1.0;
A[1][1][0][0]=-1.0; A[1][1][0][0] = -1.0;
A[1][1][1][1]=-1.0; A[1][1][1][1] = -1.0;
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A); Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);