diff --git a/dune/matrix-vector/triangularsolve.hh b/dune/matrix-vector/triangularsolve.hh index 9289992e7bd0f4d6323d2aa8ad31653adc569ba7..e41f362f71a9c41aa9f18059877ef4bc677e8da3 100644 --- a/dune/matrix-vector/triangularsolve.hh +++ b/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