diff --git a/dune/solvers/norms/energynorm.hh b/dune/solvers/norms/energynorm.hh index 7c8f1521d3d1829f589756e7028d940bd1fa5ed0..280e95f86f9605a9d86e43e23791ddd021b6c951 100644 --- a/dune/solvers/norms/energynorm.hh +++ b/dune/solvers/norms/energynorm.hh @@ -62,7 +62,8 @@ return sqrt(this->EnergyNorm<OperatorType,DiscFuncType>::normSquared(f)); } - //! Compute the square of the norm of the given vector + /** \brief Compute the square of the norm of the given vector + \todo This could be implemented without the temporary. */ virtual double normSquared(const DiscFuncType& f) const { if (iterationStep_ == NULL && matrix_ == NULL) @@ -76,7 +77,20 @@ tmp = 0; A.umv(f, tmp); - return f*tmp; + double ret = f*tmp; + + if (ret < 0) + { + if (ret < -1e-13) + { + char msg[1024]; + sprintf(msg, "Supplied linear operator is not positive (semi-)definite: (u,Au) = %e", ret); + DUNE_THROW(Dune::RangeError, msg); + } + ret = 0.0; + } + + return ret; } /** \brief Compute the squared norm for a given vector and matrix @@ -87,7 +101,22 @@ DiscFuncType tmp(u.size()); tmp = 0; A.umv(u, tmp); - return u*tmp; + + double ret = u*tmp; + + if (ret < 0) + { + if (ret < -1e-13) + { + char msg[1024]; + sprintf(msg, "Supplied linear operator is not positive (semi-)definite: (u,Au) = %e", ret); + DUNE_THROW(Dune::RangeError, msg); + } + ret = 0.0; + } + + return ret; + } DUNE_DEPRECATED; protected: