CRS.hh

Go to the documentation of this file.

 
#ifndef CRS_hh
#define CRS_hh
 
#include <sys/types.h>
#include <map>
#include "toolbox/array.hh"
#include "basics/debug.hh"
#include "operator/hashedSMatrix.hh"
 
#define CRS_rowSorting_D 0
#define CCS_rowSorting_D 0
 
namespace concepts {
 
  // ******************************************************** CRSConvertable **
 
  template<class F>
  class CRSConvertable {
  public:
    virtual uint used() const = 0;
    virtual void convertCRS(F* a, int* asub, int* xa) const = 0;
    virtual void convertCCS(F* a, int* asub, int* xa) const = 0;
    virtual void convertIJK(F* a, int* irn, int* jcn) const = 0;
 
    virtual ~CRSConvertable() {}
  };
 
  // ************************************************* sparseLineToArrays **
  
  template<typename F>
  void sparseLineToArrays(std::map<int, F>& line, F* a, int* asub) {
    uint nnz = 0;
    typename std::map<int, F>::const_iterator j = line.begin();
    for(; j != line.end(); ++j) {
      a[nnz] = j->second;       // value
      asub[nnz++] = j->first;   // column number
    }
  }
  
  // ************************************************* convertCRS_rowSorting **
 
  template<class F>
  void convertCRS_rowSorting(F& m, typename F::value_type* a, int* asub,
                             int* xa)
  {
    DEBUGL(CRS_rowSorting_D, "CRS format of " << m);
 
    typedef typename F::value_type value;
 
    int nnz = 0;
    uint r = 0;
    xa[0] = 0;
    std::map<int, value> row;  // temporary, for sorting inside row
 
    typename F::const_iterator i = m.begin();
    DEBUGL(CRS_rowSorting_D, "i = " << i);
    for(; i != m.end(); ++i) {
      DEBUGL((CRS_rowSorting_D && i != m.begin()), "i = " << i);
      if (i.row() > r) {
        sparseLineToArrays(row, a + nnz, asub + nnz);
        nnz += row.size();
        row.clear();
        for(; r < i.row();)
          xa[++r] = nnz;
      }
      row[i.col()] = *i;
    }
    sparseLineToArrays(row, a + nnz, asub + nnz);
    xa[++r] = nnz + row.size();
    for(; r < getNumberofRows(m) ; )
      xa[++r] = nnz + row.size();
  }
 
  // ************************************************* convertCCS_rowSorting **
 
  template<class F>
  void convertCCS_rowSorting(F& m, typename F::type* a, int* asub, int* xa) {
    DEBUGL(CCS_rowSorting_D, "CCS format of " << m);
 
    //typedef typename F::type value; this is not used somewhere
 
    uint nnz = 0;
    uint n = m.nofCols();   // number of columns
 
    for(uint c = 0; c < n; ++c)
      xa[c] = 0;
    // First loop over all entrances and count the number of entrances
    // per column (instead of last one).
    for(typename F::const_iterator i = m.begin(); i != m.end(); ++i) {
      if (i.col() < n-1) ++xa[i.col()+1];
      ++nnz;
    }
 
    // cumulated sum -> first index of each column
    for(uint c = 2; c < n; ++c)
      xa[c] += xa[c-1];
    xa[n] = nnz;
 
    // copy it to this array for the current index in each column
    Array<int> idx(xa, n+1);
    DEBUGL(CCS_rowSorting_D, "col_ptr = " << idx);
    for(typename F::const_iterator i = m.begin(); i != m.end(); ++i) {
      int &j = idx[i.col()];
      a[j] = *i       ;      // nzval
      asub[j] = i.row();     // rowind
      DEBUGL(CCS_rowSorting_D, "write (" << i.row() << ", " << i.col() <<
       ", " << *i << ") to idx " << j);
      ++j;
    }
  }
 
  template<class F>
  void convertIJK_unSorted(F& m, typename F::type* a, int* irn, int* jcn) {
    //uint nnz = 0; //nothing happens here?
 
  }
 
} // namespace concepts
 
#endif // CRS_hh