SuperSparseMatrixT.hh 5.34 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
//=============================================================================
//
//  CLASS SuperSparseMatrixT
//
//=============================================================================


#ifndef COMISO_SUPERSPARSEMATRIXT_HH
#define COMISO_SUPERSPARSEMATRIXT_HH


//== INCLUDES =================================================================

#include <iostream>
#include <map>
#include <math.h>
#include <Eigen/Dense>
18

19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
//== FORWARDDECLARATIONS ======================================================

//== NAMESPACES ===============================================================

namespace COMISO {

//== CLASS DEFINITION =========================================================



	      
/** \class SuperSparseMatrixT SuperSparseMatrixT.hh <COMISO/.../SuperSparseMatrixT.hh>

    Brief Description.
  
    A more elaborate description follows.
*/

template <class VType>
class SuperSparseMatrixT
{
public:

  // value type
  typedef VType                                  VT;
  typedef std::pair<unsigned int, unsigned int> PII;

  // iterate over elements
  class iterator
  {
  public:

    // default constructor
    iterator() {}
    // copy constructor
    iterator(const iterator& _it) { m_it_ = _it.m_it_; }
    // map iterator constructor
    iterator(const typename std::map<PII,VType>::iterator& _m_it) { m_it_ = _m_it; }

    // get reference to data
    VT& operator*() { return m_it_->second; }

    // get row and col of value
    unsigned int row() {return m_it_->first.first; }
    unsigned int col() {return m_it_->first.second; }

    // post-increment
    iterator& operator++()    { ++m_it_; return(*this);}
    iterator  operator++(int) { return iterator(++m_it_); }

    bool operator== (const iterator& _it) { return (m_it_ == _it.m_it_);}
    bool operator!= (const iterator& _it) { return (m_it_ != _it.m_it_);}

    // get raw iterator of map
    typename std::map<PII,VType>::iterator& map_iterator() {return m_it_;}

  private:
    typename std::map<PII,VType>::iterator m_it_;
  };

  /// Constructor
  SuperSparseMatrixT(const int _n_rows = 0, const int _n_cols = 0)
  : n_rows_(_n_rows), n_cols_(_n_cols)
  {}
 
  // iterate over non-zeros
  iterator begin() { return iterator(data_.begin()); }
  iterator end()   { return iterator(data_.end());   }
  // erase element
  void     erase( iterator _it) { data_.erase(_it.map_iterator()); }

  // element access
  VT& operator()(const unsigned int _i, const unsigned int _j)
  {  return data_[PII(_i,_j)]; }

  // const element access
  const VT& operator()(const unsigned int _i, const unsigned int _j) const
  {  return data_[PII(_i,_j)]; }

  // get number of stored elements
  unsigned int nonZeros() { return data_.size(); }

  // get dimensions
  unsigned int rows() const { return n_rows_;}
  unsigned int cols() const { return n_cols_;}


  // resize matrix -> delete invalid elements
  void resize(const unsigned int _n_rows, const unsigned int _n_cols)
  {
    n_rows_ = _n_rows;
    n_cols_ = _n_cols;

    // delete out of range elements
    typename std::map<PII, VT>::iterator m_it = data_.begin();
    for(; m_it != data_.end();)
    {
      if(m_it->first.first >=n_rows_ || m_it->first.second >= n_cols_)
        data_.erase(m_it++);
      else
        ++m_it;
    }
  }

  // clear data
  void clear()
  { data_.clear(); }


  // removes all values whose absolut value is smaller than _eps
  void prune(const VT _eps)
  {
    typename std::map<PII, VT>::iterator m_it = data_.begin();
    for(; m_it != data_.end();)
    {
      if(fabs(m_it->second) < _eps)
        data_.erase(m_it++);
      else
        ++m_it;
    }
  }

  // scale matrix by scalar
  void scale(const VT _s)
  {
    typename std::map<PII, VT>::iterator m_it = data_.begin();
    for(; m_it != data_.end(); ++m_it)
      m_it->second *=_s;
  }

  void print()
  {
    iterator it  = begin();
    iterator ite = end();

    std::cerr << "######## SuperSparseMatrix ########" << std::endl;
    std::cerr << "dimension : " << rows()   << " x " << cols() << std::endl;
    std::cerr << "#non-zeros: " << nonZeros() << std::endl;
    for(; it!=ite; ++it)
      std::cerr << "(" << it.row() << "," << it.col() << ") -> " << *it << std::endl;
  }

  void print_eigenvalues()
  {
    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n_rows_, n_cols_);

    iterator it = begin();
    for(; it != end(); ++it)
      A(it.row(),it.col()) = *it;

    Eigen::EigenSolver<Eigen::MatrixXd> eigensolver(A);
    if (eigensolver.info() != Eigen::Success) abort();
    std::cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << std::endl;
//    std::cout << "Here's a matrix whose columns are eigenvectors of A \n"
//    << "corresponding to these eigenvalues:\n"
//    << eigensolver.eigenvectors() << std::endl;
  }

  
private:

  // dimension of matrix
  unsigned int n_rows_;
  unsigned int n_cols_;
  
  typename std::map<PII, VT> data_;
};


//=============================================================================
} // namespace COMISO
//=============================================================================
#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_SUPERSPARSEMATRIXT_C)
#define COMISO_SUPERSPARSEMATRIXT_TEMPLATES
#include "SuperSparseMatrixT.cc"
#endif
//=============================================================================
#endif // COMISO_SUPERSPARSEMATRIXT_HH defined
//=============================================================================