diff --git a/dune/matrix-vector/triangularsolve.hh b/dune/matrix-vector/triangularsolve.hh index 61cffac213b1dfff6dab740eb7e37d952ce27e0b..7bf57e4a39addf76b7607ce37f8774ef7bb03b1d 100644 --- a/dune/matrix-vector/triangularsolve.hh +++ b/dune/matrix-vector/triangularsolve.hh @@ -19,27 +19,29 @@ namespace MatrixVector { if (ignore != nullptr and (*ignore)[i].all()) continue; auto cIt = it->begin(); - assert(cIt.index() == it.index()); + assert(cIt.index() == i); x[i] = b[i] / *cIt; cIt++; for (; cIt != it->end(); ++cIt) { const size_t j = cIt.index(); + assert(j > i); b[j] -= x[i] * *cIt; } } } else { for (auto it = L.begin(); it != L.end(); ++it) { - size_t i = it.index(); + const size_t i = it.index(); if (ignore != nullptr and (*ignore)[i].all()) continue; for (auto cIt = it->begin(); cIt != it->end(); ++cIt) { - size_t j = cIt.index(); - if (i == j) { - x[i] = b[i] / *cIt; - break; + const size_t j = cIt.index(); + assert(j <= i); + if (j < i) { + b[i] -= *cIt * x[j]; + continue; } - assert(j < i); - b[i] -= *cIt * x[j]; + assert(j == i); + x[i] = b[i] / *cIt; } } } @@ -56,36 +58,32 @@ namespace MatrixVector { x = 0; if (transpose) { 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; + const size_t i = it.index(); if (ignore != nullptr and (*ignore)[i].all()) continue; - x[i] = b[i] / diagonal; - + auto cIt = it->beforeEnd(); + assert(cIt.index() == i); + x[i] = b[i] / *cIt; cIt--; for (; cIt != it->beforeBegin(); --cIt) { - size_t j = cIt.index(); + const size_t j = cIt.index(); assert(j < i); b[j] -= *cIt * x[i]; } } } else { for (auto it = U.beforeEnd(); it != U.beforeBegin(); --it) { - size_t i = it.index(); + const size_t i = it.index(); if (ignore != nullptr and (*ignore)[i].all()) continue; - auto cIt = it->begin(); - assert(cIt.index() == i); - auto diagonal = *cIt; - cIt++; - for (; cIt != it->end(); ++cIt) { - size_t j = cIt.index(); + auto cIt = it->beforeEnd(); + for (; cIt != it->begin(); --cIt) { + const size_t j = cIt.index(); assert(j > i); b[i] -= *cIt * x[j]; } - x[i] = b[i] / diagonal; + assert(cIt.index() == i); + x[i] = b[i] / *cIt; } } }