今年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線程池