sym.cpp

#include "sym.h"
#include "utility.h"
#include <cassert>


// You must declare "static data members" inside the class body, 
// and you also need to define them outside the class body. 
// When the program refers to a static data member that is only declared but not defined, 
// the "unresolved external symbol" error occurs.
template <class T, class ITYPE>
int Sym<T, ITYPE>::p_fetch = 0;


template <class T, class ITYPE>
void Sym<T, ITYPE>::Init()
{
	ITYPE roundrowup = nextpoweroftwo(m);
	ITYPE roundcolup = nextpoweroftwo(n);
	ITYPE roundnnzup = nextpoweroftwo(nz);

	ITYPE rowbits = highestbitset(roundrowup);
	ITYPE colbits = highestbitset(roundcolup);
	ITYPE nnzbits = highestbitset(roundnnzup);

	if(!(rowbits > 1 && colbits > 1))
	{
		cerr << "Matrix too small for this library" << endl;
		return;
	}

	rowhighbits = static_cast<ITYPE>(rowbits/2)-1;	// # higher order bits for rows
	colhighbits = static_cast<ITYPE>(colbits/2)-1;	// # higher order bits for columns

	rowlowbits = rowbits - rowhighbits;		// rowlowbits = # lower order bits for rows
	collowbits = colbits - colhighbits;		// collowbits = # lower order bits for columns

	lowrowmask = IntPower<ITYPE>(2,rowlowbits) - 1;
	highrowmask = ((roundrowup - 1) ^ lowrowmask);
	lowcolmask = IntPower<ITYPE>(2,collowbits) - 1;
	highcolmask = ((roundcolup - 1) ^ lowcolmask);

	nbc = IntPower<ITYPE>(2,colhighbits);		// #{block columns} = #{blocks in any block row}
	nbr = IntPower<ITYPE>(2,rowhighbits);		// #{block rows)
	ntop = nbc * nbr;							// size of top array
}

// Constructing empty Sym objects (size = 0) are not allowed.
template <class T, class ITYPE>
Sym<T, ITYPE>::Sym (ITYPE size, ITYPE rows, ITYPE cols): nz(size),m(rows),n(cols)
{
	assert(nz != 0 && n != 0 && m != 0);
	Init();

	num = new T[nz];
	bot = new ITYPE[nz];
	top = new ITYPE* [nbr];	

	for(ITYPE i=0; i<nbr; ++i)
		top[i] = new ITYPE[nbc+1]; 
}


// copy constructor
template <class T, class ITYPE>
Sym<T, ITYPE>::Sym (const Sym<T,ITYPE> & rhs)
: nz(rhs.nz), m(rhs.m), n(rhs.n), ntop(rhs.ntop), nbr(rhs.nbr), nbc(rhs.nbc), rowhighbits(rhs.rowhighbits),rowlowbits(rhs.rowlowbits),
highrowmask(rhs.highrowmask), lowrowmask(rhs.lowrowmask), colhighbits(rhs.colhighbits), collowbits(rhs.collowbits),
highcolmask(rhs.highcolmask), lowcolmask(rhs.lowcolmask)
{
	if(nz > 0)
	{
		num = new T[nz];
		bot = new ITYPE[nz];

		for(ITYPE i=0; i< nz; ++i)	
			num[i]= rhs.num[i];
		for(ITYPE i=0; i< nz; ++i)	
			bot[i]= rhs.bot[i];
	}
	if ( ntop > 0)
	{
		top = new ITYPE* [nbr];

		for(ITYPE i=0; i<nbr; ++i)
			top[i] = new ITYPE[nbc+1]; 

		for(ITYPE i=0; i<nbr; ++i)
			for(ITYPE j=0; j <= nbc; ++j) 
				top[i][j] = rhs.top[i][j];
	}
}

template <class T, class ITYPE>
Sym<T, ITYPE> & Sym<T, ITYPE>::operator= (const Sym<T, ITYPE> & rhs)
{
	if(this != &rhs)		
	{
		if(nz > 0)	// if the existing object is not empty
		{
			// make it empty
			delete [] bot;
			delete [] num;
		}
		if(ntop > 0)
		{
			for(ITYPE i=0; i<nbr; ++i)
				delete [] top[i];
			delete [] top;
		}

		nz	= rhs.nz;
		n	= rhs.n;
		m   	= rhs.m;
		ntop	= rhs.ntop;
		nbr 	= rhs.nbr;
		nbc 	= rhs.nbc;

		rowhighbits = rhs.rowhighbits;
		rowlowbits	= rhs.rowlowbits;
		highrowmask = rhs.highrowmask;
		lowrowmask	= rhs.lowrowmask;

		colhighbits = rhs.colhighbits;
		collowbits = rhs.collowbits;
		highcolmask = rhs.highcolmask;
		lowcolmask	= rhs.lowcolmask;

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

			for(ITYPE i=0; i< nz; ++i)	
				num[i]= rhs.num[i];
			for(ITYPE i=0; i< nz; ++i)	
				bot[i]= rhs.bot[i];
		}
		if(ntop > 0)
		{
			top = new ITYPE* [nbr];

			for(ITYPE i=0; i<nbr; ++i)
				top[i] = new ITYPE[nbc+1]; 

			for(ITYPE i=0; i<nbr; ++i)
				for(ITYPE j=0; j <= nbc; ++j) 
					top[i][j] = rhs.top[i][j];
		}
	}
	return *this;
}

template <class T, class ITYPE>
Sym<T, ITYPE>::~Sym()
{
	if( nz > 0)
	{
		delete [] bot;
		delete [] num;
	}
	if ( ntop > 0)
	{
		for(ITYPE i=0; i<nbr; ++i)
			delete [] top[i];
		delete [] top;
	}
}


template <class T, class ITYPE>
Sym<T, ITYPE>::Sym (Csc<T, ITYPE> & csc):nz(csc.nz),m(csc.m),n(csc.n)
{
	typedef std::pair<ITYPE, ITYPE> ipair;
	typedef std::pair<ITYPE, ipair> mypair;

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

	num = new T[nz];
	bot = new ITYPE[nz];
	
	top = new ITYPE* [nbr];
	for(ITYPE i=0; i<nbr; ++i)
		top[i] = new ITYPE[nbc+1]; 

	mypair * pairarray = new mypair[nz];
	ITYPE k = 0;
	for(ITYPE j = 0; j < n; ++j)
	{
		for (ITYPE 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 
			ITYPE hindex = (((highrowmask &  csc.ir[i] ) >> rowlowbits) << colhighbits)
										| ((highcolmask & j) >> collowbits);	
			ITYPE lindex = ((lowrowmask &  csc.ir[i]) << collowbits) | (lowcolmask & 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

	// Now, elements within each block are sorted with respect to their lindex (i.e. there are in row-major order)
	// This is because the default comparison operator of pair<T1,T2> uses lexicographic comparison: 
	// within each block, hindex is the same, so they are sorted w.r.t ipair(lindex,i) lexicographically.

	ITYPE cnz = 0;

	for(ITYPE i = 0; i < nbr; ++i)
	{
		for(ITYPE j = 0; j < nbc; ++j)
		{
			top[i][j] = cnz;
			while(cnz < nz && pairarray[cnz].first == ((i*nbc)+j) )	// as long as we're in that block
			{
				bot[cnz] = (pairarray[cnz].second).first;
				num[cnz++] = csc.num[(pairarray[cnz].second).second];
			}
		}
		top[i][nbc] = cnz;
	}
		
	delete [] pairarray;
}

template <typename T, typename ITYPE>
void Sym<T, ITYPE>::BMult(ITYPE * btop, ITYPE bstart, ITYPE bend, const T * x, T * y, ITYPE ysize) const
{
	if((btop[bend] - btop[bstart] < BREAKEVEN * ysize) || (bend-bstart == 1))	
	{
		// not enough nonzeros to amortize new array formation
		SubSpMV(btop, bstart, bend, x, y);
		return;
	}
	ITYPE nnb = (bend+bstart)/2;	

	cilk_spawn BMult(btop, bstart, nnb, x, y, ysize);
	if(SYNCHED)
	{
		BMult(btop, nnb, bend, x, y, ysize);
	}
	else
	{
		T * temp = new T[ysize];
		for(ITYPE i=0; i<ysize; ++i)
			temp[i] = 0.0;

		BMult(btop, nnb, bend, x, temp, ysize);
		cilk_sync;
		
		for(ITYPE i=0; i<ysize; ++i)
			y[i] += temp[i];

		delete [] temp;
	}
}

// double* restrict a; --> No aliases for a[0], a[1], ...
// bstart: block start index
// bend: block end index
template <class T, class ITYPE>
inline void Sym<T, ITYPE>::SubSpMV(ITYPE * __restrict btop, ITYPE bstart, ITYPE bend, 
							const T * __restrict x, T * __restrict suby) const
{
	for (ITYPE j = bstart ; j < bend ; ++j)		// for all blocks inside that block row
	{
		// get higher order bits for column indices
		ITYPE chi = ((j << collowbits) & highcolmask);
		const T * subx = &x[chi];
 
#ifdef AMDPREOPT
		//BlockPrefetch(&bot[btop[j]], btop[j+1]-btop[j], sizeof(ITYPE));
		//BlockPrefetch(&num[btop[j]], btop[j+1]-btop[j], sizeof(T));
		// optimization: block prefetch
		
		ITYPE blocknz = btop[j+1] - btop[j];
		if(blocknz > 8 && blocknz < 4096)
		{
			// number of elements in one cache line
			const int botinc = CLSIZE/sizeof(ITYPE);	
			const int numinc = CLSIZE/sizeof(T);

			for (ITYPE i=0; i < blocknz; i+=botinc)				// prefetch once for every cache line
			{
				__builtin_prefetch (&bot[btop[j]+i], 0, 0);		// prefetch read-only, with no temporal locality
			}
			for (ITYPE i=0; i < blocknz; i+=numinc)
			{	
				__builtin_prefetch (&num[btop[j]+i], 0, 0);	
			}
		}
#endif

#ifdef SSEOPT
		__m128i m1 = _mm_set_epi32 (mu1, mu1, mu1, mu1);
		__m128i m2 = _mm_set_epi64x(mlong1, mlong2);	// pack two __int64 integers (for x86_64)
#endif

#ifdef PERFDEB
		ofstream stat("Perfdeb.txt", ios::app);
		stat << btop[j+1] - btop[j] << endl;
#endif

		for (ITYPE k = btop[j] ; k < btop[j+1] ; ++k)	// for all nonzeros within ith block (expected =~ nnz/n = c)
		{			
			ITYPE rli = ((bot[k] >> collowbits) & lowrowmask);
			ITYPE cli = (bot[k] & lowcolmask);

			suby [rli] += num[k] * subx [cli] ;
		}
	}
}

// This is a static member function does not have a "this" pointer
template <class T, class ITYPE>
void Sym<T, ITYPE>::BlockPrefetch(void * addr,  int total, int ssize)
{
	p_fetch = 0;
	int * __restrict a = (int*) addr;
	const int inc = CLSIZE/ssize;	// number of elements in one cache line

	// Grab every 64th address
	for(int i=0; i<total; i+=inc)
	{
		p_fetch += a[i];
	}
}



template <class T, class ITYPE>
void Sym<T, ITYPE>::Transpose()
{
	// when we jump to the next block in the same block-column, we move leaddim positions inside "top" array
	// leadim ~= sqrt(n) => number of blocks in each block-row
	ITYPE leaddim = lowcolmask+1;	
	Sym symT(nz, m, n);				// create empty transposed object

	ITYPE k = 0;
	ITYPE cnz = 0;

	for(ITYPE j = 0; j < leaddim; ++j)	// scan columns of top-level structure (~sqrt(n) iterations)
	{
		for(ITYPE i = j; i < ntop ; i += leaddim)	// iterates ~ sqrt(m) times within the block column
		{
			symT.top[k++] = cnz;
			cnz += top[i+1]-top[i];
		}
	}
	symT.top[k] = cnz;

	// Embarrassingly parallel sort of indices to get new bottom array
	// ITYPE nindex = (highmask &  csc.ir [i]) | ((highmask & bot) >> 4);
}