Add asserts. Use consistent iteration scheme.

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;
       }
     }
   }