多核中的並行字首和計算
字首和計算在平行計算中很有用，因為在處理負載平衡問題時，經常需要將若干段資料重新平分，計算字首和通常是一種有效的將資料平分的方法。例如在並行基數排序中，就會用到了字首和的計算。而下面先看看單執行緒環境中的序列字首和計算。
如果給定一個數列a[n]，令S[k] ＝ <?xml:namespace prefix = v ns = "urn:schemas-microsoft-com:vml" /><?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />a[0]+a[1]+...+a[k]，（
a[4] = {1,2,3,4};
其字首和為
S[0] = a[0] = 1;
S[1] = a[0] + a[1] = 1+ 2 = 3;
S[2] = a[0] + a[1] + a[2] = 1 + 2 + 3 = 6;
S[3] = a[0] + a[1] + a[2] + a[3] = 1 + 2 + 3 + 4 = 10;
字首和的計算非常簡單，一般地，可以用下面的函式來進行序列字首和的計算：
/**序列字首和計算函式
@paramT * pInput - 輸入資料
@paramT *pOutput -
@paramint nLen - 輸入資料長度
@returnvoid - 無
*/
template <class T>
void Sequential_PrefixSum(T * pInput, T *pOutput, int nLen)
{
int i;
pOutput[0] = pInput[0];
for ( i = 1; i < nLen; i++ )
{
pOutput[i] = pInput[i] + pOutput[i-1];
}
}
在上面的序列字首和計算程式碼中可以看出，每次迴圈都依賴於上一次迴圈的結果，因此無法直接對迴圈進行並行化，要進行並行化則必須修改計算方法，下面就來看如何進行並行字首和的計算。
字首和的平行計算方法有許多種，David Callahan的論文“Recognizing and Parallelizing Bounded Recurrences”中給出了一種適合共享儲存多處理器系統中的有界遞迴計算的通用方法，字首和計算屬於有界遞迴計算的一個特例。下面先以一個例項來講解整個平行計算的過程：
例：假設有4個處理器要計算16個整數的字首和，這16個整數如下：
12345678910111213141516
第1步，先將上面資料平分成4組，每個處理器各計算一組資料的字首和，如下所示：
（1234）（5678）（9101112）（13141516）
（13610）（5111826）（91930 42）（13274258）
第2步，選取每組資料的最後一個數據，對這幾個資料計算字首和，如下所示：
（13610）（5111826）（9193042）（13274258）
（13610）（5111836）（9193078）（132742136）
經過這步的計算後，可以很容易發現，每組的最後一個數據的值已經變成了原始資料在它所處位置之前（包含本位置）的所有資料的和。例如：
36 ＝ 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8
78 = 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12
第3步：從第2組數開始，將每組中的數（除最後一個數外）加上它的前一組數的最後一個數，即可得到所有數的字首和。如下所示：
(1 3 6 10) (5+10 11+10 18+10 36) (9+36 19+36 30+36 78)(13+78 27+78 42+78 136）
(13610) (15212836) (45556678)(91105120136）
從上面的計算過程可以看出，第1步和第3步都是很容易進行並行化計算，第2步中，由於計算量非常小，用序列計算即可，下面給出上面處理過程的程式碼實現：
#defineMIN_PRARLLEL_PREFIXSUM_COUNT8192
/**並行字首和計算函式
@paramT * pInput - 輸入資料
@paramT *pOutput - 輸出資料（字首和）
@paramint nLen - 輸入資料長度
@returnvoid - 無
*/
template<class T>
void Parallel_PrefixSum(T * pInput, T *pOutput, int nLen)
{
int i;
int nCore = omp_get_num_procs();
if ( nCore < 4 || nLen < MIN_PRARLLEL_PREFIXSUM_COUNT )
{
Sequential_PrefixSum(pInput, pOutput, nLen);
return;
}
int nStep = nLen / nCore;
if ( nStep * nCore < nLen )
{
nStep += 1;
}
#pragma omp parallel for num_threads(nCore)
for ( i = 0; i < nCore; i++ )
{
int k;
int nStart = i * nStep;
int nEnd = (i+1) * nStep;
if ( nEnd > nLen )
{
nEnd = nLen;
}
pOutput[nStart] = pInput[nStart];
for ( k = nStart+1; k < nEnd; k++ )
{
pOutput[k] = pInput[k] + pOutput[k-1];
}
}
for ( i = 2; i < nCore; i++ )
{
pOutput[i * nStep - 1] += pOutput[(i-1) * nStep - 1];
}
pOutput[nLen-1] += pOutput[(nCore-1)*nStep - 1];
#pragma omp parallel for num_threads(nCore-1)
for ( i = 1; i < nCore; i++ )
{
int k;
int nStart = i * nStep;
int nEnd = (i+1) * nStep - 1;
if ( nEnd >= nLen )
{
nEnd = nLen - 1;
}
for ( k = nStart; k < nEnd; k++ )
{
pOutput[k] += pOutput[nStart-1];
}
}
return;
}
從上面並行字首和的計算過程可以看出，它的計算量比序列字首和的計算增加了差不多一倍，如果考慮程式中的實際開銷，計算增加量還不止一倍。因此在雙核CPU機器上，使用並行字首和計算沒有任何意義，只有在超過4核CPU機器上，它才有實用價值。
Parallel_PrefixSum()函式中，先是判斷CPU核數是否小於4，並且判斷了計算量是否不足，如果兩個判斷條件中有任何一個滿足，則呼叫序列字首核計算函式進行計算，然後才進行並行字首和的計算。在平行計算時只是簡單地將計算平攤到各個CPU上，沒有考慮CPU核數較多情況下計算量平攤到各個CPU核上，執行緒粒度過小的問題，主要是為了不使程式碼看起來過於繁瑣。如有需要可以修改成自動計算出最合適的執行緒數量（可能小於CPU核數），然後平行計算時使用相應的執行緒數量即可。