diff --git a/dune/solvers/common/staticmatrixtools.hh b/dune/solvers/common/staticmatrixtools.hh index f95c5791d57593c6d4797c77b4c5440435dd2ac3..70bd60caa9eee6b690ebbe1baa25545cc364672f 100644 --- a/dune/solvers/common/staticmatrixtools.hh +++ b/dune/solvers/common/staticmatrixtools.hh @@ -3,8 +3,6 @@ #ifndef STATIC_MATRIX_TOOL_HH #define STATIC_MATRIX_TOOL_HH -#include <cassert> - #include "dune/common/diagonalmatrix.hh" #include "dune/common/fmatrix.hh" #include "dune/istl/scaledidmatrix.hh" @@ -460,87 +458,7 @@ class StaticMatrix } } - template <class Matrix, class Vector, class BitVector> - static void lowerTriangularSolve(Matrix const& A, Vector const& b, - Vector& x, BitVector const* ignore, - bool transpose = false) { - static_assert( - Matrix::block_type::rows == 1 and Matrix::block_type::cols == 1, - "Only implemented for scalar problems"); - x = 0; - Vector r = b; - if (transpose) { - for (auto it = A.begin(); it != A.end(); ++it) { - const size_t i = it.index(); - if (ignore != nullptr and (*ignore)[i].all()) - continue; - auto cIt = it->begin(); - assert(cIt.index() == it.index()); - x[i] = r[i] / *cIt; - for (; cIt != it->end(); ++cIt) { - const size_t j = cIt.index(); - r[j] -= x[i] * *cIt; - } - } - } else { - for (auto it = A.begin(); it != A.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(); - if (i == j) { - x[i] = r[i] / *cIt; - break; - } - assert(j < i); - r[i] -= *cIt * x[j]; - } - } - } - } - - template <class Matrix, class Vector, class BitVector> - static void upperTriangularSolve(Matrix const& U, Vector const& b, - Vector& x, BitVector const* ignore, - bool transpose = false) { - static_assert( - Matrix::block_type::rows == 1 and Matrix::block_type::cols == 1, - "Only implemented for scalar problems"); - x = 0; - Vector r = b; - if (transpose) { - 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(); - assert(cIt.index() == i); - x[i] = r[i] / *cIt; - cIt--; - for (; cIt != it->beforeBegin(); --cIt) { - const size_t j = cIt.index(); - assert(j < i); - r[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->begin(); - assert(cIt.index() == i); - auto diagonal = *cIt; - cIt++; - for (; cIt != it->end(); ++cIt) { - const size_t j = cIt.index(); - r[i] -= *cIt * x[j]; - } - x[i] = r[i] / diagonal; - } - } - } }; #endif +