sparseMatrix.hh

/* @file operator/sparseMatrix.hh
 
   A sparse matrix
 */
 
#ifndef sparseMatrix_hh
#define sparseMatrix_hh
 
#include <map>
#if __GNUC__ == 2
#include <float.h>
#define EPS DBL_EPSILON
#else
#include <limits>
#define EPS std::numeric_limits<double>::epsilon()
#endif
#include <memory>
#include "compositions.hh"
#include "matrix.hh"
#include "CRS.hh"
#include "hashedSMatrix.hh"
#include "matrixMult.hh"
#include "basics/debug.hh"
#include "toolbox/sequence.hh"
 
#include "bilinearForm.hh"
 
#include <type_traits>
#include <typeinfo>
 
#define SparseAddInto_D 0
 
namespace concepts {
 
  // forward declaration
  template<class F, class G>
  class BilinearForm;
 
  // forward declaration
  //template<class F, class G>
  //class BilinearFormContainer;
 
  // ********************************************************** SparseMatrix **
 
  template<class F>
  class SparseMatrix : public Matrix<F>, public CRSConvertable<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;
 
    typedef _HashedSMatrix_iterator<F, F&, F*>             iterator;
    typedef _HashedSMatrix_iterator<F, const F&, const F*> const_iterator;
 
    template<class G>
    SparseMatrix(const Space<G>& spcX, const Space<G>& spcY) 
      : Matrix<F>(spcX.dim(), spcY.dim())
      , nX_(spcX.dim())
      , nY_(spcY.dim()), 
      m_(new HashedSparseMatrix<F>(spcX.dim(), spcY.dim(), 2)) 
    { }
 
    SparseMatrix(uint nofrows, uint nofcols) 
      : Matrix<F>(nofrows, nofcols)
      , nX_(nofrows)
      , nY_(nofcols), 
      m_(new HashedSparseMatrix<F>(nofrows, nofcols, 2)) 
    { }
    SparseMatrix(uint dim = 0) 
      : Matrix<F>(dim, dim)
      , nX_(dim)
      , nY_(dim), 
      m_(new HashedSparseMatrix<F>(dim, dim, 2)) 
    { }
 
    template<class G>
    SparseMatrix(const Space<G>& spcX, const Space<G>& spcY, 
                 const BilinearForm<F,G>& bf, const Real eps = 0.0,
                 const bool single = false);
 
    template<class G>
    SparseMatrix(const Space<G>& spcX, const Space<G>& spcY,
                 const Sequence<bool> & seq,
                 const BilinearForm<F,G>& bf, const Real eps = 0.0,
                 const bool single = false);
 
    template<class G>
    SparseMatrix(const Space<G>& spcX, const Space<G>& spcY, 
                 const BilinearForm<F,G>& bf,
                 const Sequence<ElementWithCell<G> * > & seq,
                 const Real eps = 0.0);
 
    template<class G>
    SparseMatrix(const Space<G>& spc, const BilinearForm<F,G>& bf,
                 const Real eps = 0.0);
 
 
    SparseMatrix(const Space<typename Realtype<F>::type >& spc, const BilinearFormContainer<F> bf,
                 const Real eps = 0.0);
 
    template<class G>
    SparseMatrix(const Space<G>& spc, 
                 const Sequence<bool> & seq, const BilinearForm<F,G>& bf,
                 const Real eps = 0.0);
 
    template<class G>
    SparseMatrix(const Space<G>& spc, const BilinearForm<F,G>& bf,
                 const Sequence<ElementWithCell<G> * > & seq,
                 const Real eps = 0.0);
    
    SparseMatrix(const SparseMatrix<F>& m, bool t = false);
 
    template<class G>
    SparseMatrix(const Space<G>& spc, const F* const v, const int i = -1);
 
    template<class G>
    SparseMatrix(const Space<G>& spc,
                 const Vector<F>& x, const Vector<F>& y);
 
    template<class H>
    SparseMatrix(Operator<H>& A, bool slow = false);
 
    void operator=(const SparseMatrix<F>&);
 
    template<class H>
    SparseMatrix(const SparseMatrix<H>& fncX, F fnc(const H&));
    template<class H>
    SparseMatrix(const SparseMatrix<H>& fncX, const F& fnc(const H&));
 
    template<class H>
    SparseMatrix(const SparseMatrix<H>& fncX); 
 
 
    virtual ~SparseMatrix(){};
 
    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) {
      (*m_)((const H*)fncY, (I*)fncX);
    }
 
    const_iterator begin(uint row = 0) const;
    iterator begin(uint row = 0);
    const_iterator end() const;
 
    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 (const_cast<const HashedSparseMatrix<F>&>(*m_))(i,j); }
    virtual F& operator()(const uint i, const uint j) { return (*m_)(i,j); }
    
    SparseMatrix<F>& operator*=(const F factor) {
      (*m_) *= factor;
      return *this;
    }
 
    virtual void resize(uint m, uint n) {
      this->dimX_ = nX_ = m; this->dimY_ = nY_ = n;
      m_.reset(new HashedSparseMatrix<F>(m, n, 2));
    }
 
    void write(char* fname) const;
 
    bool storeMatlab(const std::string filename, const std::string name = "",
                     bool append = false) const;
 
    const HashedSparseMatrix<F>* m() const {
      return const_cast<const HashedSparseMatrix<F>*>(m_.get()); }
 
    template<class H, class I>
    void addInto(Matrix<H>& dest, const I fact,
                 const uint rowoffset = 0, const uint coloffset = 0) const;
    template<class H, class I>
    void addIntoT(Matrix<H>& dest, const I fact,
                  const uint rowoffset = 0, const uint coloffset = 0) const;
 
    template<class H, class I>
    void add(const Vector<H>& v, const I fact,
             const uint rowoffset = 0, const uint coloffset = 0);
 
    template<class H, class I>
    void addT(const Vector<H>& v, const I fact,
              const uint rowoffset = 0, const uint coloffset = 0);
 
    void multiply(const SparseMatrix<F>& fact, Matrix<F>& dest) const {
      m_->multiply(fact.m(), dest);
    }
    template<class H>
    void multiply(const H& fact, Matrix<F>& dest) const {
      matrixMultiplyRowSorting(*this, fact, dest);
    }
    virtual uint used() const { return m_->used(); }
    float memory() const { return sizeof(SparseMatrix<F>) + m_->memory(); }
    void symmetrize();
 
    void compress(Real threshold = EPS) {
      m_->compress(threshold);
    }
 
    void histogram(std::map<int, uint>& hist) const;
 
    void copy(const SparseMatrix<F>& n);
 
    virtual void convertCRS(F* a, int* asub, int* xa) const;
    virtual void convertCCS(F* a, int* asub, int* xa) const;
    virtual void convertIJK(F* , int* , int* ) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    uint              nX_;
 
    uint              nY_;
 
    std::unique_ptr<HashedSparseMatrix<F> > m_;
 
    static uint storeMatlabCounter_;
 
    void copyConstructor_(const SparseMatrix<F>& m, bool t);
 
    virtual void apply_(const Vector<F>& fncY, Vector<F>& fncX) {
      conceptsAssert(0, Assertion());
      (*m_)((const F*)fncY, (F*)fncX);
    }
  };
 
  template<class F>
  template<class H>
  SparseMatrix<F>::SparseMatrix(const SparseMatrix<H>& m, F fnc(const H&)) 
    : Matrix<F>(m.dimX(), m.dimY())
    , nX_(m.dimX()), nY_(m.dimY())
    , m_(new HashedSparseMatrix<F>(nX_, nY_, 2))
  {
    int hashBits = m.m()->HashBits();
    int p = 1 << hashBits;
    typename HashedSparseMatrix<H>::Value** matrix = m.m()->Matrix();
    
    for (uint r = 0; r < nX_; ++r)
      for (int c_mod = 0; c_mod < p; ++c_mod)
        for (typename HashedSparseMatrix<H>::Value* v =
               matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
          (*m_)(r, v->idx) = fnc(v->val);
  }
 
  template<class F>
  template<class H>
  SparseMatrix<F>::SparseMatrix(const SparseMatrix<H>& m) 
    : Matrix<F>(m.dimX(), m.dimY())
    , nX_(m.dimX()), nY_(m.dimY())
    , m_(new HashedSparseMatrix<F>(nX_, nY_, 2))
  {
    int hashBits = m.m()->HashBits();
    int p = 1 << hashBits;
    typename HashedSparseMatrix<H>::Value** matrix = m.m()->Matrix();
    
    for (uint r = 0; r < nX_; ++r)
      for (int c_mod = 0; c_mod < p; ++c_mod)
        for (typename HashedSparseMatrix<H>::Value* v =
               matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
          (*m_)(r, v->idx) = (F)(v->val);
  }
 
  template<class F>
  template<class H>
  SparseMatrix<F>::SparseMatrix(const SparseMatrix<H>& m,
                                const F& fnc(const H&))
    : Matrix<F>(m.dimX(), m.dimY())
    , nX_(m.dimX()), nY_(m.dimY())
    , m_(new HashedSparseMatrix<F>(nX_, nY_, 2))
  {
    int hashBits = m.m()->HashBits();
    int p = 1 << hashBits;
    typename HashedSparseMatrix<H>::Value** matrix = m.m()->Matrix();
    
    for (uint r = 0; r < nX_; ++r)
      for (int c_mod = 0; c_mod < p; ++c_mod)
        for (typename HashedSparseMatrix<H>::Value* v =
               matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
          (*m_)(r, v->idx) = fnc(v->val);
  }
 
 
  template<class F>
  template<class H, class I>
  void SparseMatrix<F>::addInto(Matrix<H>& dest, const I fact,
                                const uint rowoffset, 
                                const uint coloffset) const
  {
    conceptsAssert(dest.dimX() >= this->dimX() + rowoffset, Assertion());
    conceptsAssert3(dest.dimY() >= this->dimY() + coloffset, Assertion(),
                    "dest.dimY() = " << dest.dimY() << ", dimY() = " <<
                    this->dimY() << ", coloffset = " << coloffset);
    DEBUGL(SparseAddInto_D, "fact = " << fact <<
           ", rowoffset = " << rowoffset <<
           ", coloffset = " << coloffset);
    if (fact != 0.0) {
 
      int hashBits = m_->HashBits();
      int p = 1 << hashBits;
      uint n = this->dimX();
      typename HashedSparseMatrix<F>::Value** matrix = m_->Matrix();
 
      for (uint r = 0; r < n; ++r)
        for (int c_mod = 0; c_mod < p; ++c_mod)
          for (typename HashedSparseMatrix<F>::Value* v =
                 matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
            if (v->val != 0.0)
              dest(r + rowoffset, v->idx + coloffset) += (v->val * fact);
    }
  }
 
  template<class F>
  template<class H, class I>
  void SparseMatrix<F>::addIntoT(Matrix<H>& dest, const I fact,
                                 const uint rowoffset, 
                                 const uint coloffset) const
  {
    conceptsAssert3(dest.dimY() >= this->dimX() + coloffset, Assertion(),
                    "dest.dimY() = " << dest.dimX() << ", dimX() = " <<
                    this->dimX() << ", coloffset = " << coloffset);
    conceptsAssert3(dest.dimX() >= this->dimY() + rowoffset, Assertion(),
                    "dest.dimX() = " << dest.dimX() << ", dimY() = " <<
                    this->dimY() << ", rowoffset = " << rowoffset);
    DEBUGL(SparseAddInto_D, "fact = " << fact <<
           ", rowoffset = " << rowoffset <<
           ", coloffset = " << coloffset);
    if (fact != 0.0) {
      int hashBits = m_->HashBits();
      int p = 1 << hashBits;
      uint n = this->dimX();
      typename HashedSparseMatrix<F>::Value** matrix = m_->Matrix();
 
      for (uint r = 0; r < n; ++r)
        for (int c_mod = 0; c_mod < p; ++c_mod)
          for (typename HashedSparseMatrix<F>::Value* v =
                 matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
            if (v->val != 0.0)
              dest(v->idx + rowoffset, r + coloffset) += (v->val*fact);
    }
  }
 
  template<class F>
  template<class H, class I>
  void SparseMatrix<F>::add(const Vector<H>& v, const I fact,
                            const uint rowoffset, const uint coloffset) 
  {
    conceptsAssert(this->dimX() >= v.dim() + rowoffset, Assertion());
    conceptsAssert(this->dimY() >= 1       + coloffset, Assertion());
 
    if (fact != 0.0) {
      const H* h = (const H*)v;
      uint row = rowoffset;
      for(uint i = v.dim(); i--; )
        (*this)(row++, coloffset) += *h++ * fact;
    }
  }
 
  template<class F>
  template<class H, class I>
  void SparseMatrix<F>::addT(const Vector<H>& v, const I fact,
                             const uint rowoffset, const uint coloffset) 
  {
    conceptsAssert(this->dimX() >= 1       + rowoffset, Assertion());
    conceptsAssert(this->dimY() >= v.dim() + coloffset, Assertion());
 
    if (fact != 0.0) {
      const H* h = (const H*)v;
      uint col = coloffset;
      for(uint i = v.dim(); i--; )
        (*this)(rowoffset, col++) += *h++ * fact;
    }
  }
 
} // namespace concepts
 
#endif // sparseMatrix_hh