diff --git a/dune/solvers/norms/diagnorm.hh b/dune/solvers/norms/diagnorm.hh index b3c941118e09b2d79f0c29c9937e53582f404852..e438ce68944ac1cf1f5f0bc3742432cf94ab2966 100644 --- a/dune/solvers/norms/diagnorm.hh +++ b/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; };