Merge branch 'feature/triangular-solve'

dune/matrix-vector/CMakeLists.txt

@@ -16,4 +16,5 @@ install(FILES
   singlenonzerorowmatrix.hh
   tranpose.hh
   transformmatrix.hh
+  triangularsolve.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/matrix-vector)
\ No newline at end of file

dune/matrix-vector/triangularsolve.hh

+#ifndef DUNE_MATRIX_VECTOR_TRIANGULARSOLVE_HH
+#define DUNE_MATRIX_VECTOR_TRIANGULARSOLVE_HH
+
+namespace Dune {
+namespace MatrixVector {
+
+  /**
+   * \brief Solves L*x=b where L is a lower triangular matrix.
+   * \param x Solution vector. Make sure the ignored nodes are set
+   *        correctly because their respective values for x as passed to
+   *        this method may affect other lines.
+   * \param ignore Tags the lines of the system where x shall not be
+   *        computed. Note that ignored nodes might still influence the
+   *        solution of other lines if not set to zero.
+   * \param transpose Set true if L is passed as an upper triangular matrix.
+   */
+  // Note: The current assert makes sure we actually deal with triangular
+  //       matrices. It might be favorable to support matrices as well
+  //       where other entries are not required to be non-existent, e.g.
+  //       if one has a in-place decomposition and does not want to
+  //       reimplement the iterators.
+  template <class Matrix, class Vector, class BitVector>
+  static void lowerTriangularSolve(Matrix const& L, Vector b, Vector& x,
+                                   BitVector const* ignore,
+                                   bool transpose = false) {
+    if (transpose) {
+      for (auto it = L.begin(); it != L.end(); ++it) {
+        const size_t i = it.index();
+        auto cIt = it->begin();
+        assert(cIt.index() == i);
+        if (ignore == nullptr or (*ignore)[cIt.index()].none())
+          x[i] = b[i] / *cIt;
+        cIt++;
+        for (; cIt != it->end(); ++cIt) {
+          const size_t j = cIt.index();
+          assert(j > i);
+          // Note: We could drop the check for ignore nodes here bcs. b[j] will
+          // be ignored anyway due to the check above.
+          if (ignore == nullptr or (*ignore)[j].none())
+            b[j] -= x[i] * *cIt;
+        }
+      }
+    } else {
+      for (auto it = L.begin(); it != L.end(); ++it) {
+        const size_t i = it.index();
+        if (ignore != nullptr and (*ignore)[i].all())
+          continue;
+        for (auto cIt = it->begin(); cIt != it->end(); ++cIt) {
+          const size_t j = cIt.index();
+          assert(j <= i);
+          if (j < i) {
+            b[i] -= *cIt * x[j];
+            continue;
+          }
+          assert(j == i);
+          x[i] = b[i] / *cIt;
+        }
+      }
+    }
+  }
+
+  /**
+   * \brief Same as lowerTriangularSolve. Ignored nodes are set to 0.
+   */
+  template <class Matrix, class Vector, class BitVector>
+  static Vector lowerTriangularSolve(Matrix const& L, Vector b,
+                                     BitVector const* ignore,
+                                     bool transpose = false) {
+    Vector x = b;
+    x = 0;
+    lowerTriangularSolve(L, std::move(b), x, ignore, transpose);
+    return x;
+  }
+
+  /**
+   * \brief Solves U*x=b where U is an upper triangular matrix.
+   *        See @lowerTriangularSolve for details.
+   */
+  template <class Matrix, class Vector, class BitVector>
+  static void upperTriangularSolve(Matrix const& U, Vector b, Vector& x,
+                                   BitVector const* ignore,
+                                   bool transpose = false) {
+    if (transpose) {
+      for (auto it = U.beforeEnd(); it != U.beforeBegin(); --it) {
+        const size_t i = it.index();
+        auto cIt = it->beforeEnd();
+        assert(cIt.index() == i);
+        if (ignore == nullptr or (*ignore)[cIt.index()].none())
+          x[i] = b[i] / *cIt;
+        cIt--;
+        for (; cIt != it->beforeBegin(); --cIt) {
+          const size_t j = cIt.index();
+          assert(j < i);
+          // Note: We could drop the check for ignore nodes here bcs. b[j] will
+          // be ignored anyway due to the check above.
+          if (ignore == nullptr or (*ignore)[j].none())
+            b[j] -= *cIt * x[i];
+        }
+      }
+    } else {
+      for (auto it = U.beforeEnd(); it != U.beforeBegin(); --it) {
+        const size_t i = it.index();
+        if (ignore != nullptr and (*ignore)[i].all())
+          continue;
+        auto cIt = it->beforeEnd();
+        for (; cIt != it->begin(); --cIt) {
+          const size_t j = cIt.index();
+          assert(j > i);
+          b[i] -= *cIt * x[j];
+        }
+        assert(cIt.index() == i);
+        x[i] = b[i] / *cIt;
+      }
+    }
+  }
+
+  /**
+   * \brief Same as upperTriangularSolve. Ignored nodes are set to 0.
+   */
+  template <class Matrix, class Vector, class BitVector>
+  static Vector upperTriangularSolve(Matrix const& U, Vector b,
+                                     BitVector const* ignore,
+                                     bool transpose = false) {
+    Vector x = b;
+    x = 0;
+    upperTriangularSolve(U, std::move(b), x, ignore, transpose);
+    return x;
+  }
+
+} // end namespace MatrixVector
+} // end namespace Dune
+
+#endif // DUNE_MATRIX_VECTOR_TRIANGULARSOLVE_HH