快速电路仿真器（FastSPICE）中的高性能矩阵向量运算实现

　　今年10-11月份参加了EDA2020（第二届）集成电路EDA设计精英挑战赛，通过了初赛，并参加了总决赛，最后拿了一个三等奖，虽然成绩不是很好，但是想把自己做的分享一下，我所做的题目是概伦电子出的F题-快速电路仿真器（FastSPICE）中的高性能矩阵向量运算实现，下面我将给出自己的实现方案，仅供参考。

1、题目描述与分析

1.1、赛题叙述

首先先把题目写出来：

　　在晶体管级电路瞬态仿真过程中，仿真器需要根据电路连接关系并结合 KCL、KVL 定理建立微分方程，然后求解离散化的方程，中间每个步长输出电压和电流的仿真结果。在瞬态仿真中的每个步长都会涉及到大量的矩阵向量运算，例如电流输出的计算，需要计算很多条支路电流并求和。矩阵向量运算在整个仿真中占据了相当比例的仿真时间，尤其是针对存储器（memory）电路的快速电路仿真器（FastSPICE）。因此，加速矩阵向量运算可以直接提高电路的仿真速度。
　　电路方程的矩阵大小与电路的节点个数成正比，比如，大型的存储器电路（DRAM, Flash,SRAM 等等），其维度可达到 1000 万以上，如此高维度的矩阵以普通二维数组的方式存储会消耗巨大的内存。由于电路的每个节点通常只与少数节点有连接关系，因此电路方程矩阵通常是稀疏度很高的稀疏矩阵。稀疏矩阵特有的存储方式仅包含矩阵中的非零元素，可大大减少内存的消耗量，为存储高维度的矩阵提供可能。稀疏矩阵的矩阵向量运算与普通的矩阵向量运算略有不同，对整体运算速度有较明显的影响。希望参赛者能够开发一种针对电路方程的高效矩阵向量运算方案和实现，减少运算时间，从而提高电路仿真速度。
　　本赛题要求参赛队完成快速电路仿真器 (FastSPICE)中的相关支路电流计算任务，其中电路元件的电导将构成 N×N 阶稀疏矩阵，节点电压组成 N×1 阶矩阵，其计算过
程包含以下矩阵运算：
　　　　　　　　　　　　𝑅𝑛×1 = 𝐴𝑛×𝑛𝑆𝑛×1 （1）
　　　　　　　　　　　　𝐻𝑛×1 = 𝐿𝑛×1 − 𝐺𝑛×𝑛𝑆𝑛×1 (2)
　　　　　　　　　　　　𝐸𝑛×𝑛 = 𝐺𝑛×𝑛 + 𝐶𝑛×𝑛 (3)
　　其中稀疏矩阵的存储方式采用 Compressed Sparse Row(CSR)格式，由 3 个数组𝑟𝑜𝑤_𝑜𝑓𝑓𝑠𝑒𝑡、𝑐𝑜𝑙𝑢𝑚𝑛_𝑖𝑛𝑑𝑖𝑐𝑒𝑠和𝑣𝑎𝑙𝑢𝑒𝑠组成（如图 1 所示）。CSR 格式仅对稀疏矩阵中非零元进行存储，其中𝑐𝑜𝑙𝑢𝑚𝑛_𝑖𝑛𝑑𝑖𝑐𝑒𝑠和𝑣𝑎𝑙𝑢𝑒𝑠两个数组分别按行依次存储矩阵中非零元的列值和数值。𝑟𝑜𝑤_𝑜𝑓𝑓𝑠𝑒𝑡表示某一行的第一个非零元在𝑣𝑎𝑙𝑢𝑒𝑠数组里的偏移位置。𝑟𝑜𝑤_𝑜𝑓𝑓𝑠𝑒𝑡数组的长度为矩阵行数+1，𝑐𝑜𝑙𝑢𝑚𝑛_𝑖𝑛𝑑𝑖𝑐𝑒𝑠和𝑣𝑎𝑙𝑢𝑒𝑠两个数组的长度相同，为矩阵中非零元的个数。设稀疏矩阵第𝑖行（𝑖 = 0, 1, 2, … , 𝑛）的第一个零元的偏移位置为𝑟𝑜𝑤_𝑜𝑓𝑓𝑠𝑒𝑡[𝑖] = 𝑃𝑖，第𝑖 + 1行的第一个非零元的偏移位置为𝑟𝑜𝑤_𝑜𝑓𝑓𝑠𝑒𝑡[𝑖 + 1] = 𝑃𝑖+1 ， 则 第 𝑖 行 的 非 零 元 的 列 值 和 数 值 分 别 为𝑐𝑜𝑙𝑢𝑚𝑛_𝑖𝑛𝑑𝑖𝑐𝑒𝑠[𝑗], 𝑣𝑎𝑙𝑢𝑒𝑠[𝑗], 𝑗 ∈ [𝑃𝑖, 𝑃𝑖+1)，且第𝑖行的非零元的个数为𝑃𝑖+1 − 𝑃𝑖。


图1 CSR结构

 1.2 、赛题分析

　　本赛题的描述十分简单，要求完成一个矩阵-向量相乘的计算程序。选手们的程序的输入是一系列大小不一的矩阵和与矩阵尺寸匹配的向量。程序需要在保证准确性的情况下，用最短的时间把矩阵与向量相乘的结果写到内存指定的位置。

（1）.API会一次性给定多个矩阵（数量<1024），每个矩阵会有一个对应的向量。选手的任务就是计算每个矩阵和向量的乘积。

（2）.给定的每一个稀疏矩阵，其规模都不是特别大。最大的矩阵不会超过10000x10000。

　　个人思路：这道题目，需要实现double精度，且评分标准为60%看效率，30%看程序运行内存占有量，10%创新点，因此我们需要搞清楚CSR格式的结构，利用结构体指针实现double精度；而提高程序性能，可以从并行、指令集和缓存优化角度进行优化。

　　1）并行计算

　　并行计算大家都比较容易理解，所谓三个臭皮匠顶一个诸葛亮。一个大任务如果能够切分成若干个子任务，每人完成一个，那么总的时间一般都能够减少。并行计算的关键诀窍，在于【尽可能减少不同子任务之间的依赖性】。一个极端的例子，如果N个子任务之间完全不存在依赖性，各自独立同时进行。那么这N个子任务完成的总时间与一个子任务完成的时间相同。这种情况下，我们可以认为并行的加速比为N。这个加速比是相对不并行的情况而言的。另外一个极端的例子，如果第2个子任务所需要的数据需要等第1个子任务完成才能得到，第3个子任务所需要的数据需要等第2个子任务完成才能得到，依次类推。那么这种情况，哪怕再多的人来参与这个任务，大部分的人，在大部分时间都只是在等待。这样就无法从并行得到任何好处。一个良好的并行计算程序，需要去挖掘问题中的并行性，也就是设计并行计算方法，减少子任务之间的依赖性。那么这个题目中，你觉得怎样的并行是最优的呢？

　　2）指令集优化

　　我们平时写程序的时候，一般情况不会涉及到极致的性能要求。这种情况下一般不需要考虑在指令集层面进行优化。遇到像矩阵运算这一类十分底层的计算，就需要想尽一切办法，榨干CPU的最后一点算力来提高我们程序的性能。现代的CPU中，都有SIMD的指令集，例如x86体系的SSE/AVX指令集等。所谓SIMD（Single Instruction Multiple Data）指令集指的是同一条指令，同时作用在不同的数据上。这样可以实现指令集意义上的并行计算。需要注意的事情是，SIMD指令集并不能够自动的让你的程序实现加速。他需要你先将你的数据用规则的方式准备好。而稀疏矩阵的数据，本身是不规则的，也就是说非0元素有可能在任意位置。那么怎么把这些数据处理成规则的，使得SIMD指令集可以同时处理多个数据，需要各位选手自己去琢磨了。还有一点需要注意的是，将非规则的数据处理成规则的，这件事情本身也是需要额外消耗时间的。如果处理的不好，SIMD得到的收益，说不定都被这个额外消耗的时间抵消掉了。此外，除了SIMD之外，是否还有其他的指令集，可以用于稀疏矩阵运算的优化呢？这个是个开放的问题，需要大家自己调研，尝试。

　　3）缓存优化

　　除了并行计算和指令集优化这两个赛题中给了提示的方向。其实还有一个方向，各位同学有必要去考虑，那就是缓存优化。在计算机中，CPU与DRAM内存之间存取数据，速度是很慢的。因此现代的CPU一般都具有多级缓存。缓存的使用是CPU内部的事情，我们在程序中多数情况下不会直接操作缓存。然而，这并不意味着我们无法对缓存进行优化。如果能够了解CPU对缓存的存取策略，我们就可以设计良好的程序，使得缓存的命中率提高，这样就可以提高程序的性能。

1.3、实现思路

　　快速电路仿真器（FastSPICE）中的高性能矩阵向量运算实现算法主要分为两个实现算法，分别为 matrix_calc_taskA 和 matrix_calc_taskB。算法实现主要由电路元件电导构成的 N*N 阶稀疏矩阵和节点电压构成的 N*1 阶矩阵，计算过程主要由三个矩阵运算构成，分别为：
　　　　　　　　　　Rn*1=An*nSn*1(1)
　　　　　　　　　　Hn*1=Ln*1-Gn*nSn*1 (2)
　　　　　　　　　　En*n=Gn*n+Cn*n(3)
　　对于 matrix_calc_taskA 算法，主要实现针对于式（1）计算内容中的其中计算方法：
　　　　　　　　　　Idn*1=An*nSn*1(6)
　　而对于 matrix_calc_taskB 实现主要有式（1）以及式（2）、式（3）导出的计算任务:
　　　　　　　　　　IGn*1=Gn*nSn*1(4)
　　　　　　　　　　ICn*1=Cn*nSn*1(5)
　　　　　　　　　　Rn*1=Dn*1-Gn*nSn*1(7)
　　　　　　　　　　Hn*1=D’n*1-Cn*nSn*1(8)
　　　　　　　　　　An*n=Gn*n+aCn*n(9)
通过制作生成 matrix_calc_taskA、matrix_calc_taskB 动态链接库.so，给出电路方程的高效矩阵向量运算方案、具体实现及测试结果，实现了减少运算时间，提高电路仿真速度的目标。

2、技术路线和实现方式

2.1、matrix_calc_taskA 算法实现

　　外层循环指向所要计算的矩阵，通过矩阵阶数参数来划分线程数，调用 future 头文件实现多线程机制，使用特殊的 provider 进行数据的异步访问，从而实现线程间的通信。第一内层循环开始就逐渐添加到线程当中，每一个线程的开始索引地址为矩阵阶数/线程数*线程索引（第几个线程所要工作的起始地址），计算矩阵行数而不在循环当中计算索引数，第二个内层循环计算矩阵行索引，并计算索引行在偏移量数组当中的偏移数，第三层循环则通过偏移量循环与另一个目标向量计算得出结果，最终实现 Idn*1=An*nSn*1(6) 式，完成任务A，流程图如图2 所示。每一个线程都分摊一定工作量，并调用每一个线程的 wait 方法并发执行，每一个线程执 行 的 过 程 中 ， 在 计 算 (6)过 程 中 用 到 了 AVX 指 令 集 ， 通 过 使 用 AVX 指 令 集 中__mm256__loadu__pd、__mm256__mul__pd 等函数一步进行运算，加速实现 Idn*1=An*nSn*1，不需要等待普通循环中循环执行同样的操作，部分实现代码如下。 

 1     std::future<void> ft1 = async(std::launch::async,[&](){
 2     int count_maxtrix = 0; ////to store ---->size
 3 
 4 
 5     unsigned it = 0; //to store ---->rowArraySize
 6 
 7     int node = 0;    //to store ---->rowArray[]
 8     int j = 0;   //to store ---->rowOffset[]
 9 
10     // add here to enhance the speed
11     int my_rowArraysize = 0;
12     int my_rowOffset = 0;
13 
14     int maxtrix_begin  = size/16;
15     int maxtrix_end = size/16*2;
16 
17 
18     __m256d ID_m;
19     __m256d ID_mm;
20     __m256d NormalMatrix_m;
21     __m256d columnIndice_m;
22     __m256d mux_m;
23 
24     //double S_array[4];
25     int S_count;
26     double ID_array[4];
27 
28     bool ready = 0;
29 
30         for (count_maxtrix = maxtrix_begin; count_maxtrix < maxtrix_end; ++count_maxtrix)
31     {
32 
33         my_rowArraysize = (*(ptr+count_maxtrix))->rowArraySize;
34         //for (it = 0; it < (*(ptr+count_maxtrix))->rowArraySize; ++it)
35         for (it = 0; it < my_rowArraysize; ++it)
36         {
37             node = (*(ptr+count_maxtrix))->rowArray[it];
38 
39             //for (j = (*(ptr+count_maxtrix))->rowOffset[node]; j < (*(ptr+count_maxtrix))->rowOffset[node + 1]; ++j)
40                 my_rowOffset = (*(ptr+count_maxtrix))->rowOffset[node + 1];
41                
42 
43 
44                 j = (*(ptr+count_maxtrix))->rowOffset[node];
45 
46 
47 
48                 while((int)((my_rowOffset-j)/4) > 1)
49                 {   
50                     columnIndice_m =_mm256_set_pd((*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j+3]],
51                                         (*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j+2]],
52                                         (*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j+1]],
53                                         (*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j]] );
54                     //ID_m = _mm256_load_pd((*(ptr+count_maxtrix))->Id+j);
55                     NormalMatrix_m = _mm256_loadu_pd((*(ptr+count_maxtrix))->valueNormalMatrix+j);
56                     
57                 
58                 
59                     mux_m = _mm256_mul_pd(NormalMatrix_m, columnIndice_m);
60                 
61                     //ID_mm = _mm256_add_pd(ID_m, mux_m);
62                     
63                     _mm256_storeu_pd(ID_array, mux_m);
64 
65                     for(S_count = 0; S_count < 4; S_count++)
66                     {
67                         (*(ptr+count_maxtrix))->Id[node] += ID_array[S_count];
68                     }
69 
70                     ready = 1;
71 
72                     j += 4;     
73                 } 
74 
75                 if(ready == 1)
76                 {
77                     //for (; j < my_rowOffset; ++j)
78                     while(j < my_rowOffset)
79                     {
80                         (*(ptr+count_maxtrix))->Id[node] += (*(ptr+count_maxtrix))->valueNormalMatrix[j] * (*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j]];
81                         ++j;                    
82                     }
83                     ready = 0;
84                 }
85                 else
86                 {
87                     //for (; j < my_rowOffset; ++j)
88                     while(j < my_rowOffset)
89                     {
90                         (*(ptr+count_maxtrix))->Id[node] += (*(ptr+count_maxtrix))->valueNormalMatrix[j] * (*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j]];
91                         ++j;
92                     }
93                 }
94                 
95         }
96     }
97     });


图2 TaskA实现流程图 

2.2、matrix_calc_taskB 算法实现

　　与 matrix_calc_taskA 相类似，taskB 在调用 future 头文件以及 thread 头文件实现多线程机制，同样将该任务分为 16 线程数来并行计算。而在线程内，每个线程需要具体需要实现的计算为下所示：
　　　　　　　　　　IGn*1=Gn*nSn*1(4)
　　　　　　　　　　ICn*1=Cn*nSn*1(5)
　　　　　　　　　　Rn*1=Dn*1-Gn*nSn*1(7)
　　　　　　　　　　Hn*1=Dn*1-Cn*nSn*1(8)
　　　　　　　　　　An*n=Gn*n+aCn*n(9)
将指向矩阵的二级指针参数通过头文件的形式传递参数，并在每一个循环总设置后对应的引值，首先依据输入的数据数量进行 for 循环展开，取出需要计算矩阵网表的行号数组的个数
rowArraySize，展开之后取出所要计算的行号 rowArray 数组，之后通过对 TaskMatrixInfoB结构体的所需参数如 A、S、IC、IG 等进行调用，结合优化的及计算算法来实现 IGn*1、ICn*1、
Rn*1、Hn*1、An*n等式子的运算，之后再循环遍历数组即可，而再实现 IGn*1、ICn*1、Rn*1、Hn*1、An*n等式子的运算过程中，调用了 AVX 指令集，来进一步加速实现运算，流程图如图 3所示，部分TaskB实现代码如下：
  1     std::future<void> ft16 = async(std::launch::async,[&](){
  2 
  3         int count_maxtrix = 0; ////to store ---->size
  4 
  5         int it = 0;   //to store ---->rowArraySize
  6         int p = 0;    //to store ---->rowOffset
  7 
  8         int row = 0;  //to store ---->rowArray[x]
  9         int k = 0;    //to store ---->rowOffset*2
 10 
 11         int col = 0;    //to store ---->columnIndice
 12 
 13         double cond = 0; //to store ---->valueSpiceMatrix[k]
 14         double cap = 0; //to store ---->valueSpiceMatrix[k+1]
 15 
 16 
 17 
 18             /*
 19         function:Just for example ,for formula(7 8)
 20         Rnx1 = Dnx1 - Gnxn*Snx1
 21         Hnx1 = D'nx1 - Cnxn*Snx1
 22             */
 23         int kl = 0;     //to store ---->rowArray[]*2
 24         double current = 0;  //to store ---->D[kl]
 25         double charge = 0;   //to store ---->D[kl+1]
 26 
 27 
 28         // add here to enhance the speed
 29         int my_rowArraysize = 0;
 30         int my_rowOffset = 0;
 31         double my_s = 0;
 32 
 33         int maxtrix_begin  = size/16*7;
 34         int maxtrix_end = size/16*8;
 35 
 36 
 37         __m256d cond_m;
 38         __m256d cap_m;
 39         __m256d IG_Curr_m;
 40         __m256d IC_Char_m;
 41         __m256d A_m;
 42 
 43         double IG_Curr_array[4];
 44         double IC_Char_array[4];
 45         double A_m_array[4];
 46         int u_count;
 47         bool ready = 0;
 48 
 49         for (count_maxtrix = maxtrix_begin; count_maxtrix < maxtrix_end; ++count_maxtrix)
 50         {
 51 
 52                     //for (it = 0; it < (*(ptr+count_maxtrix))->rowArraySize; ++it)
 53             my_rowArraysize = (*(ptr+count_maxtrix))->rowArraySize;
 54             for (it = 0; it < my_rowArraysize; ++it)
 55             {
 56                 //…………………………………………formula(4 5)、(7 8)、(9) Share a for Loop………………………………………………………………
 57                 /*
 58                     Function:Just for example ,for formula(4 5)
 59                     IGnx1 = Gnxn*Snx1
 60                     ICnx1 = Cnxn*Snx1
 61                         */
 62                 row = (*(ptr+count_maxtrix))->rowArray[it];   //share
 63 
 64 
 65 
 66 
 67                 /*
 68                     function:Just for example ,for formula(7 8)
 69                     Rnx1 = Dnx1 - Gnxn*Snx1
 70                     Hnx1 = D'nx1 - Cnxn*Snx1
 71                     */
 72                 kl = row * 2;
 73 
 74                 current = (*(ptr+count_maxtrix))->D[kl];
 75                 charge = (*(ptr+count_maxtrix))->D[kl + 1];
 76 
 77 
 78 
 79                 //for (p = (*(ptr+count_maxtrix))->rowOffset[row]; p < (*(ptr+count_maxtrix))->rowOffset[row + 1]; ++p)
 80                 my_rowOffset = (*(ptr+count_maxtrix))->rowOffset[row + 1];
 81 
 82 
 83                 p = (*(ptr+count_maxtrix))->rowOffset[row];
 84                 //for (p = (*(ptr+count_maxtrix))->rowOffset[row]; p < my_rowOffset; ++p)
 85                 while((int)((my_rowOffset-p)/4) > 1)
 86                 {
 87                     //col = (*(ptr+count_maxtrix))->columnIndice[p];
 88                     k = p * 2;
 89                     //cond = (*(ptr+count_maxtrix))->valueSpiceMatrix[k];
 90                     //cap = (*(ptr+count_maxtrix))->valueSpiceMatrix[k + 1];
 91 
 92 
 93                     
 94                     /*
 95                     Function:Just for example ,for formula(4 5)
 96                     IGnx1 = Gnxn*Snx1
 97                     ICnx1 = Cnxn*Snx1
 98                         */
 99 
100                     //my_s = (*(ptr+count_maxtrix))->S[col];
101 
102                     //_mm256_blend_pd (__m256d a, __m256d b, const int imm8);
103 
104                     cond_m = _mm256_set_pd((*(ptr+count_maxtrix))->valueSpiceMatrix[k+6],
105                                             (*(ptr+count_maxtrix))->valueSpiceMatrix[k+4],
106                                             (*(ptr+count_maxtrix))->valueSpiceMatrix[k+2],
107                                             (*(ptr+count_maxtrix))->valueSpiceMatrix[k]);
108 
109                     
110                     cap_m = _mm256_set_pd((*(ptr+count_maxtrix))->valueSpiceMatrix[k + 7],
111                                     (*(ptr+count_maxtrix))->valueSpiceMatrix[k + 5],
112                                     (*(ptr+count_maxtrix))->valueSpiceMatrix[k + 3],
113                                     (*(ptr+count_maxtrix))->valueSpiceMatrix[k + 1]);
114                     /*
115                     alpha_m = _mm256_set_pd((*(ptr+count_maxtrix))->alpha,
116                                         (*(ptr+count_maxtrix))->alpha,
117                                         (*(ptr+count_maxtrix))->alpha,
118                                         (*(ptr+count_maxtrix))->alpha);
119 
120                     my_s_m =  _mm256_set_pd((*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p+3]]
121                                 ,(*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p+2]]
122                                 ,(*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p+1]]
123                                 ,(*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p]]); */
124 
125                     IG_Curr_m = _mm256_mul_pd(cond_m
126                             ,_mm256_set_pd((*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p+3]]
127                                 ,(*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p+2]]
128                                 ,(*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p+1]]
129                                 ,(*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p]]));
130                     IC_Char_m = _mm256_mul_pd(cap_m
131                                     ,_mm256_set_pd((*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p+3]]
132                                 ,(*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p+2]]
133                                 ,(*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p+1]]
134                                 ,(*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[p]]));
135 
136                     A_m = _mm256_add_pd(cond_m, _mm256_mul_pd(_mm256_set_pd((*(ptr+count_maxtrix))->alpha,
137                                         (*(ptr+count_maxtrix))->alpha,
138                                         (*(ptr+count_maxtrix))->alpha,
139                                         (*(ptr+count_maxtrix))->alpha)
140                                         ,cap_m));
141 
142         
143                     _mm256_storeu_pd(IG_Curr_array, IG_Curr_m);
144                     _mm256_storeu_pd(IC_Char_array, IC_Char_m);
145                     _mm256_storeu_pd(A_m_array, A_m);
146         
147 
148                     for(u_count = 0; u_count < 4; u_count++)
149                     {
150                             //(*(ptr+count_maxtrix))->Id[node] += ID_array[S_count];
151 
152                             (*(ptr+count_maxtrix))->IG[row] += IG_Curr_array[u_count];
153                             (*(ptr+count_maxtrix))->IC[row] += IC_Char_array[u_count];
154 
155                             /*                 
156                             current -= (*(ptr+count_maxtrix))->valueSpiceMatrix[k] * (*(ptr+count_maxtrix))->S[col];
157                             charge -= (*(ptr+count_maxtrix))->valueSpiceMatrix[k + 1] * (*(ptr+count_maxtrix))->S[col];
158                             */
159                             current -= IG_Curr_array[u_count];
160                             charge -= IC_Char_array[u_count];
161             
162 
163                             /*
164                             function:Just for example ,for formula(9)
165                             Anxn = Gnxn + alpha*Cnxn
166                                 */
167                             (*(ptr+count_maxtrix))->A[p+u_count] = A_m_array[u_count];
168                     }
169 
170 
171                     ready = 1;
172                     p += 4;
173 
174                 }
175 
176                 if(ready == 1)
177                 {
178                     //for (; j < my_rowOffset; ++j)
179                     while(p < my_rowOffset)
180                     {
181                         col = (*(ptr+count_maxtrix))->columnIndice[p];
182                         k = p * 2;
183                         cond = (*(ptr+count_maxtrix))->valueSpiceMatrix[k];
184                         cap = (*(ptr+count_maxtrix))->valueSpiceMatrix[k + 1];
185 
186 
187 
188                         /*
189                         Function:Just for example ,for formula(4 5)
190                         IGnx1 = Gnxn*Snx1
191                         ICnx1 = Cnxn*Snx1
192                             */
193 
194                         my_s = (*(ptr+count_maxtrix))->S[col];
195                         (*(ptr+count_maxtrix))->IG[row] += cond * my_s;
196                         (*(ptr+count_maxtrix))->IC[row] += cap * my_s;
197 
198 
199 
200 
201         /*                 current -= (*(ptr+count_maxtrix))->valueSpiceMatrix[k] * (*(ptr+count_maxtrix))->S[col];
202                         charge -= (*(ptr+count_maxtrix))->valueSpiceMatrix[k + 1] * (*(ptr+count_maxtrix))->S[col];
203         */
204                         current -= cond * my_s;
205                         charge -= cap * my_s;
206         
207 
208                         /*
209                         function:Just for example ,for formula(9)
210                         Anxn = Gnxn + alpha*Cnxn
211                             */
212                         (*(ptr+count_maxtrix))->A[p] = cond + (*(ptr+count_maxtrix))->alpha * cap;
213                         ++p;                    
214                     }
215                     ready = 0;
216                 }
217                 else
218                 {
219                     //for (; j < my_rowOffset; ++j)
220                     while(p < my_rowOffset)
221                     {
222                         col = (*(ptr+count_maxtrix))->columnIndice[p];
223                         k = p * 2;
224                         cond = (*(ptr+count_maxtrix))->valueSpiceMatrix[k];
225                         cap = (*(ptr+count_maxtrix))->valueSpiceMatrix[k + 1];
226 
227 
228 
229                         /*
230                         Function:Just for example ,for formula(4 5)
231                         IGnx1 = Gnxn*Snx1
232                         ICnx1 = Cnxn*Snx1
233                             */
234 
235                         my_s = (*(ptr+count_maxtrix))->S[col];
236                         (*(ptr+count_maxtrix))->IG[row] += cond * my_s;
237                         (*(ptr+count_maxtrix))->IC[row] += cap * my_s;
238 
239 
240 
241 
242         /*                 current -= (*(ptr+count_maxtrix))->valueSpiceMatrix[k] * (*(ptr+count_maxtrix))->S[col];
243                         charge -= (*(ptr+count_maxtrix))->valueSpiceMatrix[k + 1] * (*(ptr+count_maxtrix))->S[col];
244         */
245                         current -= cond * my_s;
246                         charge -= cap * my_s;
247         
248 
249                         /*
250                         function:Just for example ,for formula(9)
251                         Anxn = Gnxn + alpha*Cnxn
252                             */
253                         (*(ptr+count_maxtrix))->A[p] = cond + (*(ptr+count_maxtrix))->alpha * cap;
254                         
255                         ++p;
256                     }
257                 }
258 
259 
260                 /*
261                     function:Just for example ,for formula(7 8)
262                     Rnx1 = Dnx1 - Gnxn*Snx1
263                     Hnx1 = D'nx1 - Cnxn*Snx1
264                     */
265                 (*(ptr+count_maxtrix))->R[row] = current;
266                 (*(ptr+count_maxtrix))->H[row] = charge;
267 
268 
269 
270 
271             }
272         }
273     });

 

 

图3 TaskB程序实现流程图

2.3、线程加速仿真网表计算实现

　　实验通过创建多线程来加速仿真网表的计算，本项目中，使用了 C++ 11 支持多线程机制。通常想要在 C++程序中添加多线程则需要借助操作系统的 API 来实现，但是对于需要运行在不同平台上的底层程序，传统的实现并发功能的方法显得十分短板，future 头文件中声 明 了 std::promise,std::package_task 两 个 provider 类 ， 以 及 std::future 和std::shared_future 两个 future 类，另外还有一些相关的类型变量定义和函数。作为一种类模板的定义接口，std:: future 头文件中的各个类可以有效的获取每一个异步线程的结果，同时也为线程之间异步通信提供了手段。线程实现函数部分代码如图 4 所示：
图 4 线程实现函数部分代码
　　举个例子，如果 N 个线程的子任务之间，各自独立同时进行。那么这 N 个子任务完成的总时间 Ts 与一个子任务完成的时间相同。这种情况下，可以认为并行的加速比为 N。这个加速比是相对不并行的情况而言的。而另外一个串行的例子，如果第 2 个子任务所需要的数据需要等第 1 个子任务完成才能得到，第 3 个子任务所需要的数据需要等第 2 个子任务完成才能得到，依次类推，所需时间为NTs。具体如图5及图 6所示： 

图 5线程加速示意图

 

图 6 任务串行执行示意图

2.4、AVX 指令集仿真网表计算实现

　　相比于线程级并行，指令集并行则是在微观角度上更加细致的概念，计算机处理问题是通过指令集来实现的，而当每一个线程级任务进入系统运行请求 CPU 操作时，系统会将每一
个线程任务分解为多个指令，而在更接近底层的指令调用则需要流水线技术来实现指令集并行操作，如图7 所示。
图 7 SIMD 实现计算加速
　　具体实现依据了 Intel 公司的技术文档，定义对应的变量如__mm256d 等，对所需的指令集命令如_mm256_mul_pd 实现单指令多数据流乘法运算，具体实现如下所示：
 1 columnIndice_m =_mm256_set_pd((*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j+3]],
 2                                         (*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j+2]],
 3                                         (*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j+1]],
 4                                         (*(ptr+count_maxtrix))->S[(*(ptr+count_maxtrix))->columnIndice[j]] );
 5                     //ID_m = _mm256_load_pd((*(ptr+count_maxtrix))->Id+j);
 6                     NormalMatrix_m = _mm256_loadu_pd((*(ptr+count_maxtrix))->valueNormalMatrix+j);
 7                     
 8                 
 9                 
10                     mux_m = _mm256_mul_pd(NormalMatrix_m, columnIndice_m);
11                 
12                     //ID_mm = _mm256_add_pd(ID_m, mux_m);
13                     
14                     _mm256_storeu_pd(ID_array, mux_m);

3、TaskA 与 TaskB 任务 testApp 版本数据测试

3.1、matrix_calc_taskA 任务的实现与测试数据

　　首先是展示 matrix_calc_taskA 任务的 testApp 版本的测试结果（以-loop 1000-matrixnum 1024 相同条件），通过修改及提高算法来实现性能的优越性。最初实现 matrix_calc_taskA 计算式（6）时 version1 版本主要通过简单 for 循环遍历二级指针指向的矩阵向量，通过三层 for 循环计算出 Id 数组与两个参数指针指向的向量之间的数量关系，并分别计算相关任务，虽然 verison1 方法在精度上并无任何差错，但是version1 版本动态库在服务器上时 taskA 耗时 20s,见图8所示。
图8 TaskA-version1测试
　　Version5 版本修改将在之前的 for 循环上运用 C++语言的线程机制，将原本三个 for循环要实现的计算工作量裁分为 16 个独立线程工作量，充分利用线程的并发性来提升整体的计算效率，且在原有的 16 线程前提下增加了 AVX 指令集，提升了加速效果， 时间减少到 5.19s。
图9 TaskA-vesion5测试

3.2、TaskA 与 TaskB 任务 nanospice 电路仿真器和网表数据测试

　　因为服务器后期 cpu 占有率及个人电脑配置等问题，case1 及 case2 的每次测试时间较长，现以加速效果最好的 16 线程并行+AVX 指令集 version5 版本与官方的动态链接库进行
测试比对。case1 网表测试：使用 16 线程并行+AVX 指令集技术，精度无误差。（测试用例已改变，并与官方的动态链接库进行比对）
图 8 case1 网表 16 线程并行+AVX 指令集技术与官方库测试对比

 4、设计总结及扩展

　　在设计解决 matrix_calc_taskA 以及 matrix_calc_taskB 两种实现算法时，不断的在探索优质的解决算法，同时兼顾精度、时间、内存消耗等三个评判指标，在第五版本的设计中，加上了多线程及 AVX 指令集，通过 testApp 测试，可以发现 16 线程+AVX 指令集提升的速度是最快的，而在 nanospice 电路仿真器和网表数据测试中，可以看出 16 线程+AVX 指令集实现的动态链接库比官方的动态链接库有明显的加速效果，但二者在内存使用上相差无几，这是对赛题的算法实现最优的并行加速效果。
　　总的来说，自己感觉设计没有弄的很好，希望完善一下，具体可以参考F组第一名的答辩，可以吸收吸收.
　　重难点分析：
（1）随机性
　　稀疏矩阵非零值位置没有规律；矩阵向量规律庞大且不固定；运算行号与行数计算任务变化
（2）数据调度
　　矩阵与向量在内存上分布形式已确定为CSR；任务粒度小，突发性执行次数多，对模块相应要求高；不连续访问所造成的不规律访存，造成cache missing，带宽利用率低
　　解决方案：
（1）代码优化：空间局部优化；循环效率优化；访存效率优化
（2）多线程并行：固定任务非加锁多线程；跳跃任务非加锁多线程；基于任务队列的多线程；基于矩阵索引的线程池（最后采用）
（3）亲和性设置：CBT（core banding thread)（最后采用）；SBT（socket banding thread）------------从结果影响巨大，提高效能非常多
　总优化过程：基础代码--------->经过循环效率优化、空间局部优化、访存优化------>基础优化代码------>经过基于矩阵索引的线程池----->CBT------>基于矩阵索引的CBT线程池
 
 
 

 

 

 

 

 

posted @ 2021-01-31 23:45  `Konoha  阅读(1049)  评论(0)  编辑  收藏  举报