diff --git a/dune/matrix-vector/triangularsolve.hh b/dune/matrix-vector/triangularsolve.hh index fa014d17bf84aa34f652a2469b3c0012c9aa7752..76bd1334cd378ed2dbdd73efcd4f727bff245327 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();