bmcsb.h

#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