From 1799e7f2d9870044cae8762b558e37e1e132d56e Mon Sep 17 00:00:00 2001
From: Uli Sack <usack@math.fu-berlin.de>
Date: Fri, 20 Apr 2012 13:52:30 +0000
Subject: [PATCH] if for a semi definite Matrix normSquared(f) is negative zero
i.e. e.g. -1e-16 operator() will return -nan. Now normSquared returns 0 in
such a case and throws an exception for 'really' negative values
[[Imported from SVN: r6070]]
---
dune/solvers/norms/energynorm.hh | 35 +++++++++++++++++++++++++++++---
1 file changed, 32 insertions(+), 3 deletions(-)
diff --git a/dune/solvers/norms/energynorm.hh b/dune/solvers/norms/energynorm.hh
index 7c8f1521..280e95f8 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:
--
GitLab