denseMatrix.hh

Go to the documentation of this file.

 
#ifndef denseMatrix_hh
#define denseMatrix_hh
 
#include "space/element.hh"
#include "compositions.hh"
#include "matrix.hh"
#include <new>              // exception handling for function new
 
#include <type_traits>
#include <typeinfo>
 
namespace concepts {
 
  // forward declaration
  template<class F>
  class Space;
 
  template<class F, class G>
  class BilinearForm;
 
  // *********************************************************** DenseMatrix **
 
  template<class F>
  class DenseMatrix : public Matrix<F>, public ElementMatrix<F> {
  public:
    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;
 
    template<class G>
    DenseMatrix(const Space<G>& spcX, const Space<G>& spcY,
                const BilinearForm<F,G>& bf, const bool single=true);
 
    template<class G>
    DenseMatrix(const Space<G>& spc, const BilinearForm<F,G>& bf);
    template<class G>
    DenseMatrix(const Space<G>& spc, const BilinearForm<F,G>& bf,
                const Sequence<ElementWithCell<G> * > & seq);
    template<class G>
    DenseMatrix(const Space<G>& spcX, const Space<G>& spcY);
 
    DenseMatrix(int m = 0, int n = 0, bool transpose = false)
      : Matrix<F>(m, n), ElementMatrix<F>(transpose ? n : m, transpose ? m : n)
    {
      if(transpose)
        this->transpose();
    }
 
    DenseMatrix(const LiCo<F>& L);
    DenseMatrix(const Compose<F>& L);
    DenseMatrix(const DenseMatrix<F>& m);
    DenseMatrix(const Matrix<F>& m, bool t = false);
    DenseMatrix(const Vector<F>& v, bool t = false);
 
    virtual void operator()(const Function<r_type>& fncY,
                            Function<F>& fncX);
    virtual void operator()(const Function<c_type>& fncY,
                            Function<c_type>& fncX);
 
    template<class H, class I>
    void operator()(const Vector<H>& fncY, Vector<I>& fncX);
 
    virtual void transpMult(const Vector<r_type>& fncY,
                            Vector<F>& fncX);
    virtual void transpMult(const Vector<c_type>& fncY,
                            Vector<c_type>& fncX);
 
    virtual F operator()(const uint i, const uint j) const {
      return ElementMatrix<F>::operator()(i,j);
    }
 
    virtual F& operator()(const uint i, const uint j) {
      return ElementMatrix<F>::operator()(i,j);
    }
 
    virtual DenseMatrix<F>& operator=(const F& v) 
    {
      ElementMatrix<F>::operator=(v);
      return *this; 
    }
    
    virtual void resize(uint m, uint n) {
      ElementMatrix<F>::resize(this->dimX_ = m, this->dimY_ = n);
    }
 
    virtual void setTranspose(bool t) {
      if(this->isTranspose() != t) std::swap(this->dimX_,this->dimY_);
      ElementMatrix<F>::setTranspose(t);
    }
 
    void addInto(Matrix<F>& dest, const F fact,
                 const uint rowoffset = 0, const uint coloffset = 0) const;
    
    template<class T>
    void addInto(Matrix<T>& dest, const T fact,
                 const uint rowoffset = 0, const uint coloffset = 0) const
    {
      conceptsAssert(dest.dimX() >= this->dimX() + rowoffset, Assertion());
      conceptsAssert3(dest.dimY() >= this->dimY() + coloffset, Assertion(),
                      "dest.dimY() = " << dest.dimY() << ", dimY() = " <<
                      this->dimY() << ", coloffset = " << coloffset);
      for (uint i = 0; i < this->dimX(); ++i)
        for (uint j = 0; j < this->dimY(); ++j) {
          const F a = operator()(i,j);
          dest(i+rowoffset,j+coloffset) += a*fact;
        }
    };
 
    template<class T>
    void addIntoT(Matrix<T>& dest, const T fact,
                  const uint rowoffset = 0, const uint coloffset = 0) const {
      conceptsAssert(dest.dimY() >= this->dimX() + coloffset, Assertion());
      conceptsAssert3(dest.dimX() >= this->dimY() + rowoffset, Assertion(),
                      "dest.dimY() = " << dest.dimY() << ", dimY() = " <<
                      this->dimY() << ", coloffset = " << coloffset);
      for (uint i = 0; i < this->dimX(); ++i)
        for (uint j = 0; j < this->dimY(); ++j) {
          const F a = operator()(i,j);
          dest(j+rowoffset,i+coloffset) += a*fact;
        }
    };
 
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };
 
  template<class F>
  template<class H, class I>
  void DenseMatrix<F>::operator()(const Vector<H>& fncY, Vector<I>& fncX)
  {
    conceptsAssert(ElementMatrix<F>::m() == fncX.dim(), Assertion());
    conceptsAssert(ElementMatrix<F>::n() == fncY.dim(), Assertion());
    const F* m = (const F*)this->getData();
    I* x(fncX); 
    for (uint i = 0; i < ElementMatrix<F>::m(); ++i) {
      H* y(fncY); 
      I fi = 0.0;
      for (uint j = 0; j < ElementMatrix<F>::n(); ++j)
        fi += *m++ * *y++;
      *x++ = fi;
    }
  }
 
} // namespace concepts
 
#endif // senseMatrix_hh