From 6e9d5b166ac040df4115061a1b6270746a2c46bb Mon Sep 17 00:00:00 2001
From: Max Kahnt <max.kahnt@fu-berlin.de>
Date: Thu, 22 May 2014 14:53:12 +0200
Subject: [PATCH] Fix TruncatedBlockGSStep localExactSolve.
Implementation mistakenly changed the solved system for indefinite
matrices with diagonal entries of value 0.
The idea is to also be able to tackle systems that are semi-definite
due to trivial rows by simply not altering the initial x in these components.
This was accomplished by setting the corresponding diagonal entry to 1.
Updated the check for trivial rows such that it is not triggered
if only the diagonal element is zero but only the entire row.
---
.../iterationsteps/truncatedblockgsstep.hh | 27 +++++++++++++------
1 file changed, 19 insertions(+), 8 deletions(-)
diff --git a/dune/solvers/iterationsteps/truncatedblockgsstep.hh b/dune/solvers/iterationsteps/truncatedblockgsstep.hh
index 6069091..6d93b05 100644
--- a/dune/solvers/iterationsteps/truncatedblockgsstep.hh
+++ b/dune/solvers/iterationsteps/truncatedblockgsstep.hh
@@ -293,19 +293,30 @@ public:
{
for(typename MBlock::size_type i=0; i<A.N(); ++i)
{
+ assert(A.exists(i,i)); // No implementation for matrices w/o explicit diagonal elements is available here.
+
+ bool rowIsZero = true;
if (ignore[i])
{
A[i] = 0.0;
}
+ else
+ {
+ typename MBlock::row_type::Iterator inner_it = A[i].begin();
+ typename MBlock::row_type::Iterator inner_end = A[i].end();
+ for(; inner_it!=inner_end; ++inner_it)
+ if (*inner_it != 0.0)
+ {
+ rowIsZero = false;
+ break;
+ }
+ }
- typename MBlock::row_type::Iterator inner_it = A[i].begin();
- typename MBlock::row_type::Iterator inner_end = A[i].end();
- for(; inner_it!=inner_end; ++inner_it)
- if (inner_it.index()==i and *inner_it==0.0)
- {
- *inner_it = 1.0;
- b[i] = x[i];
- }
+ if (rowIsZero)
+ {
+ A[i][i] = 1.0;
+ b[i] = x[i];
+ }
}
A.solve(x,b);
--
GitLab