operator/hashedSMatrix.hh

Go to the documentation of this file.

/* SMatrix
 * sparse matrices
 */

#ifndef hSparseMatrix_hh
#define hSparseMatrix_hh

#include <iostream>
#include <fstream>
#if __GNUC__ == 2
#  include <float.h>
#  define EPS DBL_EPSILON
#else
#  include <limits>
#  define EPS std::numeric_limits<double>::epsilon()
#endif
#include "basics/typedefs.hh"
#include "basics/outputOperator.hh"
#include "toolbox/pool.hh"
#include "hashedSMatrixIterator.hh"

namespace concepts {

  // forward declaration
  template<class F>
  class Matrix;

  template<class F>
  class DenseMatrix;

  // **************************************************** HashedSparseMatrix **

  template <class T>
  class HashedSparseMatrix {
  public:
    typedef T value_type;
    typedef typename Realtype<T>::type r_type;
    typedef typename Cmplxtype<T>::type c_type;

    typedef _HashedSMatrix_iterator<T, T&, T*>             iterator;
    typedef _HashedSMatrix_iterator<T, const T&, const T*> const_iterator;

    struct Value {
      Value* lnk;
      uint   idx;
      T      val;

      Value(Value* l, int i, T v = 0.0) { lnk = l; idx = i; val = v; }
      void* operator new(size_t sz) { return new char[sz]; }
      void* operator new
      (size_t, Pool<typename HashedSparseMatrix<T>::Value>& pool) { 
        return pool.alloc(); }
      Value& operator*=(const T factor) {
        val *= factor;
        return *this;
      }
    };

    HashedSparseMatrix(uint r, uint c, uint h);
    ~HashedSparseMatrix();

    inline Value** Matrix() const { return matrix; } // new
    inline uint HashBits() const { return hashBits; } // new

    T& operator()(const uint r, const uint c);
    T operator()(const uint r, const uint c) const;

    const uint rows() const { return nofRows; }
    const uint cols() const { return nofCols; }

    template<typename F, typename G>
    void operator()(const F f[], G g[]) const;

    HashedSparseMatrix<T>& operator*=(const T factor);

    template<typename F, typename G>
    void transpMult(const F f[], G g[]) const;

    uint used() const { return _used; }
    float memory() const { 
      return sizeof(HashedSparseMatrix<T>) + pool.memory() 
        - sizeof(pool); 
    }

    void write(std::ostream& ofs) const;
    void outputMatlab(std::ostream& ofs) const;
    void outputSparseQR(std::ostream& ofs) const;

    void multiply(const HashedSparseMatrix<T>* const fact,
                  concepts::Matrix<T>& dest) const;
    void multiply(const concepts::Matrix<T>* const fact,
                  concepts::Matrix<T>& dest) const;

    const_iterator begin(uint row = 0) const;
    iterator begin(uint row = 0);
    static iterator end();

    std::ostream& info(std::ostream& os) const;

    void compress(Real threshold = EPS);
  private:
    Pool<Value> pool;
    Value**     matrix;
    uint        _used;

    uint        nofRows;

    uint        nofCols;

    uint        hashBits;
    uint        hashMsk;

    HashedSparseMatrix(const HashedSparseMatrix<T>&);
    void operator =(const HashedSparseMatrix<T>&);

    void dropEntry_(Value*& v, Value*& vPrev, const uint r, const uint c);
  };

  template<typename T>
  std::ostream& operator<<(std::ostream& os, const HashedSparseMatrix<T>& o) {
#ifdef DEBUG
    os << std::flush;
#endif
    return o.info(os);
  }

  template<class T>
  template<class F, class G>
  void HashedSparseMatrix<T>::transpMult(const F f[], G g[]) const {
    uint r = nofRows << hashBits;

    while (r--)
      for(Value* v = matrix[r]; v != NULL; v = v->lnk) 
        g[v->idx] += v->val * f[r >> hashBits];
  }

  template<class T>
  void HashedSparseMatrix<T>::multiply(const HashedSparseMatrix<T>* const fact,
                                       concepts::Matrix<T>& dest) const {
    conceptsAssert(nofCols == fact->nofRows, Assertion());
    conceptsAssert(nofRows == dest.dimX(), Assertion());
    conceptsAssert(fact->nofCols == dest.dimY(), Assertion());
    uint r = nofRows << hashBits;

    while (r--) {
      for (Value* v= matrix[r]; v != 0; v = v->lnk) {
        for (uint k = 0; k < fact->nofCols; ++k) {
          Value* w = fact->matrix[v->idx << fact->hashBits | (k & fact->hashMsk)];
          for (; w != 0 && w->idx != k; w = w->lnk);
          if (w != 0) {
            T sum = v->val * w->val;
            if (sum != (T)0)
              dest(r >> hashBits, k) += sum;
          }
        } // for k
      } // for v
    } // while
  }


  template<class T>
  void HashedSparseMatrix<T>::multiply(const concepts::Matrix<T>* const fact,
                                       concepts::Matrix<T>& dest) const {
    conceptsAssert(nofCols == fact->dimX(), Assertion());
    conceptsAssert(nofRows == dest.dimX(), Assertion());
    conceptsAssert(fact->dimY() == dest.dimY(), Assertion());
    uint r = nofRows << hashBits;
    const uint BnofCols = fact->dimY();
    while (r--) {
      for (Value* v= matrix[r]; v != 0; v = v->lnk) {
        for (uint k = 0; k < BnofCols; ++k) {
          T sum = v->val * (*fact)(v->idx,k);
          if (sum != (T)0)
            dest(r >> hashBits, k) += sum;
        } // for k
      } // for v
    } // while
  }

  template<class T>
  template<class F, class G>
  void HashedSparseMatrix<T>::operator()(const F f[], G g[]) const {
    G sum = 0.0;
    uint r = nofRows << hashBits;

    while (r--) {
      for(Value* v = matrix[r]; v != NULL; v = v->lnk) 
        sum += v->val * f[v->idx];
      if (!(r & hashMsk)) { 
        g[r >> hashBits] = sum;
        sum = 0.0;
      }
    }
  }
  
} // namespace concepts

#endif // hSparseMatrix_hh