matrix.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MATRIX_HH
#define DUNE_MATRIX_HH
 
#include <vector>
#include <memory>
 
#include <dune/istl/vbvector.hh>
#include <dune/common/ftraits.hh>
 
namespace Dune {
 
  template<class T, class A=std::allocator<T> >
  class Matrix
  {
  public:
 
    typedef typename T::field_type field_type;
 
    typedef T block_type;
 
    typedef A allocator_type;
 
    typedef typename VariableBlockVector<T,A>::window_type row_type;
 
    typedef typename A::size_type size_type;
 
    typedef typename VariableBlockVector<T,A>::Iterator RowIterator;
 
    typedef typename row_type::iterator ColIterator;
 
    typedef typename VariableBlockVector<T,A>::ConstIterator ConstRowIterator;
 
    typedef typename row_type::const_iterator ConstColIterator;
 
    enum {
      blocklevel = T::blocklevel+1
    };
 
    Matrix() : data_(0), cols_(0)
    {}
 
    Matrix(size_type rows, size_type cols) : data_(rows,cols), cols_(cols)
    {}
 
    void setSize(size_type rows, size_type cols) {
      data_.resize(rows,cols);
      cols_ = cols;
    }
 
    RowIterator begin()
    {
      return data_.begin();
    }
 
    RowIterator end()
    {
      return data_.end();
    }
 
    RowIterator beforeEnd ()
    {
      return data_.beforeEnd();
    }
 
    RowIterator beforeBegin ()
    {
      return data_.beforeBegin();
    }
 
    ConstRowIterator begin() const
    {
      return data_.begin();
    }
 
    ConstRowIterator end() const
    {
      return data_.end();
    }
 
    ConstRowIterator beforeEnd() const
    {
      return data_.beforeEnd();
    }
 
    ConstRowIterator beforeBegin () const
    {
      return data_.beforeBegin();
    }
 
    Matrix& operator= (const field_type& t)
    {
      data_ = t;
      return *this;
    }
 
    row_type& operator[](size_type row) {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (row<0)
        DUNE_THROW(ISTLError, "Can't access negative rows!");
      if (row>=N())
        DUNE_THROW(ISTLError, "Row index out of range!");
#endif
      return data_[row];
    }
 
    const row_type& operator[](size_type row) const {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (row<0)
        DUNE_THROW(ISTLError, "Can't access negative rows!");
      if (row>=N())
        DUNE_THROW(ISTLError, "Row index out of range!");
#endif
      return data_[row];
    }
 
    size_type N() const {
      return data_.N();
    }
 
    size_type M() const {
      return cols_;
    }
 
    size_type rowdim() const {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (M()==0)
        DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!");
#endif
      size_type dim = 0;
      for (size_type i=0; i<data_.N(); i++)
        dim += data_[i][0].rowdim();
 
      return dim;
    }
 
    size_type coldim() const {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (N()==0)
        DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!");
#endif
      size_type dim = 0;
      for (size_type i=0; i<data_[0].size(); i++)
        dim += data_[0][i].coldim();
 
      return dim;
    }
 
    size_type rowdim(size_type r) const {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (r<0 || r>=N())
        DUNE_THROW(ISTLError, "Rowdim for nonexisting row " << r << " requested!");
      if (M()==0)
        DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!");
#endif
      return data_[r][0].rowdim();
    }
 
    size_type coldim(size_type c) const {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (c<0 || c>=M())
        DUNE_THROW(ISTLError, "Coldim for nonexisting column " << c << " requested!");
      if (N()==0)
        DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!");
#endif
      return data_[0][c].coldim();
    }
 
    Matrix<T>& operator*=(const field_type& scalar) {
      data_ *= scalar;
      return (*this);
    }
 
    Matrix<T>& operator/=(const field_type& scalar) {
      data_ /= scalar;
      return (*this);
    }
 
    Matrix& operator+= (const Matrix& b) {
#ifdef DUNE_ISTL_WITH_CHECKING
      if(N()!=b.N() || M() != b.M())
        DUNE_THROW(RangeError, "Matrix sizes do not match!");
#endif
      data_ += b.data_;
      return (*this);
    }
 
    Matrix& operator-= (const Matrix& b) {
#ifdef DUNE_ISTL_WITH_CHECKING
      if(N()!=b.N() || M() != b.M())
        DUNE_THROW(RangeError, "Matrix sizes do not match!");
#endif
      data_ -= b.data_;
      return (*this);
    }
 
    Matrix transpose() const {
      Matrix out(N(), M());
      for (size_type i=0; i<M(); i++)
        for (size_type j=0; j<N(); j++)
          out[j][i] = (*this)[i][j];
 
      return out;
    }
 
    friend Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) {
      Matrix<T> out(m1.N(), m2.M());
      out = 0;
 
      for (size_type i=0; i<out.N(); i++ ) {
        for ( size_type j=0; j<out.M(); j++ )
          for (size_type k=0; k<m1.M(); k++)
            out[i][j] += m1[i][k]*m2[k][j];
      }
 
      return out;
    }
 
    template <class X, class Y>
    friend Y operator*(const Matrix<T>& m, const X& vec) {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (m.M()!=vec.size())
        DUNE_THROW(ISTLError, "Vector size doesn't match the number of matrix columns!");
#endif
      Y out(m.N());
      out = 0;
 
      for (size_type i=0; i<out.size(); i++ ) {
        for ( size_type j=0; j<vec.size(); j++ )
          out[i] += m[i][j]*vec[j];
      }
 
      return out;
    }
 
    template <class X, class Y>
    void mv(const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
 
      for (size_type i=0; i<data_.N(); i++) {
        y[i]=0;
        for (size_type j=0; j<cols_; j++)
          (*this)[i][j].umv(x[j], y[i]);
 
      }
 
    }
 
    template<class X, class Y>
    void mtv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
      if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
#endif
 
      for(size_type i=0; i<y.N(); ++i)
        y[i]=0;
      umtv(x,y);
    }
 
    template <class X, class Y>
    void umv(const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
 
      for (size_type i=0; i<data_.N(); i++) {
 
        for (size_type j=0; j<cols_; j++)
          (*this)[i][j].umv(x[j], y[i]);
 
      }
 
    }
 
    template<class X, class Y>
    void mmv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).mmv(x[j.index()],y[i.index()]);
      }
    }
 
    template <class X, class Y>
    void usmv(const field_type& alpha, const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
 
      for (size_type i=0; i<data_.N(); i++) {
 
        for (size_type j=0; j<cols_; j++)
          (*this)[i][j].usmv(alpha, x[j], y[i]);
 
      }
 
    }
 
    template<class X, class Y>
    void umtv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).umtv(x[i.index()],y[j.index()]);
      }
    }
 
    template<class X, class Y>
    void mmtv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).mmtv(x[i.index()],y[j.index()]);
      }
    }
 
    template<class X, class Y>
    void usmtv (const field_type& alpha, const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).usmtv(alpha,x[i.index()],y[j.index()]);
      }
    }
 
    template<class X, class Y>
    void umhv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).umhv(x[i.index()],y[j.index()]);
      }
    }
 
    template<class X, class Y>
    void mmhv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).mmhv(x[i.index()],y[j.index()]);
      }
    }
 
    template<class X, class Y>
    void usmhv (const field_type& alpha, const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
      if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).usmhv(alpha,x[i.index()],y[j.index()]);
      }
    }
 
    //===== norms
 
    typename FieldTraits<field_type>::real_type frobenius_norm () const
    {
      return std::sqrt(frobenius_norm2());
    }
 
    typename FieldTraits<field_type>::real_type frobenius_norm2 () const
    {
      double sum=0;
      for (size_type i=0; i<N(); ++i)
        for (size_type j=0; j<M(); ++j)
          sum += data_[i][j].frobenius_norm2();
      return sum;
    }
 
    typename FieldTraits<field_type>::real_type infinity_norm () const
    {
      double max=0;
      for (size_type i=0; i<N(); ++i) {
        double sum=0;
        for (size_type j=0; j<M(); j++)
          sum += data_[i][j].infinity_norm();
        max = std::max(max,sum);
      }
      return max;
    }
 
    typename FieldTraits<field_type>::real_type infinity_norm_real () const
    {
      double max=0;
      for (size_type i=0; i<N(); ++i) {
        double sum=0;
        for (size_type j=0; j<M(); j++)
          sum += data_[i][j].infinity_norm_real();
        max = std::max(max,sum);
      }
      return max;
    }
 
    //===== query
 
    bool exists (size_type i, size_type j) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (i<0 || i>=N()) DUNE_THROW(ISTLError,"row index out of range");
      if (j<0 || i>=M()) DUNE_THROW(ISTLError,"column index out of range");
#endif
      return true;
    }
 
  protected:
 
    VariableBlockVector<T,A> data_;
 
    size_type cols_;
  };
 
} // end namespace Dune
 
#endif