Fix docu. Add convenience methods.

dune/matrix-vector/triangularsolve.hh

@@ -6,8 +6,19 @@ 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,
@@ -23,6 +34,8 @@ namespace MatrixVector {
         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;
         }
@@ -46,9 +59,21 @@ namespace MatrixVector {
     }
   }
 
+  /**
+   * \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 = 0;
+    lowerTriangularSolve(L, std::move(b), x, ignore, transpose);
+    return x;
+  }
+
   /**
    * \brief Solves U*x=b where U is an upper triangular matrix.
-   * \param transpose Set true if U is passed as a lower triangular matrix.
+   *        See @lowerTriangularSolve for details.
    */
   template <class Matrix, class Vector, class BitVector>
   static void upperTriangularSolve(Matrix const& U, Vector b, Vector& x,
@@ -65,6 +90,8 @@ namespace MatrixVector {
         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];
         }
@@ -85,7 +112,20 @@ namespace MatrixVector {
       }
     }
   }
-}
-}
 
-#endif
+  /**
+   * \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 = 0;
+    upperTriangularSolve(U, std::move(b), x, ignore, transpose);
+    return x;
+  }
+
+} // end namespace MatrixVector
+} // end namespace Dune
+
+#endif // DUNE_MATRIX_VECTOR_TRIANGULARSOLVE_HH