diff --git a/dune/solvers/norms/fullnorm.hh b/dune/solvers/norms/fullnorm.hh index d8b8fb86df737918e6541b890f40e778db17c634..4ef80d4077db8bc9ecb8d8872507d9fe0ceac0ff 100644 --- a/dune/solvers/norms/fullnorm.hh +++ b/dune/solvers/norms/fullnorm.hh @@ -1,15 +1,70 @@ -#ifndef __FULLNORM__HH__ -#define __FULLNORM__HH__ +#ifndef FULLNORM_HH +#define FULLNORM_HH #include <dune/common/fvector.hh> #include <dune/istl/bvector.hh> #include "norm.hh" -typedef Dune::BlockVector<Dune::FieldVector<double,1> > Vector; - +/** \brief Vector norm induced by filled in low rank linear operator M + * + * \f$M = m^t*m\f$ + * \f$\Vert u \Vert_M = (u, Mu)^{1/2} = (mu,mu)^{1/2}\f$ + * + * \tparam LowRankFactor the type of the factor used to represent the low rank operator + * \tparam Vector the vector type the norm may be applied to + */ +template <class LowRankFactor, class Vector> class FullNorm: public Norm<Vector> { + public: + FullNorm(double alpha, const LowRankFactor &lowRankFactor) : + lowRankFactor_(lowRankFactor), + alpha(alpha) + {} + + //! Compute the norm of the given vector + double operator()(const Vector &v) const + { + Vector r(lowRankFactor_.N()); + r = 0.0; + + lowRankFactor_.umv(v,r); +// for(int row = 0; row < v.size(); ++row) +// r += m[row] * v[row]; + + return sqrt(fabs(alpha*(r*r))); + } + + //! Compute the norm of the difference of two vectors + double diff(const Vector &v1, const Vector &v2) const + { + Vector r(lowRankFactor_.N()); + r = 0.0; + +// Vector temp(v1); +// v1 -= v2; +// lowRankFactor_.umv(temp,r); + + 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]); + + return sqrt(fabs(alpha*(r*r))); + } + + private: + const LowRankFactor& lowRankFactor_; + + const double alpha; + +}; + + +class FullNorm<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> > Vector; public: FullNorm(double alpha, const Vector &m) : m(m), @@ -44,5 +99,4 @@ class FullNorm: public Norm<Vector> const double alpha; }; - #endif