corrected bugs introduced in yesterweeks changes

[[Imported from SVN: r12159]]

dune/solvers/norms/diagnorm.hh

@@ -15,43 +15,53 @@
   *  \tparam Diagonal some container type having the random access operator[]; used to store the diagonal entries of the matrix
   *  \tparam VectorType the vector type the norm may be applied to
   */
-template <class Diagonal=Dune::BlockVector<Dune::FieldVector <double,1> >, class VectorType=Dune::BlockVector<Dune::FieldVector <double,1> > >
+template <class Diagonal=Dune::BlockVector<Dune::FieldVector <double,1> >, class V=Dune::BlockVector<Dune::FieldVector <double,1> > >
 class DiagNorm:
-    public Norm<VectorType>
+    public Norm<V>
 {
+    public:
+        typedef V VectorType;
+        typedef typename VectorType::field_type field_type;
+    private:
         typedef typename VectorType::size_type SizeType;
 
     public:
-        DiagNorm(const double alpha, const Diagonal &d) :
+        DiagNorm(const field_type alpha, const Diagonal &d) :
             d(d),
             alpha(alpha)
         {}
 
         //! Compute the norm of the given vector
-        double operator()(const VectorType &v) const
+        field_type operator()(const VectorType &f) const
+        {
+            return std::sqrt(this->DiagNorm<Diagonal,VectorType>::normSquared(f));
+        }
+
+        //! Compute the square of the norm of the given vector
+        virtual field_type normSquared(const VectorType& f) const
         {
-            double r = 0.0;
+            field_type r = 0.0;
 
             typename VectorType::value_type Dv(0.0);
-            for(SizeType row = 0; row < v.size(); ++row)
+            for(SizeType row = 0; row < f.size(); ++row)
             {
-                d[row].mv(v[row],Dv);
-                r += Dv*v[row];
+                d[row].mv(f[row],Dv);
+                r += Dv*f[row];
             }
 
-            return std::sqrt(std::abs(alpha * r));
+            return std::abs(alpha * r);
         }
 
         //! Compute the norm of the difference of two vectors
-        double diff(const VectorType &v1, const VectorType &v2) const
+        field_type diff(const VectorType &f1, const VectorType &f2) const
         {
-            double r = 0.0;
+            field_type r = 0.0;
 
             typename VectorType::value_type Dv(0.0);
-            for(SizeType row = 0; row < v1.size(); ++row)
+            for(SizeType row = 0; row < f1.size(); ++row)
             {
-                d[row].mv(v1[row]-v2[row],Dv);
-                r += Dv*(v1[row]-v2[row]);
+                d[row].mv(f1[row]-f2[row],Dv);
+                r += Dv*(f1[row]-f2[row]);
             }
 
             return std::sqrt(std::abs(alpha * r));
@@ -60,7 +70,7 @@ class DiagNorm:
     private:
         const Diagonal& d;
 
-        const double alpha;
+        const field_type alpha;
 
 };
 
@@ -68,33 +78,42 @@ template<>
 class DiagNorm<Dune::BlockVector<Dune::FieldVector <double,1> >, Dune::BlockVector<Dune::FieldVector <double,1> > >:
     public Norm<Dune::BlockVector<Dune::FieldVector <double,1> > >
 {
-        typedef Dune::BlockVector<Dune::FieldVector <double,1> > VectorType;
+    public:
+        typedef double field_type;
+        typedef Dune::BlockVector<Dune::FieldVector <field_type,1> > VectorType;
+    private:
         typedef VectorType::size_type SizeType;
 
     public:
-        DiagNorm(const double alpha, const VectorType &d) :
+        DiagNorm(const field_type alpha, const VectorType &d) :
             d(d),
             alpha(alpha)
         {}
 
         //! Compute the norm of the given vector
-        double operator()(const VectorType &v) const
+        field_type operator()(const VectorType &f) const
+        {
+            return std::sqrt(this->DiagNorm<>::normSquared(f));
+        }
+
+        //! Compute the square of the norm of the given vector
+        virtual field_type normSquared(const VectorType& f) const
         {
-            double r = 0.0;
+            field_type r = 0.0;
 
-            for(SizeType row = 0; row < v.size(); ++row)
-                r += d[row] * v[row] * v[row];
+            for(SizeType row = 0; row < f.size(); ++row)
+                r += d[row] * f[row] * f[row];
 
-            return std::sqrt(std::abs(alpha * r));
+            return std::abs(alpha * r);
         }
 
         //! Compute the norm of the difference of two vectors
-        double diff(const VectorType &v1, const VectorType &v2) const
+        field_type diff(const VectorType &f1, const VectorType &f2) const
         {
-            double r = 0.0;
+            field_type r = 0.0;
 
-            for (SizeType row = 0; row < v1.size(); ++row)
-                r += (double)d[row] * (v1[row]-v2[row]) * (v1[row] - v2[row]);
+            for (SizeType row = 0; row < f1.size(); ++row)
+                r += (field_type)d[row] * (f1[row]-f2[row]) * (f1[row] - f2[row]);
 
             return std::sqrt(std::abs(alpha * r));
         }
@@ -102,7 +121,7 @@ class DiagNorm<Dune::BlockVector<Dune::FieldVector <double,1> >, Dune::BlockVect
     private:
         const VectorType &d;
 
-        const double alpha;
+        const field_type alpha;
 
 };