bdmatrix.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_ISTL_BDMATRIX_HH
#define DUNE_ISTL_BDMATRIX_HH
 
#include <memory>
 
#include <dune/common/rangeutilities.hh>
#include <dune/common/scalarmatrixview.hh>
 
#include <dune/istl/bcrsmatrix.hh>
 
namespace Dune {
  template <class B, class A=std::allocator<B> >
  class BDMatrix : public BCRSMatrix<B,A>
  {
  public:
 
    //===== type definitions and constants
 
    using field_type = typename Imp::BlockTraits<B>::field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    //typedef BCRSMatrix<B,A>::row_type row_type;
 
    typedef typename A::size_type size_type;
 
    static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
 
    BDMatrix() : BCRSMatrix<B,A>() {}
 
    explicit BDMatrix(int size)
      : BCRSMatrix<B,A>(size, size, BCRSMatrix<B,A>::random) {
 
      for (int i=0; i<size; i++)
        this->BCRSMatrix<B,A>::setrowsize(i, 1);
 
      this->BCRSMatrix<B,A>::endrowsizes();
 
      for (int i=0; i<size; i++)
        this->BCRSMatrix<B,A>::addindex(i, i);
 
      this->BCRSMatrix<B,A>::endindices();
 
    }
 
    BDMatrix (std::initializer_list<B> const &list)
    : BDMatrix(list.size())
    {
      size_t i=0;
      for (auto it = list.begin(); it != list.end(); ++it, ++i)
        (*this)[i][i] = *it;
    }
 
    void setSize(size_type size)
    {
      this->BCRSMatrix<B,A>::setSize(size,   // rows
                                     size,   // columns
                                     size);  // nonzeros
 
      for (auto i : range(size))
        this->BCRSMatrix<B,A>::setrowsize(i, 1);
 
      this->BCRSMatrix<B,A>::endrowsizes();
 
      for (auto i : range(size))
        this->BCRSMatrix<B,A>::addindex(i, i);
 
      this->BCRSMatrix<B,A>::endindices();
    }
 
    BDMatrix& operator= (const BDMatrix& other) {
      this->BCRSMatrix<B,A>::operator=(other);
      return *this;
    }
 
    BDMatrix& operator= (const field_type& k) {
      this->BCRSMatrix<B,A>::operator=(k);
      return *this;
    }
 
    template <class V>
    void solve (V& x, const V& rhs) const {
      for (size_type i=0; i<this->N(); i++)
      {
        auto&& xv = Impl::asVector(x[i]);
        auto&& rhsv = Impl::asVector(rhs[i]);
        Impl::asMatrix((*this)[i][i]).solve(xv,rhsv);
      }
    }
 
    void invert() {
      for (size_type i=0; i<this->N(); i++)
        Impl::asMatrix((*this)[i][i]).invert();
    }
 
  private:
 
    // ////////////////////////////////////////////////////////////////////////////
    //   The following methods from the base class should now actually be called
    // ////////////////////////////////////////////////////////////////////////////
 
    // createbegin and createend should be in there, too, but I can't get it to compile
    //     BCRSMatrix<B,A>::CreateIterator createbegin () {}
    //     BCRSMatrix<B,A>::CreateIterator createend () {}
    void setrowsize (size_type i, size_type s) {}
    void incrementrowsize (size_type i) {}
    void endrowsizes () {}
    void addindex (size_type row, size_type col) {}
    void endindices () {}
  };
}  // end namespace Dune
 
#endif