operator/CRS.hh

Go to the documentation of this file.

/* CRS format
 */

#ifndef CRS_hh
#define CRS_hh

#include <sys/types.h>
#include <map>
#include "toolbox/array.hh"
#include "basics/debug.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 ~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+1] = 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;

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



} // namespace concepts

#endif // CRS_hh