From de1a2cd9af20f438d356e63e409c49ed8b8e3b65 Mon Sep 17 00:00:00 2001
From: Max Kahnt <max.kahnt@fu-berlin.de>
Date: Tue, 26 Jul 2016 11:53:27 +0200
Subject: [PATCH] Add upper solver w/o additional copy
---
dune/matrix-vector/triangularsolve.hh | 17 ++++++++++++++++-
1 file changed, 16 insertions(+), 1 deletion(-)
diff --git a/dune/matrix-vector/triangularsolve.hh b/dune/matrix-vector/triangularsolve.hh
index fa014d1..76bd133 100644
--- a/dune/matrix-vector/triangularsolve.hh
+++ b/dune/matrix-vector/triangularsolve.hh
@@ -36,7 +36,22 @@ namespace MatrixVector {
x = 0;
Vector r = b;
if (transpose) {
- // TODO: not yet handled
+ for (auto it = U.beforeEnd(); it != U.beforeBegin(); --it) {
+ size_t i = it.index();
+ auto cIt = it->beforeEnd();
+ assert(cIt.index() == i);
+ auto diagonal = *cIt;
+ if (ignore != nullptr and (*ignore)[i].all())
+ continue;
+ x[i] = r[i] / diagonal;
+
+ cIt--;
+ for (; cIt != it->beforeBegin(); --cIt) {
+ size_t j = cIt.index();
+ assert(j < i);
+ r[j] -= *cIt * x[i];
+ }
+ }
} else {
for (auto it = U.beforeEnd(); it != U.beforeBegin(); --it) {
size_t i = it.index();
--
GitLab