if for a semi definite Matrix normSquared(f) is negative zero i.e. e.g. -1e-16...

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

@@ -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: