英特尔® 线程挑战赛—数独
收藏
去年参加Intel的线程挑战赛，一直想写个博客，看见Intel在搞多核编程征文，就发上来了。如果您觉得文章对您有用，请投我一票，谢谢！投票地址：http://intel.csdn.net/multicoreblog/show.aspx
数独(Sudoku)是一种逻辑谜题游戏,解答方法是将数字放到栅格里,但同行同列或者同一个子块里不能有相同的数字。通常栅格一般取 9x9。这样的话，每一行每一列以及这 9 个 3x3 的不重叠子块的每一个，都会包括整数 1-9 的一个具体实例。除了 9x9 的栅格外，也有可能会用 16x16 或者 6x6 的，还有一些变化形式并不使用方形的子块。 问题：写一个分析程序判断 6x6 的数独(Sudoku)谜题是否有唯一解。正确的解应该将数字 1-6 不重复的放到每一行每一列。同样的，不相重叠的 2x3 子块(2行，3列)必须包含这六个整数的一个唯一的实例。例如： 2 3 1 | 4 5 6 输入文件描述：当程序开始执行时，将输入文件的名字送入，通常是以命令行参数的形式。输入文件会包含一定数量的行，每行有 36 个非空子符。每行代表了在一个可能的 6×6 数独(Sudoku)谜题的初始状态下，全部字符的逐行排列。迷题初始状态中的空格将由星号字符代替。文件的结尾表示输入的结束。 输出：输出要标准化。应该指出每一个输入谜题是有唯一解，无解还是有多解。 一次输入 3 个谜题的输入示例： *314*******1**356**621**3*******561* ***4****41*2*4321**6534*5*16****6*** 5******5******5******5******5******5 输出举例： 谜题 # 1 有唯一解 谜题 # 2 无解 谜题 # 3 有多解 计时：以时钟时间(Wall-clock time)为准(包括输入输出时间，即I/O时间) 数独由m x n 个 m x n 的子块，形成m x n x m x n = 36的栅格，(本题的m为2，n为3)，根据规则单元格中可填写的数字为 1到m x n，同时每个单元格Cell[i,j]必须遵守以下三个规则： 1. Cell[i,j]值在第i行是唯一的。 2. Cell[i,j]值在第j行是唯一的。 3. Cell[i,j]值在它所属的子块中是唯一的。 当每个Cell[i,j]都存在满足以上三个规则的值，我们便认为该数独有解。 换而言之如果两个单元格满足以下条件则这两个单元格存在关联关系,不能填入相同的值。 1. 单元格的横坐标相同。 2. 单元格的纵坐标相同。 3. 单元格所在子块相同。 与Cell[i,j]存在互斥关系的单元格一共有m x n x 3 - m – n - 1 个。 栅格的每个单元格用一个int值表示，初始填为(1 << XGRID_BSIZE) – 1【低XGRID_BSIZ个bit的值为1】，表示该单元格可以填的值为[1, XGRID_BSIZE]。 当一个单元格填入值t，就把它的关联单元格第t个bit置为0，表示该这些关联单元格不能填入t。再将这个填入值t的单元格高bit置为1【或上XFILLEDMASK】。例如： 单元格的值为 0x0003 表示该单元格填入1和2是合法的 单元格的值为 0x0005 表示该单元格填入1和3是合法的 单元格的值为 0x8008 表示该单元格一填入4。 加载迷题的时候，直接对有初始值的单元格进行填值，填值前对合法性进行检查，如果不合法直接返回无解。 迷题加载完成后，对所有未填值的单元格进行检测，找到一个可能性最少的【为1的bit位最少的单元格】单元格Cell[i,j]进行尝试填入合法的值【这样做可以降低堆栈的深度，减少分支的数量。】，如果在处理该Cell[i,j]的关联单元格时发现该关联单元格的值变为为0，表示该关联单元格没有了合法值，对Cell[i,j]尝试填入其他合法值。当所有单元格都填值完成，表示找到了该迷题的一组解。使用简单的递归调用即可完成求解。为追求效率，可以将递归展开。 首先定义一些宏、类型及结构。 // 填值失败 #define XFILL_INVALID 0x7FFFFFFF // 已填值标记 #define XFILLEDMASK 0x8000 // 每个小块的宽度 #define XGRID_BW 3 // 每个小块的高度 #define XGRID_BH 2 // 每个小块的大小 #define XGRID_BSIZE (XGRID_BW * XGRID_BH) // 整个网格的大小 #define XGRID_AREA (XGRID_BSIZE * XGRID_BSIZE) // 整个网格的大小SSE对齐 #define XGRID_AREASSE ((XGRID_AREA + 7) >> 3 << 3) // 每个单元格关联的单元格数量 #define XGRID_MUTUALITY (XGRID_BSIZE * 3 - XGRID_BW - XGRID_BH - 1) // 每个单元格关联的单元格数量SSE对齐 #define XGRID_MUTUALITYSSE ((XGRID_MUTUALITY + 7) >> 3 << 3) // 修改关联单元格的条件 #define XGRID_CHANGE_MUTUALITY_SIZE (XGRID_AREA * 30 / 100) // 堆栈的最大深度 #define XGRID_STACK_SIZE (XGRID_AREA - XGRID_CHANGE_MUTUALITY_SIZE) // 当数独的子块大于等于9的时候可以动态清除关联以提高效率 #if (XGRID_BSIZE >= 9) #define XGRID_MUTUALITYSSE_CLEAR #endif // 用于保存单元格的值 typedef int XCells[XGRID_AREASSE]; // 用于保存网格的数据结构 typedef struct tagXGrid { __declspec(align(16)) XCells xCells; // 单元格的值 #ifdef XGRID_MUTUALITYSSE_CLEAR // 单元格的关联单元格索引数组 __declspec(align(16)) uint8 pMutuality[XGRID_AREA][XGRID_MUTUALITYSSE]; // 单元格的关联单元格数量 __declspec(align(16)) uint8 nMutuality[XGRID_AREA]; #endif // 已经填充的单元格位置 __declspec(align(16)) int pFilledPos[XGRID_AREA]; // 已填充的单元格个数 int nFilled; // 运行模式和填充方案的数量 int *pRunMode,*pTotal; }XGrid; XGrid结构用于表示整个栅格。 初始化部分：
1，初始化PopCnt数组 g_xPopCnt，用于寻找填值方案最少的单元格，用查表法获取bit位为1的个数。 g_xPopCnt[0] = 0; g_xPopCnt[1] = 1; g_xPopCnt[2] = 1; g_xPopCnt[3] = 2; for(i = 0x04; i < 0x10; ++i) g_xPopCnt[i] = g_xPopCnt[i >> 2] + g_xPopCnt[i & 0x03]; for(i = 0x10; i < 0x100; ++i) g_xPopCnt[i] = g_xPopCnt[i >> 4] + g_xPopCnt[i & 0x0F]; #pragma omp parallel for for(i = 0x100; i < XFILLEDMASK; ++i) g_xPopCnt[i] = g_xPopCnt[i >> 8] + g_xPopCnt[i & 0xFF]; #pragma omp parallel for for(i = XFILLEDMASK; i < 0x10000; ++i) g_xPopCnt[i] = 0xFF; 2，初始化每个单元格的关联单元格数组。 // 计算每个单元格的关联单元格数组 for(i = 0; i < XGRID_AREA; ++i) { g_nMutuality[i] = 0; for(j = 0; j < XGRID_AREA; ++j) { if(i == j) continue; // 如果两个单元格满足以下条件则这两个单元格存在关联关系 // 1. 单元格的横坐标相同 // 2. 单元格的纵坐标相同 // 3. 单元格所在块的横坐标相同且单元格所在块的纵坐标相同[两个单元格在同一个块中] if( ( i % XGRID_BSIZE == j % XGRID_BSIZE ) || ( i / XGRID_BSIZE == j / XGRID_BSIZE ) || ( ( i / XGRID_BSIZE / XGRID_BH == j / XGRID_BSIZE / XGRID_BH ) && ( ( i % XGRID_BSIZE ) / XGRID_BW ) == ( ( j % XGRID_BSIZE ) / XGRID_BW ) ) ) { g_xMutuality[i][g_nMutuality[i]++] = j; } } for(j = XGRID_MUTUALITY; j < XGRID_MUTUALITYSSE; ++j) g_xMutuality[i][j] = XGRID_AREA; } 单元格填值部分：
1，填充单元格 前提为nPos位置的单元格置上nMask表示的值是合法的。首先将pCells[nPos]的值改为nMask + XFILLEDMASK，表示该单元格填上nMask表示的值，然后处理该单元格的关联单元格，去掉该单元格填充nMask的可能性，当某个关联单元格清除后仅剩下一种可能性则继续递归调用XFillCell填充掉该单元格，最后将填充的单元格位置保存到pGrid->pFilledPos中 __inline int XFillCell(XGrid* pGrid, XCells pCells, int nPos, int nMask) { int nFilled = pGrid->nFilled; for(;;) { #ifdef XGRID_MUTUALITYSSE_CLEAR uint8* const pMutuality = pGrid->pMutuality[nPos]; const uint8 nMutuality = pGrid->nMutuality[nPos]; #else uint8* const pMutuality = g_xMutuality[nPos]; const uint8 nMutuality = g_nMutuality[nPos]; #endif pCells[nPos] = nMask + XFILLEDMASK; pGrid->pFilledPos[nFilled++] = nPos; nPos = -1; #pragma unroll(16) for(int i = 0; i < nMutuality; ++i) { const uint8 nTmp = pMutuality[i]; // 去掉该单元格填充nMask的可能性 pCells[nTmp] &= ~nMask; // 如果单元格没有可选择的填值方案返回失败 if(pCells[nTmp] == 0) return XFILL_INVALID; // 判断该单元格是否只有一种填值可能,如果是将调用XFillCell给该单元格填值 if(g_xPopCnt[pCells[nTmp]] == 1) nPos = nTmp; } if(nPos < 0) break; nMask = pCells[nPos]; } pGrid->nFilled = nFilled; return pGrid->nFilled; } 迷题加载部分： 1，加载时填充单元格，与XFillCell不同，迷题加载的时候需要判断填值的合法性。 __inline int XFillCellLoad(XGrid* pGrid, XCells pCells, int nPos, int nMask) { int nNewPos = -1; #ifdef XGRID_MUTUALITYSSE_CLEAR uint8* const pMutuality = pGrid->pMutuality[nPos]; const uint8 nMutuality = pGrid->nMutuality[nPos]; #else uint8* const pMutuality = g_xMutuality[nPos]; const uint8 nMutuality = g_nMutuality[nPos]; #endif // 如果该单元格已经填充了nMask 直接返回 if(pCells[nPos] == nMask + XFILLEDMASK) return pGrid->nFilled; // 判断当前单元格能填nMask表示的值 if((pCells[nPos] & nMask) == 0) return XFILL_INVALID; // 为第nPos个单元格填上nMask表示的值 pCells[nPos] = nMask + XFILLEDMASK; for(int i = 0; i < nMutuality; ++i) { const uint8 nTmp = pMutuality[i]; // 去掉该单元格填充nMask的可能性 pCells[nTmp] &= ~nMask; // 如果单元格没有可选择的填值方案返回失败 if(pCells[nTmp] == 0) return XFILL_INVALID; // 判断该单元格是否只有一种填值可能,如果是将调用XFillCell给该单元格填值 if(g_xPopCnt[pCells[nTmp]] == 1) nNewPos = nTmp; } // 记录当前填充的单元格 pGrid->pFilledPos[pGrid->nFilled++] = nPos; if(nNewPos >= 0) { if(XFillCell(pGrid, pCells, nNewPos, pCells[nNewPos]) == XFILL_INVALID) { // 恢复已经填值的单元格 --pGrid->nFilled; return XFILL_INVALID; } } return pGrid->nFilled; } 2，加载迷题，将pInupt表示的迷题加载到pGrid中。 // 加载测试数据 __inline int XLoadGrid(XGrid* pGrid, const char* pInput) { int nPos = 0; int nMask = (1 << XGRID_BSIZE) - 1; #ifdef XGRID_MUTUALITYSSE_CLEAR XMEMCOPY(pGrid->pMutuality, g_xMutuality, sizeof(g_xMutuality)); XMEMCOPY(pGrid->nMutuality, g_nMutuality, sizeof(g_nMutuality)); #endif for(nPos = 0; nPos < XGRID_AREA; ++nPos) pGrid->xCells[nPos] = nMask; for(nPos = 0; nPos < XGRID_AREA; ++nPos) { if(pInput[nPos] >= '1' && pInput[nPos] <= '0' + XGRID_BSIZE) { // 填充nPos的值为1 << ( pInput[nPos] - '1' ) if(XFillCellLoad(pGrid, pGrid->xCells, nPos, 1 << ( pInput[nPos] - '1' )) == XFILL_INVALID) { return XFILL_INVALID; } } } #ifdef XGRID_MUTUALITYSSE_CLEAR // [XGRID_BSIZE>=9且需要查找所有情况时清除关联可加速] for(int i = 0; i < pGrid->nFilled; ++i) { XRemoveMutuality(pGrid, pGrid->pFilledPos[i]); } #endif return XGRID_AREA - pGrid->nFilled; }
迷题寻解部分：
1，查找可填的数字最少的单元格 __inline int XFindMinFeasibilityPos(XCells pCells) { int nMinCnt = 0xFF; int nMinPos = -1; // 在所有未填值的单元格中查找可填值方案最少的单元格 for(int i = 0; i < XGRID_AREA && nMinCnt > 1; ++i) { if(g_xPopCnt[pCells[i]] < nMinCnt) { nMinPos = i; nMinCnt = g_xPopCnt[pCells[i]]; } } return nMinPos; } 2，迷题求解，将递归展开。 __inline void XFillCellStack(XGrid* pGrid) { __declspec(align(16)) XCells xCellsStack[XGRID_AREA + 1]; __declspec(align(16)) int nMinPosStack[XGRID_AREA]; __declspec(align(16)) int nCurValStack[XGRID_AREA]; __declspec(align(16)) int nFilledStack[XGRID_AREA]; int nIndex = 0; XMEMCOPY(xCellsStack[0], pGrid->xCells, sizeof(xCellsStack[0])); // 查找当前填值可能性最少的单元格 nMinPosStack[nIndex] = XFindMinFeasibilityPos(xCellsStack[nIndex]); nCurValStack[nIndex] = xCellsStack[nIndex][nMinPosStack[nIndex]]; nFilledStack[nIndex] = pGrid->nFilled; while( (*pGrid->pRunMode & XRUN_FIND_FLAG) && nIndex >= 0 ) { nFilledStack[nIndex + 1] = XFILL_INVALID; // 尝试填充当前单元格 while(nFilledStack[nIndex + 1] == XFILL_INVALID && nCurValStack[nIndex]) { XMEMCOPY(xCellsStack[nIndex + 1], xCellsStack[nIndex], sizeof(xCellsStack[0])); unsigned int nMask = nCurValStack[nIndex] & -nCurValStack[nIndex]; nCurValStack[nIndex] -= nMask; nFilledStack[nIndex + 1] = XFillCell(pGrid, xCellsStack[nIndex + 1], nMinPosStack[nIndex], nMask); } if(nFilledStack[nIndex + 1] > XGRID_AREA) { // 没有合法的填充值 pGrid->nFilled = nFilledStack[--nIndex]; } else if(nFilledStack[nIndex + 1] != XGRID_AREA) { // 是合法的填充值, 还存在未填值的单元格 ++nIndex; // 查找当前填值可能性最少的单元格 nMinPosStack[nIndex] = XFindMinFeasibilityPos(xCellsStack[nIndex]); nCurValStack[nIndex] = xCellsStack[nIndex][nMinPosStack[nIndex]]; } else { // 所有单元格均已填值 XAddSolution(pGrid->pRunMode, pGrid->pTotal, xCellsStack[nIndex + 1]); pGrid->nFilled = nFilledStack[nIndex]; } } } 3，迷题求解。 // 对pInput表示的迷题求解 int XFindGrid(int nRunMode, const char* pInput) { int nTotal = 0; XGrid xGrid = {0}; xGrid.pRunMode = &nRunMode; xGrid.pTotal = &nTotal; // 加载网格数据 switch(XLoadGrid(&xGrid, pInput)) { case 0: nTotal = 1; break; case XFILL_INVALID: break; default: // 填充单元格 XFillCellStack(&xGrid); } return nTotal; } 并行优化： 串行模式 for(int i = 0; i < nCount; ++i) { pTotal[i] = XFindGrid(nRunMode, pInput + i * nLineSize); } 使用OpenMP优化 #pragma omp parallel for schedule(guided, 1) for(int i = 0; i < nCount; ++i) { pTotal[i] = XFindGrid(nRunMode, pInput + i * nLineSize); } 使用TBB优化 class CXSudokuTask { int m_nRunMode; int m_nLineSize; const char* m_pInput; int* m_pTotal; public: CXSudokuTask(int nRunMode, int nLineSize, int* pTotal, const char* pInput) : m_nRunMode(nRunMode), m_nLineSize(nLineSize), m_pTotal(pTotal),m_pInput(pInput){ } void operator () (const blocked_range<size_t> & r) const { for (size_t i = r.begin(); i != r.end(); ++ i) { m_pTotal[i] = XFindGrid(m_nRunMode, m_pInput + i * m_nLineSize); } } }; parallel_for(blocked_range<size_t>(0, nCount, max(1, nCount/32)), CXSudokuTask(nRunMode, nLineSize, pTotal, pInput)); Windows平台: 使用ICC编译.不使用PGO【使用PGO后发现效率降低】 Linux平台: 使用Icpc编译 1. 上传压缩包种的Src和Linux两个目录到服务器上. 2. 进入Linux目录执行make 3. 进入Bin目录 执行文件为XSoduku. 其他: 如果输出结果需要使用英文 请打开XSudoku.h第86行的XUSE_ENGLISH宏定义,默认使用汉语. 需要使用Intel编译器，TBB，请下载并正确安装。 测试平台： Time of this report: 4/9/2009, 14:16:02 Machine name: X-DELL Operating System: Windows XP Professional (5.1, Build 2600) Service Pack 3 (2600.xpsp_sp3_gdr.080814-1236) Language: Chinese (Regional Setting: Chinese) System Manufacturer: Dell Inc. System Model: Vostro 1400 BIOS: Phoenix ROM BIOS PLUS Version 1.10 A03 Processor: Intel(R) Core(TM)2 Duo CPU T5270 @ 1.40GHz (2 CPUs), ~550MHz Memory: 1014MB RAM Page File: 812MB used, 1628MB available Windows Dir: C:\WINDOWS DirectX Version: DirectX 9.0c (4.09.0000.0904) DX Setup Parameters: Not found DxDiag Version: 5.03.0001.0904 32bit Unicode 测试数据为100k个迷题。 使用串行模式测试 XSudoku.exe d100k.in os > ret.txt 计算时间:0.411790秒 输出时间:0.032550秒 总时间:0.444407秒 使用OpenMP并行模式测试 XSudoku.exe d100k.in opo > ret.txt 计算时间:0.209383秒 输出时间:0.034386秒 总时间:0.243842秒 使用TBB并行模式测试 XSudoku.exe d100k.in opt > ret.txt 计算时间:0.216393秒 输出时间:0.035294秒 总时间:0.251751秒
6 5 4 | 3 2 1
-------------
1 4 3 | 5 6 2
5 6 2 | 1 3 4
-------------
3 1 6 | 2 4 5
4 2 5 | 6 1 3
(第一行对应了上面已解决的那个谜题在初始状态下的字符排列)算法分析：
编译说明:
测试结果:
函数性能分析:
根据Vtune采集的数据分析:
目前时间开销最大的是XFillCell该函数功能是为某个单元格填值.Debug版本CPI为0.869.效率比较理想
其次是XFillCellLoad该函数功能是按输入的数据为某个单元格填值.Debug版本CPI为0.70.效率比较理想.
根据Thread Profiler采集数据的分析: 算法并行的时间达到94%,比较完美.
发表于 @ 2009年04月09日 22:20:00 | 评论( loading... ) | 举报| 收藏
新一篇:2009 英特尔® 线程挑战赛 第一题 基数排序
- 查看最新精华文章 请访问博客首页相关文章
-
- 发表评论
-
- 表 情：
-
- 评论内容：
- 用 户 名：