utility.h

#ifndef _UTILITY_H
#define _UTILITY_H

#define __int64 long long
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <climits>
#include <iostream>
#include <cmath>
#include <vector>
#include <mmintrin.h>  //  MMX
#include <xmmintrin.h> //  SSE
#include <emmintrin.h> //  SSE 2  
#include <pmmintrin.h> //  SSE 3

using namespace std;

#include <cilk/cilk_api.h>
#include <cilk/cilk.h>
#define SYNCHED __cilkrts_synched()
#define DETECT __cilkscreen_enable_checking()
#define ENDDETECT __cilkscreen_disable_checking()
#define WORKERS __cilkrts_get_nworkers()

#ifdef BWTEST
	#define UNROLL 100
#else
	#define UNROLL 1
#endif

#ifndef CILK_STUB
#ifdef __cplusplus
extern "C" {
#endif
/*
 * __cilkrts_synched
 *
 * Allows an application to determine if there are any outstanding
 * children at this instant. This function will examine the current
 * full frame to determine this.
 */

CILK_EXPORT __CILKRTS_NOTHROW
int __cilkrts_synched(void);

#ifdef __cplusplus
} // extern "C"
#endif
#else /* CILK_STUB */
/* Stubs for the api functions */
#define __cilkrts_synched() (1)
#endif /* CILK_STUB */

#ifdef STATS
	#include <cilk/reducer_opadd.h>
	cilk::reducer_opadd<__int64> blockparcalls;
	cilk::reducer_opadd<__int64> subspmvcalls;
	cilk::reducer_opadd<__int64> atomicflops;
#endif

void * address;
void * base;

using namespace std;

//  convert category to type 
  template< int Category > struct int_least_helper {}; // default is empty
  template<> struct int_least_helper<8> { typedef uint64_t least; };		// 8x8 blocks require 64-bit bitmasks
  template<> struct int_least_helper<4> { typedef unsigned short least; };	// 4x4 blocks require 16-bit bitmasks
  template<> struct int_least_helper<2> { typedef unsigned char least; };	// 2x2 blocks require 4-bit bitmasks, so we waste half of the array here

const uint64_t masktable64[64] = {0x8000000000000000, 0x4000000000000000, 0x2000000000000000, 0x1000000000000000,
				0x0800000000000000, 0x0400000000000000, 0x0200000000000000, 0x0100000000000000,
				0x0080000000000000, 0x0040000000000000, 0x0020000000000000, 0x0010000000000000,
				0x0008000000000000, 0x0004000000000000, 0x0002000000000000, 0x0001000000000000,
				0x0000800000000000, 0x0000400000000000, 0x0000200000000000, 0x0000100000000000,
				0x0000080000000000, 0x0000040000000000, 0x0000020000000000, 0x0000010000000000,
				0x0000008000000000, 0x0000004000000000, 0x0000002000000000, 0x0000001000000000,
				0x0000000800000000, 0x0000000400000000, 0x0000000200000000, 0x0000000100000000,
				0x0000000080000000, 0x0000000040000000, 0x0000000020000000, 0x0000000010000000,
				0x0000000008000000, 0x0000000004000000, 0x0000000002000000, 0x0000000001000000,
				0x0000000000800000, 0x0000000000400000, 0x0000000000200000, 0x0000000000100000,
				0x0000000000080000, 0x0000000000040000, 0x0000000000020000, 0x0000000000010000,
				0x0000000000008000, 0x0000000000004000, 0x0000000000002000, 0x0000000000001000,
				0x0000000000000800, 0x0000000000000400, 0x0000000000000200, 0x0000000000000100,
				0x0000000000000080, 0x0000000000000040, 0x0000000000000020, 0x0000000000000010,
				0x0000000000000008, 0x0000000000000004, 0x0000000000000002, 0x0000000000000001 };


const unsigned short masktable16[16] = {0x8000, 0x4000, 0x2000, 0x1000, 0x0800, 0x0400, 0x0200, 0x0100, 
					0x0080, 0x0040, 0x0020, 0x0010, 0x0008, 0x0004, 0x0002, 0x0001 };


const unsigned char masktable4[4] = { 0x08, 0x04, 0x02, 0x01 };	// mask for 2x2 register blocks


template <typename MTYPE>
MTYPE GetMaskTable(unsigned int index)
{
	return 0;
} 


template <>
uint64_t GetMaskTable<uint64_t>(unsigned int index)
{
	return masktable64[index];
} 

template <>
unsigned short GetMaskTable<unsigned short>(unsigned int index)
{
	return masktable16[index];
} 


template <>
unsigned char GetMaskTable<unsigned char>(unsigned int index)
{
	return masktable4[index];
} 

#define BALANCETH 2
#define RBDIM 8
#define RBSIZE (RBDIM*RBDIM)		// size of a register block (8x8 in this case)
#define SLACKNESS 8
#define KBYTE 1024
#define L2SIZE (256*KBYTE / RHSDIM)	// less than half of the L2 Cache (L2 should hold x & y at the same time) - scaled back by RHSDIM
#define CLSIZE 64			// cache line size

/* Tuning Parameters */
#define BREAKEVEN 4		// A block (or subblock) with less than (BREAKEVEN * dimension) nonzeros won't be parallelized
#define MINNNZTOPAR 128		// A block (or subblock) with less than MINNNZTOPAR nonzeros won't be parallelized
#define BREAKNRB (8/RBDIM)	// register blocked version of BREAKEVEN
#define MINNRBTOPAR (256/RBDIM)	// register blocked version of MINNNZPAR
#define LOGSERIAL 15
#define ROLLING 20

#define EPSILON 0.0001
#define REPEAT 10

// "absolute" difference macro that has no possibility of unsigned wrap
#define absdiff(x,y)   ( (x) > (y) ? (x-y) : (y-x))


unsigned rmasks[32] = { 0x00000001, 0x00000002, 0x00000004, 0x00000008,
			0x00000010, 0x00000020, 0x00000040, 0x00000080,  
			0x00000100, 0x00000200, 0x00000400, 0x00000800,  
			0x00001000, 0x00002000, 0x00004000, 0x00008000,  
			0x00010000, 0x00020000, 0x00040000, 0x00080000,  
			0x00100000, 0x00200000, 0x00400000, 0x00800000,  
			0x01000000, 0x02000000, 0x04000000, 0x08000000,  
			0x10000000, 0x20000000, 0x40000000, 0x80000000 };  


void popcountall(const uint64_t * __restrict M, unsigned * __restrict count, size_t size);
void popcountall(const unsigned short * __restrict M, unsigned * __restrict count, size_t size);
void popcountall(const unsigned char * __restrict M, unsigned * __restrict count, size_t size);


template <typename T>
void printhistogram(const T * scansum, size_t size, unsigned bins)
{
	ofstream outfile;
	outfile.open("hist.csv");
	vector<T> hist(bins);	// an STD-vector is zero initialized
	for(size_t i=0; i< size; ++i)
		hist[scansum[i]]++;

	outfile << "Fill_ratio" << "," << "count" << endl;
	for(size_t i=0; i< bins; ++i)
	{
		outfile << static_cast<float>(i) / bins  << "," << hist[i] << "\n";
	}
}

struct thread_data
{
	unsigned sum;
	unsigned * beg;
	unsigned * end;
};
	
unsigned int highestbitset(unsigned __int64 v);

template <typename MTYPE>
unsigned prescan(unsigned * a, MTYPE * const M, int n)	
{
	unsigned * end = a+n;
	unsigned * _a = a;	
	MTYPE * __restrict _M = M;
	unsigned int lgn;
	unsigned sum = 0;
	while ((lgn = highestbitset(n)) > LOGSERIAL)
	{
		unsigned _n = rmasks[lgn];	// _n: biggest power of two that is less than n
		int numthreads = SLACKNESS*WORKERS;
		thread_data * thdatas = new thread_data[numthreads];
		unsigned share = _n/numthreads;
		cilk_for(int t=0; t < numthreads; ++t)
		{
			popcountall(_M+t*share, _a+t*share, ((t+1)==numthreads)?(_n-t*share):share);
			thdatas[t].sum = 0;
			thdatas[t].beg = _a + t*share;
			thdatas[t].end = _a + (((t+1)==numthreads)?_n:((t+1)*share));
			thdatas[t].sum = accumulate(thdatas[t].beg, thdatas[t].end, thdatas[t].sum);
		}
		for(int t=0; t<numthreads; ++t)
		{
			unsigned temp = thdatas[t].sum;
			thdatas[t].sum = sum;
			sum += temp;
		}
		cilk_for(int tt=0; tt<numthreads; ++tt)
		{				
			unsigned * beg = thdatas[tt].beg;
			unsigned * end = thdatas[tt].end;	
			unsigned locsum = thdatas[tt].sum;

			while(beg != end)
			{
				unsigned temp = *beg;
				*beg++ = locsum;   // changing the value of (*beg) changes the corresponding aliased pointer _a as well
				locsum += temp; 
			}
		}
		_a += _n;	// move the pointer on a
		_M += _n;	// move the pointer on M
		n  &=  ~_n;	// clear the highest bit
		delete [] thdatas;
	}
	popcountall(_M, _a, end-(_a));
	while(_a != end)
	{
		unsigned temp = *_a;
		*_a = sum;
		sum += temp; 
		_a++;
	}
	return sum;
}

extern "C"
unsigned char *aligned_malloc( uint64_t size ) {
  unsigned char *ret_ptr = (unsigned char *)malloc( size + 16 );
  int temp = (unsigned long)ret_ptr & 0xF;
  int shift = 16 - temp;
  ret_ptr += shift;
  ret_ptr[ -1 ] = shift;
  return( ret_ptr );
}

extern "C"
void aligned_free( unsigned char *ptr ) {
  ptr -= ptr[ -1 ];
  free( ptr );
}


template <typename ITYPE>
ITYPE CumulativeSum (ITYPE * arr, ITYPE size)
{
    ITYPE prev;
    ITYPE tempnz = 0 ;
    for (ITYPE i = 0 ; i < size ; ++i)
    {
		prev = arr[i];
		arr[i] = tempnz;
		tempnz += prev ;	
    }
    return (tempnz) ;		    // return sum
}


template <typename T>
T machineEpsilon()
{	
	T machEps = 1.0;
 	do {
       		machEps /= static_cast<T>(2.0);
       		// If next epsilon yields 1, then break, because current
       		// epsilon is the machine epsilon.
    	}
    	while ((T)(static_cast<T>(1.0) + (machEps/static_cast<T>(2.0))) != 1.0);
 
    	return machEps;
}


template<typename _ForwardIter, typename T>
void iota(_ForwardIter __first, _ForwardIter __last, T __value)
{
	while (__first != __last)
     		*__first++ = __value++;
}
	
template<typename T, typename I>
T ** allocate2D(I m, I n)
{
	T ** array = new T*[m];
	for(I i = 0; i<m; ++i) 
		array[i] = new T[n]();
	return array;
}

template<typename T, typename I>
void deallocate2D(T ** array, I m)
{
	for(I i = 0; i<m; ++i) 
		delete [] array[i];
	delete [] array;
}


template < typename T >
struct absdiff : binary_function<T, T, T>
{
        T operator () ( T const &arg1, T const &arg2 ) const
        {
                using std::abs;
                return abs( arg1 - arg2 );
        }
};



template <int D>
void MultAdd(double & a, const double & b, const double & c)
{
	for(int i=0; i<D; i++)
	{
		a += b * c;
	}	
	
}

// bit interleave x and y, and return result
// only the lower order bits of x and y are assumed valid
template <typename ITYPE>
ITYPE BitInterleaveLow(ITYPE x, ITYPE y)
{
	ITYPE z = 0; // z gets the resulting Morton Number.
	int ite = sizeof(z) * CHAR_BIT / 2;

	for (int i = 0; i < ite; ++i) 
	{
		// bitwise shift operations have precedence over bitwise OR and AND
  		z |= (x & (1 << i)) << i | (y & (1 << i)) << (i + 1);
	}
	return z;
}

// bit interleave x and y, and return result z (which is twice in size)
template <typename ITYPE, typename OTYPE>
OTYPE BitInterleave(ITYPE x, ITYPE y)
{
	OTYPE z = 0; // z gets the resulting Morton Number.
	int ite = sizeof(x) * CHAR_BIT;

	for (int i = 0; i < ite; ++i) 
	{
		// bitwise shift operations have precedence over bitwise OR and AND
  		z |= (x & (1 << i)) << i | (y & (1 << i)) << (i + 1);
	}
	return z;
}

template <unsigned BASE>
inline unsigned IntPower(unsigned exponent)
{
	unsigned i = 1; 
	unsigned power = 1;

	while ( i <= exponent ) 
	{
		power *= BASE;
		i++;
	}
	return power;
}

template <>
inline unsigned IntPower<2>(unsigned exponent)
{	
	return rmasks[exponent];
}



// T should be uint32, uint64, int32 or int64; force concept requirement
template <typename T>
bool IsPower2(T x)
{
	return ( (x>0) && ((x & (x-1)) == 0));
}

unsigned int nextpoweroftwo(unsigned int v)
{
	// compute the next highest power of 2 of 32(or 64)-bit n
	// essentially does 1 << (lg(n - 1)+1).

	unsigned int n = v-1;

	// any "0" that is immediately right to a "1" becomes "1" (post: any zero has at least two "1"s to its left) 
	n |= n >> 1;

	// turn two more adjacent "0" to "1" (post: any zero has at least four "1"s to its left)
	n |= n >> 2;
	n |= n >> 4;	// post: any zero has at least 8 "1"s to its left
	n |= n >> 8;	// post: any zero has at least 16 "1"s to its left
	n |= n >> 16;	// post: any zero has at least 32 "1"s to its left

	return ++n;
}

// 64-bit version
// note: least significant bit is the "zeroth" bit
// pre: v > 0
unsigned int highestbitset(unsigned __int64 v)
{
	// b in binary is {10,1100, 11110000, 1111111100000000 ...}  
	const unsigned __int64 b[] = {0x2ULL, 0xCULL, 0xF0ULL, 0xFF00ULL, 0xFFFF0000ULL, 0xFFFFFFFF00000000ULL};
	const unsigned int S[] = {1, 2, 4, 8, 16, 32};
	int i;

	unsigned int r = 0; // result of log2(v) will go here
	for (i = 5; i >= 0; i--) 
	{
		if (v & b[i])	// highestbitset is on the left half (i.e. v > S[i] for sure)
		{
			v >>= S[i];
			r |= S[i];
		} 
	}
	return r;
}

__int64 highestbitset(__int64 v)
{
	if(v < 0)
	{
		cerr << "Indices can not be negative, aborting..." << endl;
		return -1;
	}
	else
	{
		unsigned __int64 uv = static_cast< unsigned __int64 >(v);
		unsigned __int64 ur = highestbitset(uv);
		return static_cast< __int64 > (ur);
	}
}
		
// 32-bit version 
// note: least significant bit is the "zeroth" bit
// pre: v > 0
unsigned int highestbitset(unsigned int v)
{
	// b in binary is {10,1100, 11110000, 1111111100000000 ...}  
	const unsigned int b[] = {0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000};
	const unsigned int S[] = {1, 2, 4, 8, 16};
	int i;

	unsigned int r = 0; 
	for (i = 4; i >= 0; i--) 
	{
		if (v & b[i])	// highestbitset is on the left half (i.e. v > S[i] for sure)
		{
			v >>= S[i];
			r |= S[i];
		} 
	}
	return r;
}

int highestbitset(int v)
{
	if(v < 0)
	{
		cerr << "Indices can not be negative, aborting..." << endl;
		return -1;
	}
	else
	{	
		unsigned int uv = static_cast< unsigned int> (v);
		unsigned int ur = highestbitset(uv);
		return static_cast< int > (ur); 
	}
}

/* This function will return n % d.
   d must be one of: 1, 2, 4, 8, 16, 32, … */
inline unsigned int getModulo(unsigned int n, unsigned int d)
{
	return ( n & (d-1) );
} 

// Same requirement (d=2^k) here as well
inline unsigned int getDivident(unsigned int n, unsigned int d)
{
	while((d = d >> 1))
		n = n >> 1;
	return n;
}

#endif