diff --git a/dune/solvers/norms/fullnorm.hh b/dune/solvers/norms/fullnorm.hh index 51d5dd1c9d0945312b0d771077b508ec209a458f..6a34d891cb1c1c57b05310190387d9699669d942 100644 --- a/dune/solvers/norms/fullnorm.hh +++ b/dune/solvers/norms/fullnorm.hh @@ -22,31 +22,39 @@ template <class LowRankFactor=Dune::BlockVector<Dune::FieldVector <double,1> >, class FullNorm: public Norm<VectorType> { public: - FullNorm(const double alpha, const LowRankFactor &lowRankFactor) : + typedef typename VectorType::field_type field_type; + + FullNorm(const field_type alpha, const LowRankFactor &lowRankFactor) : lowRankFactor_(lowRankFactor), alpha(alpha) {} - //! Compute the norm of the given vector - double operator()(const VectorType &v) const + //! Compute the norm of the given vector + field_type operator()(const VectorType &f) const + { + return std::sqrt(this->FullNorm<LowRankFactor,VectorType>::normSquared(f)); + } + + //! Compute the square of the norm of the given vector + virtual field_type normSquared(const VectorType& f) const { VectorType r(lowRankFactor_.N()); r = 0.0; - lowRankFactor_.umv(v,r); + lowRankFactor_.umv(f,r); - return std::sqrt(std::abs(alpha*(r*r))); + return std::abs(alpha*(r*r)); } - //! Compute the norm of the difference of two vectors - double diff(const VectorType &v1, const VectorType &v2) const + //! Compute the norm of the difference of two vectors + field_type diff(const VectorType &f1, const VectorType &f2) const { VectorType r(lowRankFactor_.N()); r = 0.0; for (typename LowRankFactor::size_type k=0; k<lowRankFactor_.N(); ++k) for (typename LowRankFactor::size_type i=0; i<lowRankFactor_.M(); ++i) - lowRankFactor_[k][i].umv(v1[i]-v2[i],r[k]); + lowRankFactor_[k][i].umv(f1[i]-f2[i],r[k]); return std::sqrt(std::abs(alpha*(r*r))); } @@ -54,7 +62,7 @@ class FullNorm: public Norm<VectorType> private: const LowRankFactor& lowRankFactor_; - const double alpha; + const field_type alpha; }; @@ -65,30 +73,38 @@ class FullNorm<Dune::BlockVector<Dune::FieldVector<double,1> >, Dune::BlockVecto typedef Dune::BlockVector<Dune::FieldVector<double,1> > VectorType; typedef VectorType::size_type SizeType; + typedef double field_type; + public: - FullNorm(const double alpha, const VectorType &m) : + FullNorm(const field_type alpha, const VectorType &m) : m(m), 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->FullNorm<>::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 += m[row] * v[row]; + for (SizeType row = 0; row < f.size(); ++row) + r += m[row] * f[row]; - return std::sqrt(std::abs(alpha*r*r)); + return std::abs(alpha*r*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)m[row] * (v1[row] - v2[row]); + for (SizeType row = 0; row < f1.size(); ++row) + r += (field_type)m[row] * (f1[row] - f2[row]); return std::sqrt(std::abs(alpha*r*r)); } @@ -96,7 +112,7 @@ class FullNorm<Dune::BlockVector<Dune::FieldVector<double,1> >, Dune::BlockVecto private: const VectorType &m; - const double alpha; + const field_type alpha; }; #endif