operator/matrix.hh

Go to the documentation of this file.

/* interface for a matrix
 */

#ifndef matrix_hh
#define matrix_hh

#include "compositions.hh"
#include "matrixIterator.hh"
#include "space/elementPairs.hh"

#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 _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(*this, r); }
    const_iterator end() const { 
      return const_iterator(*this, nofRows(), 0);
    }

    template<class G>
    static void assembly(Matrix<F>& dest, const Space<G>& spc,
                         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,
                         BilinearForm<F,G>& bf, const Real threshold = 0.0);
    template<class G>
    static void assembly(Matrix<F>& dest, BilinearForm<F,G>& bf,
                         const ElementPairList<G>& pairs);

    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_;
  };

  // **************************************************************** 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