operator/denseMatrix.hh

Go to the documentation of this file.

/* A dense matrix
 */

#ifndef denseMatrix_hh
#define denseMatrix_hh

#include "space/element.hh"
#include "compositions.hh"
#include "matrix.hh"
#include <new>              // exception handling for function new

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;

    template<class G>
    DenseMatrix(const Space<G>& spcX, const Space<G>& spcY,
                BilinearForm<F,G>& bf);

    template<class G>
    DenseMatrix(const Space<G>& spc, BilinearForm<F,G>& bf);
    template<class G>
    DenseMatrix(const Space<G>& spcX, const Space<G>& spcY);

    DenseMatrix(int m, int n, 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);
    }

    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);
    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