2009 英特尔® 线程挑战赛 第七题 矩阵乘法

2009 英特尔® 线程挑战赛—线段求交

邓辉   denghui0815@hotmail.com

问题描述

问题：编写一段多线程代码，使用 Strassen 算法将两个随机矩阵相乘。应用程序将生成两个矩阵：A(M,P) 和 B(P,N)，您需要使用 (1) 串行算法和 (2) Strassen 算法将它们相乘，得到乘积 C(M,N)。然后比较这两种算法的计算结果，确保使用 Strassen 算法得到的结果与使用串行算法得到的结果一致。

 

通过命令行输入应用程序参数值。输入的值将是 3 个整数（M、N 和 P），用来描述要使用的矩阵的大小。

 

代码限制：我们会预先提供一个用 C 语言编写的非常简单的串行计算应用程序。您应该在此源文件的基础上添加多线程代码，保留原来 main 函数、矩阵生成函数、串行乘法代码以及矩阵乘积结果比较函数的主体部分。如果要用其他语言实现，可以对这些代码做相应的更改。另外还可以更改内存分配和其他代码，以便实现 Strassen 计算的线程化。（建议您在编写代码时记下所做的更改以及进行这些更改的理由。）完成更改后，您提交的解决方案必须使用某种形式的 Strassen 算法计算出第二个矩阵乘积。

 

串行算法

原始的矩阵乘法为复杂度均为O(n^3），当矩阵足够大时，计算效率会非常低下，1969年提出的Strassen 矩阵乘法采用分治法，将A和B各分为四个子矩阵进行运算，降低了矩阵乘法的复杂度，算法如下：

S1 = A00 + A11   M1 = S1  x S6   C00 = M1 + M4 – M5 + M7

S2 = A10 + A11   M2 = S2  x B11  C01 = M3 + M5

S3 = A00 + A01   M3 = A00 x S7   C10 = M2 + M4

S4 = A10 – A00   M4 = A11 x S8   C11 = M1 + M3 – M2 + M6

S5 = A01 – A11   M5 = S3  x B11

S6 = B00 + B11   M6 = S4  x S9

S7 = B01 – B11   M7 = S5  x S10

S8 = B10 – B00

S9 = B00 + B01

S10 = B10 +B11

共采用7次乘法和18次加减法运算。

 

 

 

 

 

Winograd Strassen是一种变形的矩阵乘法，其原理相同，算法如下：

S1 = A10 + A11   M1 = S2  x S6   T1 = M1 + M2  C00 = M2 + M3

S2 = S1  - A00   M2 = A00 x B00  T2 = T1 + M4  C01 = T1 + M5 + M6

S3 = A00 - A10   M3 = A01 x B10                C10 = T2 - M7

S4 = A01 - S2    M4 = S3  x S7                 C11 = T2 + M5

S5 = B01 - B00   M5 = S1  x S5

S6 = B11 - S5    M6 = S4  x B11

S7 = B11 - B01   M7 = A11 x S8

S8 = S6  - B10

共采用7次乘法和15次加减法运算。

 

由于少3次加减法运算，所以Winograd效率会高于标准的Strassen 矩阵算法，但这两种算法的复杂度均为O(n^2.81）。

 

为提高效率，采用SSE2指令优化小块矩阵乘法运算C = A x B及C += A x B,见函数XSM_Mul_Small和XSM_MulAdd_Small，数据要求：矩阵的数据指针pA,pB,pC均为16字节对齐，n必须为16的倍数，p必须为2的倍数.

 

由于矩阵运算的数据量较大，cache命中率显得非常重要，所以对于中等尺寸的矩阵乘法采用先计算A00 x B00，A00 x B01，A10 x B01，A10 x B00再计算A01 x B10，A01 x B11，A11 x B11，A11 x B10的方法提高cache命中率。

 

Winograd Strassen矩阵乘法是一个递归算法，而且需要在递归函数内部使用内存保存中间结果S1—S8，M1—M7，T1，T2。如果在递归函数内存动态分配的方式，时间开销会比较大，所以采用在外部根据M，N，P的值计算整个递归需要的临时数据空间，动态分配后，传入递归函数内使用，降低时间开销。

并行算法

    通过使用TBB的Task，可以很方便的将递归算法并行化，在tbb::task* execute()中分裂7个子Task分别计算M1—M7即可实现Strassen矩阵乘法的并行优化。主要代码如下：

tbb::task* execute()

{

     if (m_nDepth == m_stUseOneMemoryDepth)

     {    // 进入串行算法内，降低内存分配和多线程调度开销

         XSM_Mul_Serial(m_nM, m_nN, m_nP, m_pA, m_nWidthA, m_pB, m_nWidthB, m_pC, m_nWidthC);

     }

     else

     {

         tbb::task_list list;

         ……

        

// M1 = S2  x S6

         list.push_back( *new( allocate_child() ) CXSMMulTask(…));

         // M2 = A00 x B00

         list.push_back( *new( allocate_child() ) CXSMMulTask(…));

         // M3 = A01 x B10

         list.push_back( *new( allocate_child() ) CXSMMulTask(…));

         // M4 = S3  x S7

         list.push_back( *new( allocate_child() ) CXSMMulTask(…));

         // M5 = S1  x S5

         list.push_back( *new( allocate_child() ) CXSMMulTask(…));

         // M6 = S4  x B11

         list.push_back( *new( allocate_child() ) CXSMMulTask(…));

         // M7 = A11 * S8

         list.push_back( *new( allocate_child() ) CXSMMulTask(…));

 

         set_ref_count(8);

 

         spawn_and_wait_for_all(list);

 

         ……

     }

 

     return NULL;

}

 

// Strassen矩阵乘法，并行版本

void XSM_Mul_Parallel(int m, int n, int p, XMType* pA, int nWidthA,

 XMType* pB, int nWidthB, XMType* pC, int nWidthC)

{

     tbb::task_scheduler_init init;

     // 获取深度

     int nDepth = XSM_Mul_GetDepth(m, n, p);

     // 控制Task数量应为7 * 7 = 49 次调用串行算法

     CXSMMulTask::SetUseOneMemoryDepth(max(0, nDepth - 2));

 

CXSMMulTask& xtask = *new(tbb::task::allocate_root()) CXSMMulTask(m, n, p, pA, nWidthA, pB, nWidthB, pC, nWidthC, nDepth);

     tbb::task::spawn_root_and_wait(xtask);

}

在Task处理前先判断递归深度是否小于m_stUseOneMemoryDepth，如果是可直接进入串行计算，减少Task的数量，降低TBB维护Task的开销，同时减少动态内存分配的时间消耗。由于测试机位16核，所以的m_stUseOneMemoryDepth是根据为nDepth – 2，最低层的Task数量为7 * 7 = 49个。

优化工具

Hotspots检测

使用Intel Amplifier的Hotspots检测功能查找二分查找算法的热点函数，结果如下：

检测结果显示热点函数在XSM_MulAdd_Small中，该函数已经采用SSE2指令优化。

 

Concurrency检测

使用Intel Amplifier的Concurrency检测功能查找可进行并行优化的代码，结果如下：

 

检测结果显示，main函数是完全串行的，占整个运行时间的比例较小，XSM_MulAdd_Small函数并行度较好，执行时间也比较长。

 

 

Locks and Waits检测

使用Intel Amplifier的Locks and Waits检测功能查找两种算法的锁和同步等待消耗，结果如下：

检测结果显示使用Task有不到100ms的锁和等待消耗。

其他优化

通过XSM_Mul_AdjustSize函数扩展M,N,P的值为MM，NN，PP以支持任意尺寸的矩阵乘法。调整后的MM，NN，PP在满足XSM_MulAdd_Small的数据要求的同时保证MM-M，NN-N，PP-P最小。

性能测试

操作系统：    32bit的测试在32位XP

CPU：         Intel(R) Core(TM)2  5670  @ 1.80GHz

    内存：       1G

   时间单位：    秒

 

M，N，P	原始串行乘法	串行Strassen	并行Strassen
512，512，512	1.88	0.09	0.08
1024，1024，1024	39.16	0.63	0.47
2048，2048，2048	324.44	4.45	3.27
4096，2048，2048	N	8.70	6.33

   

 

 

 

编译说明

Windows平台:

    使用VS2008和Intel Parallel Studio

1. 用VS2008打开本项目.

2. 选择Relase编译.

3. 进入Bin目录执行文件为XStrassen.exe.

Linux平台:

       使用ICC和TBB

1.     上传压缩包种的Src和Linux两个目录到服务器上.

2.     进入XStrassen/Linux目录 执行make

3.     进入XStrassen/Bin目录 执行文件为XStrassen.

其他：

       主办方请使用Win32平台Release版本测试，谢谢！

优化结论

通过Intel Parallel Studio和TBB，实现了Strassen矩阵乘法的并行优化，由于串行算法采用了单一内存分配，效率较高，而并行算法依然存在动态内存分配部分，导致效率降低，所以在双核PC上加速比较低，希望在16核的测试机上能有更好的结果。

致谢

感谢Clay Breshears所做的解答，感谢Mu，Pryce为本文章发表到ISN所做的工作，感谢Xia, JeffX P为我的解决方案进行了认真细致的翻译。