matrix.hh

Go to the documentation of this file.

 
#ifndef matrix_hh
#define matrix_hh
 
#include "compositions.hh"
#include "matrixIterator.hh"
#include "space/elementPairs.hh"
#include "toolbox/sequence.hh"
 
#include <type_traits>
#include <typeinfo>
 
#define ApplyMatrix_D 0
 
namespace concepts {
 
  // forward declarations
  template<typename F, typename G>
  class BilinearForm;
 
  class InOutParameters;
 
  // **************************************************************** Matrix **
 
  template<class F>
  class Matrix : public Operator<F> {
  public:
    typedef F                                       value_type;
    typedef typename Realtype<F>::type              r_type;
    typedef typename Cmplxtype<F>::type             c_type;
    typedef typename std::conditional<std::is_same<typename Realtype<F>::type, F>::value ,
      typename Realtype<F>::type, typename Cmplxtype<F>::type >::type d_type;
 
    typedef _Matrix_iterator<F, F&, F*>             iterator;
    typedef _Matrix_iterator<F, const F&, const F*> const_iterator;
 
    Matrix(uint nofRows, uint nofCols) : Operator<F>(nofRows, nofCols) {}
 
    const uint nofRows() const { return this->dimX(); }
    const uint nofCols() const { return this->dimY(); }
 
    iterator begin(uint r = 0) { return iterator(*this, r); }
    iterator end() { return iterator(*this, nofRows(), 0); }
    const_iterator begin(uint r = 0) const { return const_iterator(const_cast<Matrix<F>&>(*this), r); }
    const_iterator end() const { 
      return const_iterator(const_cast<Matrix<F>&>(*this), nofRows(), 0);
    }
 
    virtual void set(const uint i, const uint j,
                     const F value,
                     const bool use_threshold = false,
                     const Real threshold_value = 1e-8);
 
    virtual void add(const uint i, const uint j,
                     const F value,
                     const bool use_threshold = false,
                     const Real threshold_value = 1e-8);
    
    template<class G>
    static void assembly(Matrix<F>& dest, const Space<G>& spc,
                         const BilinearForm<F,G>& bf, const Real threshold = 0.0);
    template<class G>
    static void assembly(Matrix<F>& dest, Scan<Element<G> > * sc,
                         const BilinearForm<F,G>& bf, const Real threshold = 0.0);
    template<class G>
    static void assembly(Matrix<F>& dest,
                         const Sequence<ElementWithCell<G> * > seq,
                         const BilinearForm<F,G>& bf, const Real threshold = 0.0);
    template<class G>
    static void assembly(Matrix<F>& dest,
                         const Space<G>& spc, const Sequence<bool>& seq, 
                         const BilinearForm<F,G>& bf, const Real threshold = 0.0);
 
    template<class G>
    static void assembly(Matrix<F>& dest,
                         const Space<G>& spcX, const Space<G>& spcY,
                         const BilinearForm<F,G>& bf, const Real threshold = 0.0, const bool single = false);
 
    template<class G>
    static void assembly(Matrix<F>& dest,
                         const Sequence<ElementWithCell<G> * > seqX,
                         const Space<G>& spcY,
                         const BilinearForm<F,G>& bf, const Real threshold = 0.0);
 
    template<class G>
    static void assembly(Matrix<F>& dest, const BilinearForm<F,G>& bf,
                         const ElementPairList<G>& pairs);
    template<class G>
    static void assembly(Matrix<F>& dest,
                         const Space<G>& spcX, const Space<G>& spcY, const Sequence<bool>& seq, 
                         const BilinearForm<F,G>& bf, const Real threshold = 0.0, const bool single = false);
 
 
    virtual F operator()(const uint i, const uint j) const = 0;
    virtual F& operator()(const uint i, const uint j) = 0;
 
    virtual void operator()(const Function<r_type>& fncY,
                            Function<F>& fncX) = 0;
    virtual void operator()(const Function<c_type>& fncY,
                            Function<c_type>& fncX) = 0;
 
    virtual bool operator==(const Matrix<F>& otherMat) const {
      uint rows = nofRows();
      uint cols = nofCols();
 
      if( rows != otherMat.nofRows() || cols != otherMat.nofCols() )
        return false;
 
      for(uint i=0; i<rows; i++)
        for(uint j=0; j<cols; j++)
          if( (*this)(i,j)!=otherMat(i,j) )
            return false;
 
      return true;
    }
 
    virtual void transpMult(const Vector<r_type>& fncY,
                            Vector<F>& fncX) = 0;
    virtual void transpMult(const Vector<c_type>& fncY,
                            Vector<c_type>& fncX) = 0;
 
    static void setTimings(InOutParameters* timings);
    static bool timings();
  private:
    static InOutParameters* timings_;
    static uint timeCntr_;
  };
  
  // ****************************************************** getNumberofRows **
  
  template<class F>
  uint getNumberofRows(F& m) { return m.nofRows(); }
  
 
  // **************************************************************** apply **
 
  template<class F, class H, class I>
  void apply(Operator<F>& op, const Matrix<H>& mX, Matrix<I>& mY) {
    // number of rows
    const uint n = op.dimX();
 
    // number of rows should be the same as for the solver
    conceptsAssert(mX.dimX() == n, Assertion());
    conceptsAssert(mY.dimX() == n, Assertion());
    // the resulting matrix should have the same number of columns
    // then the applicated matrix
    conceptsAssert(mX.dimY() == mY.dimY(), Assertion());
 
    // loop over columns
    Vector<H> fncX(mX.dimX());
    Vector<I> fncY(mY.dimX());
    I& f = fncY(0);
    for (uint c = 0; c < mX.dimY(); ++c) {
      // Write column of applicated matrix to vector.
      // Optimal would be a replacement by a method of the matrix.
      for (uint r = 0; r < n; ++r)
        fncX(r) = mX(r, c);
 
      DEBUGL(ApplyMatrix_D, "fncX = " << fncX);
      op(fncX, fncY);
      DEBUGL(ApplyMatrix_D, "fncY = " << fncY);
 
      // Add vector to column of resulting matrix, could be
      // Optimal would be a replacement by a method of the matrix.
      for (uint r = 0; r < n; ++r)
        if ((f = fncY(r)) != 0.0)
          mY(r, c) += f;
      DEBUGL(ApplyMatrix_D, "finished column " << c);
    }
  }
 
} // namespace concepts
 
#endif // matrix_hh