xshortestpub.cpp

最短路径的并行算法 采用TBB 需要intel编译器 速度很快

/*/////////////////////////////////////////
// 版权所有(C)  2000-2008 邓辉           //
// Email：      [e-mail]  //
// 说明：       Intel线程优化大赛参赛作品//
/////////////////////////////////////////*/

#include "XShortestPub.h"
#include "XFileMap.h"

// 读取一个名称
#define XREADNAME(pRead, pDst)													\
{																				\
	XMEMCOPY(pDst, pRead, XNODE_NAME_SIZE);										\
	pRead += XNODE_NAME_SIZE + 1;												\
}

// 读取一个数字
#define XREADNUM(pRead, pEnd, nRead)											\
{																				\
	while(*pTmp < '0' || *pTmp > '9') ++pTmp;									\
	nRead = *pRead++ - '0';														\
	while(*pRead >= '0' && *pRead <= '9')										\
	{																			\
		nRead = nRead * 10 + *pRead++ - '0';									\
	}																			\
	while(pRead < pEnd && *pRead != '\n') ++pRead;								\
	++pRead;																	\
}

void Error (char *s)
{
	fprintf(stderr, "ERROR: %s\n",s);
	exit(1);
}

void Warning (char *s)
{
	fprintf(stderr, "WARNING: %s\n",s);
}

// 加载测试数据
XGraph<uint32>* XDataFileRead(const char* szFile, uint32& nMode)
{
#ifdef XCODE_WIN
	// 获取系统信息
	SYSTEM_INFO SysInfo;
	GetSystemInfo(&SysInfo);
	const uint32 nGran = SysInfo.dwAllocationGranularity;
#else
	const uint32 nGran = 1 << 16;
#endif

	// 打开内存映射
	uint32 nFileSize = 0;

	XFILEMAPHANDLE hFileMap = XFileMapOpen(szFile, XFILEMAP_READONLY, nFileSize);

	if(hFileMap == XFILEMAPINVALID) Error("XFileMapOpen Error!");

	// 读取测试数据的大小
	const char *pInput = (const char*)XFileMapView(hFileMap, XFILEMAP_READONLY, 0, nFileSize);

	if(pInput == XMAPMEMINVALID) Error("XFileMapView Error!");

	const char * pTmp = pInput;
	const char * pEnd = pTmp + nFileSize;

	uint32 nNode = 0;

	XREADNUM(pTmp, pEnd, nNode);

	XGraph<uint32>* pGraph = XGraphCreate<uint32>(nNode);
	vector<string>&	xIDToNameTab = pGraph->xIDToNameTab;
	XNameToIDMap xNameToIDMap;
	XNameToIDMap::iterator	mapIter;
	XEdge<uint32>*	pEdgeCur = pGraph->pEdge + 1;
	XEdge<uint32>*	pEdgeEnd = pGraph->pEdge + pGraph->nNode * 2;
	XEdge<uint32>**	pHead = pGraph->pHead;
	uint32** ppDistance = pGraph->ppDistance;

	uint32 nIDA, nIDB, nDistance,nEdge = 0;
	char  szName[XNODE_NAME_SIZE + 2] = {0};

	// 自动分析使用Dijkstra算法还是Floyed算法
	int nGuessEdge = nFileSize / (XNODE_NAME_SIZE * 2 + 8);
	float fGuessScale = 0;
	if(nNode <= 500)
		fGuessScale = 2.0 * log(nNode) * (nNode + nGuessEdge) / nNode / nNode;
	else
		fGuessScale = (nNode + nGuessEdge) / nNode / 150.0;

	printf("Node %d GuessEdge %d Auto: %f\n", nNode, nGuessEdge, fGuessScale);

	if(nMode <= XRUN_AUTO_SERIAL)
	{
		if(fGuessScale > 1.0)
		{
			nMode +=  XRUN_FLOYED_PARALLEL - XRUN_AUTO_PARALLEL;
		}
		else
		{
			nMode +=  XRUN_DIJKSTRA_PARALLEL - XRUN_AUTO_PARALLEL;
		}
	}
	
	BOOL bFloyed = FALSE;

	if(nMode == XRUN_FLOYED_SERIAL || nMode == XRUN_FLOYED_PARALLEL) bFloyed = TRUE;

	nNode = 0;

	// 读取边的信息
	while(pTmp < pEnd) 
	{
		XREADNAME(pTmp, szName);

		// 通过名字查找ID
		mapIter = xNameToIDMap.find(szName);

		if(mapIter != xNameToIDMap.end()) 
		{
			nIDA = (*mapIter).second;
		}
		else
		{
			xNameToIDMap.insert(XNameToIDMap::value_type(szName, nNode));
			xIDToNameTab.push_back(szName);
			nIDA = nNode++;
		}

		XREADNAME(pTmp, szName);

		// 通过名字查找ID
		mapIter = xNameToIDMap.find(szName);

		if(mapIter != xNameToIDMap.end()) 
		{
			nIDB = (*mapIter).second;
		}
		else
		{
			xNameToIDMap.insert(XNameToIDMap::value_type(szName, nNode));
			xIDToNameTab.push_back(szName);
			nIDB = nNode++;
		}

		XREADNUM(pTmp, pEnd, nDistance);

		if(bFloyed)
		{
			// 初始化距离
			ppDistance[nIDA][nIDB] = nDistance;
			ppDistance[nIDB][nIDA] = nDistance;
		}
		else
		{
			// 如果边的空间不足,进行扩展
			if(pEdgeCur + 2 > pEdgeEnd)
			{
				XGraphExtEdgeMemory<uint32>(pGraph);
				pEdgeCur = pGraph->pEdge + 1;
				pEdgeEnd = pGraph->pEdge + pGraph->nNode * 2;
			}

			// 为两个节点各增加一条边
			pEdgeCur->nLink = nIDB;
			pEdgeCur->xKey  = nDistance;
			pEdgeCur->pNext = pHead[nIDA];
			pHead[nIDA] = pEdgeCur++;

			pEdgeCur->nLink = nIDA;
			pEdgeCur->xKey  = nDistance;
			pEdgeCur->pNext = pHead[nIDB];
			pHead[nIDB] = pEdgeCur++;
		}
		
		++nEdge;
	}

	pGraph->nEdge = nEdge;

	if( pGraph->nNode != xIDToNameTab.size() ) 
	{
		Warning("没有足够的节点!");
	}

	XFileMapUnView((void*)pInput, nFileSize);

	return pGraph;
}

uint32 *g_pDisMin;
int*	g_pDisMax;

// 比较函数 用于查找最小的距离
template <class Type>
BOOL XCompareMin(const Type& lhs, const Type& rhs)
{
	return g_pDisMin[lhs] < g_pDisMin[rhs];
}

// 比较函数 用于查找最大的距离
template <class Type>
BOOL XCompareMax(const Type& lhs, const Type& rhs)
{
	if(g_pDisMax[lhs] == XDISTANCE_INFINITY) return FALSE;
	if(g_pDisMax[rhs] == XDISTANCE_INFINITY) return TRUE;
	return g_pDisMax[lhs] > g_pDisMax[rhs];
}

// 输出计算结果
void XDataFileSave(const char* szFile, XGraph<uint32>* pGraph)
{
	int i,j,nOutPut,nCount,nNoPath = 0,nThread = 8;
	
	FILE* fp = fopen(szFile, "w");

	if(fp == NULL) Error("Open File Error!");

	g_pDisMax = (int*)pGraph->pDistance;
	g_pDisMin = pGraph->pDistance;


	// 分配索引空间
	int* pIndex = (int*)scalable_malloc(pGraph->nNode * (pGraph->nNode - 1) / 2 * sizeof(int));
	
	// 初始化索引空间 并统计无法联通的节点对
	for(i = 0,nCount = 0; i < pGraph->nNode; ++i) 
	{
		int nCur = i * pGraph->nAlign + i + 1;
		int nEnd = i * pGraph->nAlign + pGraph->nNode;
		for(j = nCur; j < nEnd; ++j)
		{
			pIndex[nCount++] = j;
			if(g_pDisMin[j] == XDISTANCE_INFINITY) ++nNoPath;
		}
	}

	nOutPut = min(20, nCount);
	
	int nSubSize = (nCount - nOutPut) / nThread;

	// 查找距离最小的20个节点对
	if(pGraph->nNode > 1000)
	{
		#pragma omp parallel for schedule(guided, 1)
		for(i = 0; i < nThread; ++i)
		{
			int nBeg = i * nSubSize;
			int nEnd = i == nThread -  1 ? nCount : nBeg + nSubSize;
			partial_sort(pIndex + nBeg, pIndex + nBeg + nOutPut, pIndex + nEnd, XCompareMin<int>);
		}

		#pragma omp parallel for
		for(i = 1; i < nThread; ++i)
		{
			int* pDst = pIndex + i * nOutPut;
			int* pSrc = pIndex + i * nSubSize;
			for(int j = 0; j < nOutPut; ++j)
			{
				int tmp = pDst[j]; pDst[j] = pSrc[j]; pSrc[j] = tmp;
			}
		}

		partial_sort(pIndex, pIndex + nOutPut, pIndex + nThread * nOutPut, XCompareMin<int>);
	}
	else
	{
		partial_sort(pIndex, pIndex + nOutPut, pIndex + nCount, XCompareMin<int>);
	}

	fprintf(fp, "Smallest:\n------------------------------\n");
	for(i = 0; i < nOutPut; ++i)
	{
		fprintf(fp, "%s %s %u\n", 
			pGraph->xIDToNameTab[pIndex[i] / pGraph->nAlign].c_str(),
			pGraph->xIDToNameTab[pIndex[i] % pGraph->nAlign].c_str(),
			pGraph->pDistance[pIndex[i]]);
	}

	// 查找距离最大的20个节点对
	if(pGraph->nNode > 1000)
	{
		#pragma omp parallel for schedule(guided, 1)
		for(i = 0; i < nThread; ++i)
		{
			int nBeg = i * nSubSize;
			int nEnd = i == nThread -  1 ? nCount : nBeg + nSubSize;
			partial_sort(pIndex + nBeg, pIndex + nBeg + nOutPut, pIndex + nEnd, XCompareMax<int>);
		}

		#pragma omp parallel for
		for(i = 1; i < nThread; ++i)
		{
			int* pDst = pIndex + i * nOutPut;
			int* pSrc = pIndex + i * nSubSize;
			for(int j = 0; j < nOutPut; ++j)
			{
				int tmp = pDst[j]; pDst[j] = pSrc[j]; pSrc[j] = tmp;
			}
		}

		partial_sort(pIndex, pIndex + nOutPut, pIndex + nThread * nOutPut, XCompareMax<int>);
	}
	else
	{
		partial_sort(pIndex, pIndex + nOutPut, pIndex + nCount, XCompareMax<int>);
	}

	fprintf(fp, "\nLargest:\n------------------------------\n");
	for(i = 0; i < nOutPut; ++i)
	{
		fprintf(fp, "%s %s %u\n", 
			pGraph->xIDToNameTab[pIndex[i] / pGraph->nAlign].c_str(),
			pGraph->xIDToNameTab[pIndex[i] % pGraph->nAlign].c_str(),
			pGraph->pDistance[pIndex[i]]);
	}

	fprintf(fp, "\nThere are %d node pairs that have no path.", nNoPath);
	
	scalable_free(pIndex);

	fclose(fp);
}


void XCreateTestData(const char* szFile, int nNode, double dDensity)
{
	XGraph<uint32>* pGraph = XGraphCreate<uint32>(nNode);
	vector<string>&	xIDToNameTab = pGraph->xIDToNameTab;
	uint32** ppDistance = pGraph->ppDistance;
	int i,j,k,nEdge = nNode * dDensity;

	for(i = 0; i < nEdge; ++i)
	{
		j = rand() % nNode;
		k = (rand() % (nNode -j)) + j;
		ppDistance[j][k] = ppDistance[k][j] =  (rand() % (XDISTANCE_INFINITY >> 2)) + 1;
	}

	for(i = 0; i < nNode; ++i)
	{
		for(j = i + 1; j < nNode ; ++j)
		{
			if(ppDistance[i][j] != XDISTANCE_INFINITY) break;
		}

		if(j == nNode)
		{
			k = (rand() % (nNode - i)) + i;
			ppDistance[i][k] = ppDistance[k][i] = (rand() % (XDISTANCE_INFINITY >> 2)) + 1;
			++nEdge;
		}
	}

	FILE* fp = fopen(szFile, "w");

	if(fp != NULL)
	{
		fprintf(fp, "%d\n", nNode);

		for(i = 0, nEdge = 0; i < nNode; ++i)
		{
			for(j = i + 1; j < nNode; ++j)
			{
				if(ppDistance[j][i] != XDISTANCE_INFINITY)
				{
					fprintf(fp, "%010d %010d %d\n", i, j, ppDistance[j][i]);
					++nEdge;
				}
			}
		}

		printf("Node %d Edge %d\n", nNode, nEdge);
	}

	fclose(fp);
}