csbsym.cpp

#include "csbsym.h"
#include <iterator>
#include "utility.h"

// Choose block size as big as possible given the following constraints
// 1) The bot array is addressible by IT
// 2) The parts of x & y vectors that a block touches fits into L2 cache [assuming a saxpy() operation]
// 3) There's enough parallel slackness for block rows (at least SLACKNESS * CILK_NPROC)
template <class NT, class IT>
void CsbSym<NT, IT>::Init(int workers, IT forcelogbeta)
{
	ispar = (workers > 1);
	IT roundup = nextpoweroftwo(n);

	// if indices are negative, highestbitset returns -1, 
	// but that will be caught by the sizereq below
	IT nbits = highestbitset(roundup);
	bool sizereq;
	if (ispar)
	{
		sizereq = (IntPower<2>(nbits) > SLACKNESS * workers);
	}
	else
	{
		sizereq = (nbits > 1);
	}
	if(!sizereq)
	{
		cerr << "Matrix too small for this library" << endl;
		return;
	}

	nlowbits = nbits-1;	
	IT inf = numeric_limits<IT>::max();
	IT maxbits = highestbitset(inf);

	nhighbits = nbits-nlowbits;	// # higher order bits for rows (has at least one bit)
	if(ispar)
	{
		while(IntPower<2>(nhighbits) < SLACKNESS * workers)
		{
			nhighbits++;
			nlowbits--;
		}
	}

	// calculate the space that suby and subx occupy in L2 cache
	IT yL2 = IntPower<2>(nlowbits) * sizeof(NT);
	while(yL2 > L2SIZE)
	{
		yL2 /= 2;
		nhighbits++;
		nlowbits--;
	}

	lowmask = IntPower<2>(nlowbits) - 1;
	if(forcelogbeta != 0)
	{
		IT candlowmask  = IntPower<2>(forcelogbeta) -1;
		cout << "Forcing beta to "<< (candlowmask+1) << " instead of the chosen " << (lowmask+1) << endl;
		cout << "Warning : No checks are performed on the beta you have forced, anything can happen !" << endl;
		lowmask = candlowmask;
		nlowbits = forcelogbeta;
		nhighbits = nbits-nlowbits; 
	}
	else 
	{
		double sqrtn = sqrt(static_cast<double>(n));
		IT logbeta = static_cast<IT>(ceil(log2(sqrtn))) + 2;
		if(nlowbits > logbeta)
		{
			nlowbits = logbeta;
			lowmask = IntPower<2>(logbeta) -1;
			nhighbits = nbits-nlowbits;
		}
		cout << "Beta chosen to be "<< (lowmask+1) << endl;
	}
	highmask = ((roundup - 1) ^ lowmask);
	
	IT blcdim = lowmask + 1;
        ncsb = static_cast<IT>(ceil(static_cast<double>(n) / static_cast<double>(blcdim)));
	
	blcrange = (lowmask+1) * (lowmask+1);	// range indexed by one block
	mortoncmp = MortCompSym<IT>(nlowbits, lowmask);
}


// copy constructor
template <class NT, class IT>
CsbSym<NT, IT>::CsbSym (const CsbSym<NT,IT> & rhs)
: nz(rhs.nz), n(rhs.n), blcrange(rhs.blcrange), ncsb(rhs.ncsb), nhighbits(rhs.nhighbits), nlowbits(rhs.nlowbits), 
highmask(rhs.highmask), lowmask(rhs.lowmask), mortoncmp(rhs.mortoncmp), ispar(rhs.ispar), diagonal(rhs.diagonal)
{
	if(nz > 0)  // nz > 0 iff nrb > 0
	{
		num = new NT[nz]();
		bot = new IT[nz];

		copy ( rhs.num, rhs.num+nz, num);
		copy ( rhs.bot, rhs.bot+nz, bot );
	}
	if ( ncsb > 0)
	{
		top = new IT* [ncsb];
		for(IT i=0; i<ncsb; ++i)
			top[i] = new IT[ncsb-i+1]; 
		for(IT i=0; i<ncsb; ++i)
			for(IT j=0; j <= (ncsb-i); ++j) 
				top[i][j] = rhs.top[i][j];
	}
}

template <class NT, class IT>
CsbSym<NT, IT> & CsbSym<NT, IT>::operator= (const CsbSym<NT, IT> & rhs)
{
	if(this != &rhs)		
	{
		if(nz > 0)	// if the existing object is not empty
		{
			// make it empty
			delete [] bot;
			delete [] num;
		}
		if(ncsb > 0)
		{
			for(IT i=0; i<ncsb; ++i)
				delete [] top[i];
			delete [] top;
		}
		ispar 	= rhs.ispar;
		nz	= rhs.nz;
		n	= rhs.n;
		ncsb 	= rhs.ncsb;
		blcrange = rhs.blcrange;
		mortoncmp = rhs.mortoncmp;
		
		nhighbits = rhs.nhighbits;
		nlowbits = rhs.nlowbits;
		highmask = rhs.highmask;
		lowmask	 = rhs.lowmask;
		diagonal = rhs.diagonal;	// copy the whole sparse vector

		if(nz > 0)	// if the copied object is not empty
		{
			num = new NT[nz]();
			bot = new IT[nz];

			copy ( rhs.num, rhs.num+nz, num);
			copy ( rhs.bot, rhs.bot+nz, bot );
		}
		if(ncsb > 0)
		{
			top = new IT* [ncsb];
			for(IT i=0; i<ncsb; ++i)
				top[i] = new IT[ncsb-i+1]; 
			for(IT i=0; i<ncsb; ++i)
				for(IT j=0; j <= (ncsb-i); ++j) 
					top[i][j] = rhs.top[i][j];
		}
	}
	return *this;
}

template <class NT, class IT>
CsbSym<NT, IT>::~CsbSym()
{
	if( nz > 0)
	{
		delete [] bot;
		delete [] num;
	}
	if ( ncsb > 0)
	{
		for(IT i=0; i<ncsb; ++i)
			delete [] top[i];
		delete [] top;
	}
}

template <class NT, class IT>
CsbSym<NT, IT>::CsbSym (Csc<NT, IT> & csc, int workers):nz(csc.nz), n(csc.n)
{
        typedef std::pair<IT, IT> ipair;
        typedef std::pair<IT, ipair> mypair;

        assert(nz != 0 && n != 0);
        Init(workers);

	top = new IT* [ncsb];
	for(IT i=0; i<ncsb; ++i)
		top[i] = new IT[ncsb-i+1]; 

        mypair * pairarray = new mypair[nz];
        IT k = 0;
        for(IT j = 0; j < n; ++j)
        {
                for (IT i = csc.jc [j] ; i < csc.jc[j+1] ; ++i) // scan the jth column
                {
                        // concatenate the higher/lower order half of both row (first) index and col (second) index bits 
                        IT hindex = (((highmask &  csc.ir[i] ) >> nlowbits) << nhighbits) | ((highmask & j) >> nlowbits);
                        IT lindex = ((lowmask &  csc.ir[i]) << nlowbits) | (lowmask & j) ;

                        // i => location of that nonzero in csc.ir and csc.num arrays
                        pairarray[k++] = mypair(hindex, ipair(lindex,i));
                }
        }
        sort(pairarray, pairarray+nz);  // sort according to hindex
        SortBlocks(pairarray, csc.num);
        delete [] pairarray;
}

template <class NT, class IT>
void CsbSym<NT, IT>::SortBlocks(pair<IT, pair<IT,IT> > * pairarray, NT * val)
{
        typedef pair<IT, pair<IT, IT> > mypair;
        IT cnz = 0;
	vector<IT> tempbot;
	vector<NT> tempnum;
        IT ldim = IntPower<2>(nhighbits);  // leading dimension (not always equal to ncsb)
        for(IT i = 0; i < ncsb; ++i)
        {
                for(IT j = 0; j < (ncsb-i); ++j)
                {
                        top[i][j] = tempbot.size();
                        IT prevcnz = cnz;
                        std::vector<mypair> blocknz;
                        while(cnz < nz && pairarray[cnz].first == ((i*ldim)+(j+i)) )        // as long as we're in this block
                        {
                                IT interlowbits = pairarray[cnz].second.first;
                                IT rlowbits = ((interlowbits >> nlowbits) & lowmask);
                                IT clowbits = (interlowbits & lowmask);
                                IT bikey = BitInterleaveLow(rlowbits, clowbits);

				if(j == 0 && rlowbits == clowbits)
				{
					diagonal.push_back(make_pair((i << nlowbits)+rlowbits, val[pairarray[cnz++].second.second]));
				}
				else
				{
                                	blocknz.push_back(mypair(bikey, pairarray[cnz++].second));
				}
                        }
                        // sort the block into bitinterleaved order
                        sort(blocknz.begin(), blocknz.end());
			typename vector<mypair>::iterator itr; 
	
			for( itr = blocknz.begin(); itr != blocknz.end(); ++itr) 
                        {
                                tempbot.push_back( itr->second.first );
                                tempnum.push_back( val[itr->second.second] );
                        }
                }
                top[i][ncsb-i] = tempbot.size();
        }

	assert (cnz == (tempbot.size() + diagonal.size()));
	nz = tempbot.size();	// update the number of off-diagonal nonzeros
	bot = new IT[nz];
	num = new NT[nz];

	copy(tempbot.begin(), tempbot.end(), bot);
	copy(tempnum.begin(), tempnum.end(), num);
	sort(diagonal.begin(), diagonal.end());
}

template<class NT, class IT>
void CsbSym<NT, IT>::DivideIterationSpace(IT * & lspace, IT * & rspace, IT & lsize, IT & rsize, IT size, IT d) const
{
	if(d == 1)
	{
		lsize = size-size/2;
		rsize = size/2;
		lspace = new IT[lsize];
		rspace = new IT[rsize]; 
		for(IT i=0; i<rsize; ++i)	// we alternate indices
		{
			lspace[i] = 2*i;
			rspace[i] = 2*i+1;
		}
		if(lsize > rsize)
		{
			lspace[lsize-1] = size-1;
		}
	}
	else 
	{
		IT chunks = size / (2*d);
		int rest = size - (2*d*chunks);	// rest is modulus 2d
		lsize = d*chunks;	// initial estimates
		rsize = d*chunks; 
		if(rest > d)	// first d goes to lsize, rest goes to rsize
		{
			rsize += (rest-d);
			lsize += d;
		}
		else	// all goes to lsize
		{
			lsize += rest;
		}
		lspace = new IT[lsize];
		rspace = new IT[rsize];
		int remrest = (int) rest;	// needs to be signed integer since we're looping it until negative
		if(d == 2)
		{
			for(IT i=0; i<chunks; ++i)	// we alternate indices
			{
				lspace[2*i+0] = 4*i+0;
				lspace[2*i+1] = 4*i+1;
				rspace[2*i+0] = 4*i+2;
				rspace[2*i+1] = 4*i+3; 
			}
			if(remrest-- > 0)	lspace[2*chunks+0] = 4*chunks+0;		
			if(remrest-- > 0)	lspace[2*chunks+1] = 4*chunks+1;
			if(remrest-- > 0)	rspace[2*chunks+0] = 4*chunks+2;
		}
		else if(d == 3)
		{
			for(IT i=0; i<chunks; ++i)
			{
				lspace[3*i+0] = 6*i+0;
				lspace[3*i+1] = 6*i+1;
				lspace[3*i+2] = 6*i+2;
				rspace[3*i+0] = 6*i+3;
				rspace[3*i+1] = 6*i+4;
				rspace[3*i+2] = 6*i+5;
			}
			if(remrest-- > 0)	lspace[3*chunks+0] = 6*chunks+0;
			if(remrest-- > 0)	lspace[3*chunks+1] = 6*chunks+1;
			if(remrest-- > 0)	lspace[3*chunks+2] = 6*chunks+2;
			if(remrest-- > 0)	rspace[3*chunks+0] = 6*chunks+3;
			if(remrest-- > 0)	rspace[3*chunks+1] = 6*chunks+4;
		}		
		else if(d == 4)
		{	
			
			for(IT i=0; i<chunks; ++i)
			{
				lspace[4*i+0] = 8*i+0;
				lspace[4*i+1] = 8*i+1;
				lspace[4*i+2] = 8*i+2;
				lspace[4*i+3] = 8*i+3;
				rspace[4*i+0] = 8*i+4;
				rspace[4*i+1] = 8*i+5;
				rspace[4*i+2] = 8*i+6;
				rspace[4*i+3] = 8*i+7;
			}
			if(remrest-- > 0)	lspace[4*chunks+0] = 8*chunks+0;
			if(remrest-- > 0)	lspace[4*chunks+1] = 8*chunks+1;
			if(remrest-- > 0)	lspace[4*chunks+2] = 8*chunks+2;
			if(remrest-- > 0)	lspace[4*chunks+3] = 8*chunks+3;
			if(remrest-- > 0)	rspace[4*chunks+0] = 8*chunks+4;
			if(remrest-- > 0)	rspace[4*chunks+1] = 8*chunks+5;
			if(remrest-- > 0)	rspace[4*chunks+2] = 8*chunks+6;
		}	
		else
		{
			cout << "Diagonal d = " << d << " is not yet supported" << endl;
		}	
	}
}

template<class NT, class IT>
void CsbSym<NT, IT>::MultAddAtomics(NT * __restrict y, const NT * __restrict x, const IT d) const
{
	cilk_for(IT i=0; i< ncsb-d; ++i)	// all blocks at the dth diagonal and beyond
	{
		IT rhi = (i << nlowbits);
                NT * __restrict suby = &y[rhi];
		const NT * __restrict subx_mirror = &x[rhi];
		
		cilk_for(IT j=d; j < (ncsb-i); ++j)
		{
			IT chi = ((j+i) << nlowbits);
			const NT * __restrict subx = &x[chi];
			NT * __restrict suby_mirror = &y[chi];

			IT * __restrict r_bot = bot;
        		NT * __restrict r_num = num;
			for(IT k=top[i][j]; k<top[i][j+1]; ++k)
			{
				IT ind = r_bot[k];
				NT val = r_num[k];

				IT rli = ((r_bot[k] >> nlowbits) & lowmask);
				IT cli = (r_bot[k] & lowmask);

				atomicallyIncrementDouble(&suby[rli], val * subx[cli]);
				atomicallyIncrementDouble(&suby_mirror[cli], val * subx_mirror[rli]);
#ifdef STATS
			        atomicflops += 2;
#endif
			}
		}
	}
}


template <class NT, class IT>
void CsbSym<NT, IT>::MultMainDiag(NT * __restrict y, const NT * __restrict x) const
{
	if(Imbalance(0) > 2 * BALANCETH)
	{
		cilk_for(IT i=0; i< ncsb; ++i)	// in main diagonal, j = i
		{
			IT hi = (i << nlowbits);
                	NT * __restrict suby = &y[hi];
			const NT * __restrict subx = &x[hi];
	
			if(i == (ncsb-1) && (n-hi) <= lowmask)	// last iteration and it's irregular (can't parallelize)
			{
				IT * __restrict r_bot = bot;
       		 		NT * __restrict r_num = num;
				for(IT k=top[i][0]; k<top[i][1]; ++k)
				{
					IT ind = r_bot[k];
					NT val = r_num[k];

					IT rli = ((ind >> nlowbits) & lowmask);
					IT cli = (ind & lowmask);
					
					suby[rli] += val * subx[cli];
					suby[cli] += val * subx[rli];	// symmetric update
				}
			}
			else
			{
				BlockTriPar(top[i][0], top[i][1], subx, suby, 0, blcrange, BREAKEVEN * (nlowbits+1));
			}
		}
	}
	else	// No need for block parallelization
	{
		cilk_for(IT i=0; i< ncsb; ++i)	// in main diagonal, j = i
		{
			IT hi = (i << nlowbits);
                	NT * __restrict suby = &y[hi];
			const NT * __restrict subx = &x[hi];
	
			IT * __restrict r_bot = bot;
       		 	NT * __restrict r_num = num;
			for(IT k=top[i][0]; k<top[i][1]; ++k)
			{
				IT ind = r_bot[k];
				NT val = r_num[k];
		
				IT rli = ((ind >> nlowbits) & lowmask);
				IT cli = (ind & lowmask);

				suby[rli] += val * subx[cli];
				suby[cli] += val * subx[rli];	// symmetric update
			}
		}
	}
	const IT diagsize = diagonal.size();
	cilk_for(IT i=0; i < diagsize; ++i)
	{
		y[diagonal[i].first] += diagonal[i].second * x[diagonal[i].first];	// process the diagonal
	}
}


// Multiply the dth block diagonal
// which is composed of blocks A[i][i+n]
template <class NT, class IT>
void CsbSym<NT, IT>::MultDiag(NT * __restrict y, const NT * __restrict x, const IT d) const
{
	if(d == 0)
	{
		MultMainDiag(y, x);
		return;
	}
	IT * lspace;
	IT * rspace;
	IT lsize, rsize;
	DivideIterationSpace(lspace, rspace, lsize, rsize, ncsb-d, d);
	IT lsum = 0;
	IT rsum = 0;
	for(IT k=0; k<lsize; ++k)
	{
		lsum += top[lspace[k]][d+1] - top[lspace[k]][d];
	}
	for(IT k=0; k<rsize; ++k)
	{
		rsum += top[rspace[k]][d+1] - top[rspace[k]][d];
	}
	float lave = lsum / lsize;
	float rave = rsum / rsize;

	cilk_for(IT i=0; i< lsize; ++i)	// in the dth diagonal, j = i+d
	{
		IT rhi = (lspace[i] << nlowbits) ;
		IT chi = ((lspace[i]+d) << nlowbits);
               	NT * __restrict suby = &y[rhi];
                NT * __restrict suby_mirror = &y[chi];
		const NT * __restrict subx = &x[chi];
		const NT * __restrict subx_mirror = &x[rhi];
		
		if((top[lspace[i]][d+1] - top[lspace[i]][d] > BALANCETH * lave)	// relative denser block
			&& (!(lspace[i] == (ncsb-d-1) && (n-chi) <= lowmask)))	// and parallelizable
		{
			BlockPar(top[lspace[i]][d], top[lspace[i]][d+1], subx, subx_mirror, suby, suby_mirror, 0, blcrange, BREAKEVEN * (nlowbits+1));
		}
		else
		{	
			IT * __restrict r_bot = bot;
	       	 	NT * __restrict r_num = num;
			for(IT k=top[lspace[i]][d]; k<top[lspace[i]][d+1]; ++k)
			{
				IT rli = ((r_bot[k] >> nlowbits) & lowmask);
				IT cli = (r_bot[k] & lowmask);
				suby[rli] += r_num[k] * subx[cli];
				suby_mirror[cli] += r_num[k] * subx_mirror[rli];	// symmetric update
			}
		}
	}
	
	cilk_for(IT j=0; j< rsize; ++j)
	{
		IT rhi = (rspace[j] << nlowbits) ;
		IT chi = ((rspace[j]+d) << nlowbits);
       		NT * __restrict suby = &y[rhi];
                NT * __restrict suby_mirror = &y[chi];
		const NT * __restrict subx = &x[chi];
		const NT * __restrict subx_mirror = &x[rhi];

		if((top[rspace[j]][d+1] - top[rspace[j]][d] > BALANCETH * rave)	// relative denser block
			&& (!(rspace[j] == (ncsb-d-1) && (n-chi) <= lowmask))) // and parallelizable
		{
			BlockPar(top[rspace[j]][d], top[rspace[j]][d+1], subx, subx_mirror, suby, suby_mirror, 0, blcrange, BREAKEVEN * (nlowbits+1));
		}
		else
		{
			IT * __restrict r_bot = bot;
        		NT * __restrict r_num = num;
			for(IT k=top[rspace[j]][d]; k<top[rspace[j]][d+1]; ++k)
			{
				IT rli = ((r_bot[k] >> nlowbits) & lowmask);
				IT cli = (r_bot[k] & lowmask);
				suby[rli] += r_num[k] * subx[cli];
				suby_mirror[cli] += r_num[k] * subx_mirror[rli];	// symmetric update
			}
		}
	}
	delete [] lspace;
	delete [] rspace;
}

// Block parallelization for upper triangular compressed sparse blocks
// start/end: element start/end positions (indices to the bot array)
// bot[start...end] always fall in the `same block
// PRECONDITION: rangeend-rangebeg is a power of two 
template <class NT, class IT>
void CsbSym<NT, IT>::BlockTriPar(IT start, IT end, const NT * __restrict subx, NT * __restrict suby, 
				IT rangebeg, IT rangeend, IT cutoff) const
{
	assert(IsPower2(rangeend-rangebeg));
	if(end - start < cutoff)
	{
		IT * __restrict r_bot = bot;
        	NT * __restrict r_num = num;
		for(IT k=start; k<end; ++k)
		{
			IT ind = r_bot[k];
			NT val = r_num[k];
			
			IT rli = ((ind >> nlowbits) & lowmask);
			IT cli = (ind & lowmask);

			suby[rli] += val * subx[cli];
			suby[cli] += val * subx[rli];	// symmetric update
		}
	}
	else
	{
		// Lower_bound is a version of binary search: it attempts to find the element value in an ordered range [first, last) 
		// Specifically, it returns the first position where value could be inserted without violating the ordering
		IT halfrange = (rangebeg+rangeend)/2;
		IT qrt1range = (rangebeg+halfrange)/2;
		IT qrt3range = (halfrange+rangeend)/2;

		IT * mid = std::lower_bound(&bot[start], &bot[end], halfrange, mortoncmp);	// divides in mid column
		IT * right = std::lower_bound(mid, &bot[end], qrt3range, mortoncmp);

		/* -------
		   | 0 2 |
		   | 1 3 |
		   ------- */
		// subtracting two pointers pointing to the same array gives you the # of elements separating them
		// In the symmetric case, quadrant "1" doesn't exist (size1 = 0)
		IT size0 = static_cast<IT> (mid - &bot[start]);
		IT size2 = static_cast<IT> (right - mid);
		IT size3 = static_cast<IT> (&bot[end] - right);

		IT ncutoff = std::max<IT>(cutoff/2, MINNNZTOPAR);
	    
		cilk_spawn BlockTriPar(start, start+size0, subx, suby, rangebeg, qrt1range, ncutoff);	// multiply subblock_0
		BlockTriPar(end-size3, end, subx, suby, qrt3range, rangeend, ncutoff);			// multiply subblock_3
		cilk_sync;

		BlockPar(start+size0, end-size3, subx, subx, suby, suby, halfrange, qrt3range, ncutoff); // multiply subblock_2
	}
}

// Parallelize the block itself
// start/end: element start/end positions (indices to the bot array)
// bot[start...end] always fall in the same block
// PRECONDITION: rangeend-rangebeg is a power of two 
// TODO: we rely on the particular implementation of lower_bound for correctness, which is dangerous !
//		 what if lhs (instead of rhs) parameter to the comparison object is the splitter?
template <class NT, class IT>
void CsbSym<NT, IT>::BlockPar(IT start, IT end, const NT * __restrict subx, const NT * __restrict subx_mirror, 
			NT * __restrict suby, NT * __restrict suby_mirror, IT rangebeg, IT rangeend, IT cutoff) const
{
	assert(IsPower2(rangeend-rangebeg));
	if(end - start < cutoff)
	{
		IT * __restrict r_bot = bot;
        	NT * __restrict r_num = num;
		for(IT k=start; k<end; ++k)
		{
			IT ind = r_bot[k];
			NT val = r_num[k];
					
			IT rli = ((ind >> nlowbits) & lowmask);
			IT cli = (ind & lowmask);

			suby[rli] += val * subx[cli];
			suby_mirror[cli] += val * subx_mirror[rli];	// symmetric update
		}
	}
	else
	{
		// Lower_bound is a version of binary search: it attempts to find the element value in an ordered range [first, last) 
		// Specifically, it returns the first position where value could be inserted without violating the ordering
		IT halfrange = (rangebeg+rangeend)/2;
		IT qrt1range = (rangebeg+halfrange)/2;
		IT qrt3range = (halfrange+rangeend)/2;

		IT * mid = std::lower_bound(&bot[start], &bot[end], halfrange, mortoncmp);
		IT * left = std::lower_bound(&bot[start], mid, qrt1range, mortoncmp);
		IT * right = std::lower_bound(mid, &bot[end], qrt3range, mortoncmp);

		/* -------
		   | 0 2 |
		   | 1 3 |
		   ------- */
		// subtracting two pointers pointing to the same array gives you the # of elements separating them
		// we're *sure* that the differences are 1) non-negative, 2) small enough to be indexed by an IT
		IT size0 = static_cast<IT> (left - &bot[start]);
		IT size1 = static_cast<IT> (mid - left);
		IT size2 = static_cast<IT> (right - mid);
		IT size3 = static_cast<IT> (&bot[end] - right);

		IT ncutoff = std::max<IT>(cutoff/2, MINNNZTOPAR);
	    
		// We only perform [0,3] in parallel and then [1,2] in parallel because the symmetric update causes races when 
		// performing [0,1] in parallel (as it would perform [0,2] in the fictitious lower triangular part)
		cilk_spawn BlockPar(start, start+size0, subx, subx_mirror, suby, suby_mirror, rangebeg, qrt1range, ncutoff);	// multiply subblock_0
		BlockPar(end-size3, end, subx, subx_mirror, suby, suby_mirror, qrt3range, rangeend, ncutoff);			// multiply subblock_3
		cilk_sync;

		cilk_spawn BlockPar(start+size0, start+size0+size1, subx, subx_mirror, suby, suby_mirror, qrt1range, halfrange, ncutoff);	// multiply subblock_1
		BlockPar(start+size0+size1, end-size3, subx, subx_mirror, suby, suby_mirror, halfrange, qrt3range, ncutoff);		// multiply subblock_2
		cilk_sync;
	}
}


// double* restrict a; --> No aliases for a[0], a[1], ...
// bstart/bend: block start/end index (to the top array)
template <class NT, class IT>
void CsbSym<NT, IT>::SeqSpMV(const NT * __restrict x, NT * __restrict y) const
{
	const IT diagsize = diagonal.size();
	for(IT i=0; i < diagsize; ++i)
	{
		y[diagonal[i].first] += diagonal[i].second * x[diagonal[i].first];	// process the diagonal
	}
	for (IT i = 0 ; i < ncsb ; ++i)    // for all block rows of A 
	{
		IT rhi = (i << nlowbits);
		NT * suby = &y[rhi];
		const NT * subx_mirror = &x[rhi];

		IT * __restrict r_bot = bot;
        	NT * __restrict r_num = num;
		for (IT j = 0 ; j < (ncsb-i) ; ++j)		// for all blocks inside that block row
		{
                	// get higher order bits for column indices
                	IT chi = ((j+i) << nlowbits);
                	const NT * __restrict subx = &x[chi];
			NT * __restrict suby_mirror = &y[chi];

			for(IT k=top[i][j]; k<top[i][j+1]; ++k)
			{
				IT rli = ((r_bot[k] >> nlowbits) & lowmask);
				IT cli = (r_bot[k] & lowmask);
				NT val = r_num[k];
				suby[rli] += val * subx[cli];
				suby_mirror[cli] += val * subx_mirror[rli];	// symmetric update
			}
		}
	}
}

// Imbalance in the dth block diagonal (the main diagonal is the 0th) 
template <class NT, class IT>
float CsbSym<NT, IT>::Imbalance(IT d) const
{
	if (ncsb <= d+1)
	{
		return 0.0; //pointless
	}

	IT size = ncsb-d-1;
	IT * sums = new IT[size];
        for(size_t i=0; i< size; ++i)
        {
		sums[i] = top[i][d+1] - top[i][d]; 
        }
	IT max = *max_element(sums, sums+size);
	IT mean = accumulate(sums, sums+size, 0.0) / size;
	delete [] sums;

        return static_cast<float>(max) / mean;
}


// Total number of nonzeros in the dth block diagonal (the main diagonal is the 0th) 
template <class NT, class IT>
IT CsbSym<NT, IT>::nzsum(IT d) const
{
        IT sum = 0; 
	for(size_t i=0; i< ncsb-d; ++i)
        {
		sum += (top[i][d+1] - top[i][d]); 
        }
        return sum;
}

// Print stats to an ofstream object
template <class NT, class IT>
ofstream & CsbSym<NT, IT>::PrintStats(ofstream & outfile) const 
{
	if(nz == 0)
	{
		outfile << "## Matrix Doesn't have any nonzeros" <<endl;
		return outfile;
	}
	const IT ntop = ncsb * ncsb; 	

	outfile << "## Average block is of dimensions "<< lowmask+1 << "-by-" << lowmask+1 << endl;
	outfile << "## Number of real blocks is "<< ntop << endl;
	outfile << "## Main (0th) block diagonal imbalance: " << Imbalance(0) << endl;
	outfile << "## 1st block diagonal imbalance: " << Imbalance(1) << endl;
	outfile << "## 2nd block diagonal imbalance: " << Imbalance(2) << endl;

	outfile << "## nnz ratios (block diagonal 0,1,2): " << static_cast<float>(nzsum(0)) / nz << ", " 
		<< static_cast<float>(nzsum(1)) / nz << ", " << static_cast<float>(nzsum(2)) / nz << endl; 
	outfile << "## atomics ratio: " << static_cast<float>(nz-nzsum(0)-nzsum(1)-nzsum(2))/nz << endl;
	
	std::vector<int> blocksizes;
	for(IT i=0; i<ncsb; ++i)
	{
		for(IT j=0; j < (ncsb-i); ++j) 
		{
			blocksizes.push_back(static_cast<int> (top[i][j+1]-top[i][j]));
		}
	}	
	sort(blocksizes.begin(), blocksizes.end());
	outfile<< "## Total number of nonzeros: " << 2*nz +diagonal.size()<< endl;
	outfile<< "## Total number of stored nonzeros: "<< nz+diagonal.size() << endl;
	outfile<< "## Size of diagonal: " << diagonal.size() << endl;

	outfile << "## Nonzero distribution (sorted) of blocks follows: \n" ;
	std::copy(blocksizes.begin(), blocksizes.end(), ostream_iterator<int>(outfile,"\n"));
	outfile << endl;
	return outfile;
}


template <class NT, class IT>
ofstream & CsbSym<NT, IT>::Dump(ofstream & outfile) const
{
	for(typename vector< pair<IT,NT> >::const_iterator itr = diagonal.begin(); itr != diagonal.end(); ++itr)
	{	
		outfile << itr->first << " " << itr->second << "\n";
	}
	for(IT i =0; i<ncsb; ++i)
	{
		for(IT j=0; j< (ncsb-i); ++j)
		{
			outfile << "Dumping A.top(" << i << "," << j << ")" << endl;
			for(IT k=top[i][j]; k< top[i][j+1]; ++k)
			{
				IT rli = ((bot[k] >> nlowbits) & lowmask);
				IT cli = bot[k] & lowmask;
				outfile << "A(" << rli << "," << cli << ")=" << num[k] << endl;
			}
		}
	}
	return outfile;
}