text: #ifndef _BMCSB_H
  #define _BMCSB_H
  
  #include <cmath>
  #include <vector>
  #include <algorithm>
  #include <numeric>		// for std:accumulate()
  #include <limits>		// C++ style numeric_limits<T>
  #include "csc.h"
  #include "mortoncompare.h"
  
  using namespace std;
  
  void SSEspmv(const double * __restrict V, const uint64_t * __restrict M, const unsigned * __restrict bot, const unsigned nrb, const double * __restrict X, double * Y, unsigned lcmask, unsigned lrmask, unsigned clbits);
  
  void SSEspmv(const double * __restrict V, const unsigned short * __restrict M, const unsigned * __restrict bot, const unsigned nrb, const double * __restrict X, double * Y, unsigned lcmask, unsigned lrmask, unsigned clbits);
  
  void SSEspmv(const double * __restrict V, const unsigned char * __restrict M, const unsigned * __restrict bot, const unsigned nrb, const double * __restrict X, double * Y, unsigned lcmask, unsigned lrmask, unsigned clbits);
  
  template <class NT, class IT, unsigned TTDIM>
  class BmCsb
  {
  public:
  	BmCsb ():nz(0), m(0), n(0), nbc(0), nbr(0) {}	// default constructor (dummy)
  
  	BmCsb (const BmCsb<NT, IT, TTDIM> & rhs);			// copy constructor
  	~BmCsb();
  	BmCsb<NT,IT,TTDIM> & operator=(const BmCsb<NT,IT,TTDIM> & rhs);	// assignment operator
  	BmCsb (Csc<NT, IT> & csc, int workers);
  	
  	ofstream & PrintStats(ofstream & outfile) const;
  	IT colsize() const { return n;} 
  	IT rowsize() const { return m;} 
  	IT numregb() const { return nrb;}
  	bool isPar() const { return ispar; }
  
  private:
  	typedef typename int_least_helper<TTDIM>::least MTYPE;
  
  	void Init(int workers, IT forcelogbeta = 0);
  
  	void SubSpMV(IT * btop, IT bstart, IT bend, const NT * __restrict x, NT * __restrict suby, IT * __restrict sumscan) const;
  
  	void BMult(IT** chunks, IT start, IT end, const NT * __restrict x, NT * __restrict y, IT ysize, IT * __restrict sumscan) const;
  
  
  	void BlockPar(IT start, IT end, const NT * __restrict subx, NT * __restrict suby, 
  					IT rangebeg, IT rangeend, IT cutoff, IT * __restrict sumscan) const;
  
  	void SortBlocks(pair<IT, pair<IT,IT> > * pairarray, NT * val);
  
  	IT ** top ;	// pointers array (indexed by higher-order bits of the coordinate index), size = nbr*(nbc+1)
  	IT * bot;	// contains lower-order bits of the coordinate index, size nrb
  	MTYPE * masks;	// array of masks, size nrb
  	NT * num;	// contains numerical values, size nnz
  
  	bool ispar;
  	IT nz;		// # nonzeros 
  	IT m;		// # rows
  	IT n;		// # columns
  	IT blcrange;	// range indexed by one block
  
  	IT nbc;		// #{column blocks} = #{blocks in any block row}
  	IT nbr; 	// #{block rows)
  	IT nrb;		// #{register blocks}
  	
  	IT rowlowbits;	// # lower order bits for rows
  	IT rowhighbits;
  	IT highrowmask; // mask with the first log(m)/2 bits = 1 and the other bits = 0  
  	IT lowrowmask;
  
  	IT collowbits;	// # lower order bits for columns
  	IT colhighbits;
  	IT highcolmask; // mask with the first log(n)/2 bits = 1 and the other bits = 0  
  	IT lowcolmask;
  
  	MortonCompare<IT> mortoncmp;	// comparison operator w.r.t. the (inverted N)-morton layout
  
  	template <typename NU, typename IU, unsigned UUDIM>
  	friend void bmcsb_gespmv (const BmCsb<NU, IU, UUDIM> & A, const NU * x, NU * y);
  
  	template <class CSB>
  	friend float RowImbalance(const CSB & A);	// befriend any CSB instantiation	
  };
  
  
  #include "friends.h"
  #include "bmcsb.cpp"	
  #endif
  