OpenMP程序（k均值++）不結垢

Question

目前，我編程的k-均值++使用OpenMP和C.並行版本直到現在，我實施質心的初始化。如果您對此過程不熟悉，則其工作原理與follows大致相同。給定一個dataset（基體）與n點，k質心是使用「概率函數」，也被稱爲輪盤賭選擇initilized。

假設你有n=4點，距離下面的數組一些重心：

distances = [2, 4, 6, 8] 
dist_sum = 20 

從這些，除以distances每個條目由dist_sum，並添加以前的結果，就像定義一個累積概率數組這：

probs  = [0.1, 0.2, 0.3, 0.4] = [2/20, 4/20, 6/20, 8/20] 
acc_probs = [0.1, 0.3, 0.6, 1.0] 

然後，執行輪盤賭選擇。給定一個隨機數，說r=0.5，選擇使用r和acc_probs，下一個點迭代acc_probs直到r < acc_probs[i]。在這個例子中，所選擇的點是因爲i=2r < acc_probs[2]。

問題 在這種情況下，我非常大的矩陣（約點n=16 000 000）工作。儘管該程序給出了正確的答案（即重心的良好初始化），但它並沒有如預期的那樣成比例。該函數按照此算法計算初始質心。

double **parallel_init_centroids (double **dataset, int n, int d, int k, RngStream randomizer, long int *total_ops) { 

    double dist=0, error=0, dist_sum=0, r=0, partial_sum=0, mindist=0; 
    int cn=0, cd=0, ck = 0, cck = 0, idx = 0; 
    ck = 0; 

    double probs_sum = 0; // debug 

    int mink=0, id=0, cp=0; 

    for (ck = 0; ck < k; ck++) { 

     if (ck == 0) { 

      // 1. choose an initial centroid c_0 from dataset randomly 
      idx = RngStream_RandInt (randomizer, 0, n-1); 

     } 
     else { 

      // 2. choose a successive centroid c_{ck} using roulette selection 
      r = RngStream_RandU01 (randomizer); 
      idx = 0; 
      partial_sum = 0; 
      for (cn=0; cn<n; cn++) { 

       partial_sum = partial_sum + distances[cn]/dist_sum; 

       if (r < partial_sum) { 
        idx = cn; 
        break; 
       } 
      } 
     } 

     // 3. copy centroid from dataset 
     for (cd=0; cd<d; cd++) 
      centroids[ck][cd] = dataset[idx][cd]; 

     // reset before parallel region 
     dist_sum = 0; 

     // -- parallel region -- 
     # pragma omp parallel shared(distances, clusters, centroids, dataset, chunk, dist_sum_threads, total_ops_threads) private(id, cn, cck, cd, cp, error, dist, mindist, mink) 
     { 
      id = omp_get_thread_num(); 
      dist_sum_threads[id] = 0;    // each thread reset its entry 

      // parallel loop 
      // 4. recompute distances against centroids 
      # pragma omp for schedule(static,chunk) 
      for (cn=0; cn<n; cn++) { 

       mindist = DMAX; 
       mink = 0; 

       for (cck=0; cck<=ck; cck++) { 

        dist = 0; 

        for (cd=0; cd<d; cd++) { 

         error = dataset[cn][cd] - centroids[ck][cd]; 
         dist = dist + (error * error);     total_ops_threads[id]++; 
        } 

        if (dist < mindist) { 
         mindist = dist; 
         mink = ck; 
        } 
       } 

       distances[cn]   = mindist; 
       clusters[cn]   = mink; 
       dist_sum_threads[id] += mindist;  // each thread contributes before reduction 
      } 
     } 
     // -- parallel region -- 

     // 5. sequential reduction 
     dist_sum = 0; 
     for (cp=0; cp<p; cp++) 
      dist_sum += dist_sum_threads[cp]; 
    } 


    // stats 
    *(total_ops) = 0; 
    for (cp=0; cp<p; cp++) 
     *(total_ops) += total_ops_threads[cp]; 

    // free it later 
    return centroids; 
} 

正如你所看到的，並行區域計算中nd維點對kd維質心的距離。這項工作在p線程（最多32個）之間共享。並行區域結束後，兩個陣列被填充：distances和dist_sum_threads。第一個數組與前一個示例相同，而第二個數組包含每個線程收集的累積距離。考慮前面的例子，如果p=2線程是可用的，那麼該陣列被定義如下：

dist_sum_threads[0] = 6 ([2, 4]) # filled by thread 0 
dist_sum_threads[1] = 14 ([6, 8]) # filled by thread 1 

dist_sum通過加入的dist_sum_threads每個條目定義。該函數按預期工作，但是當線程數量增加時，執行時間會增加。這figure顯示了一些性能指標。

出了什麼問題我的執行，特別是使用OpenMP？綜上所述，採用了只有兩個雜注：

# pragma omp parallel ... 
{ 
    get thread id 

    # pragma omp for schedule(static,chunk) 
    { 
     compute distances ... 
    } 

    fill distances and dist_sum_threads[id] 
} 

換句話說，我消除了障礙，互斥訪問，以及其他的東西，可能會導致額外的開銷。但是，執行時間隨着線程數量的增加而變得最差。

更新

• 上面的代碼已經被更改爲mcve。 This snippet與我以前的代碼類似。在這種情況下，計算n=100000點和k=16質心之間的距離。
• 執行時間在並行區域之前和之後使用omp_get_wtime進行測量。總時間存儲在wtime_spent中。
• 我包含一個減少計算dist_sum。然而，它並不像預期的那樣工作（它在下面被評論爲不好的並行減少）。 dist_sum的正確值是999857108020.0，但是當使用p線程來計算它時，結果爲999857108020.0 * p，這是錯誤的。
• 性能情節updated

• 這是主要的並行功能，完整的代碼位於here：

  double **parallel_compute_distances (double **dataset, int n, int d, int k, long int *total_ops) { 
  
      double dist=0, error=0, mindist=0; 
  
      int cn, cd, ck, mink, id, cp; 
  
      // reset before parallel region 
      dist_sum = 0;    
  
      // -- start time -- 
      wtime_start = omp_get_wtime(); 
  
      // parallel loop 
      # pragma omp parallel shared(distances, clusters, centroids, dataset, chunk, dist_sum, dist_sum_threads) private(id, cn, ck, cd, cp, error, dist, mindist, mink) 
      { 
       id = omp_get_thread_num(); 
       dist_sum_threads[id] = 0;    // reset 
  
       // 2. recompute distances against centroids 
       # pragma omp for schedule(static,chunk) 
       for (cn=0; cn<n; cn++) { 
  
        mindist = DMAX; 
        mink = 0; 
  
        for (ck=0; ck<k; ck++) { 
  
         dist = 0; 
  
         for (cd=0; cd<d; cd++) { 
  
          error = dataset[cn][cd] - centroids[ck][cd]; 
          dist = dist + (error * error);        total_ops_threads[id]++; 
         } 
  
         if (dist < mindist) { 
          mindist = dist; 
          mink = ck; 
         } 
        } 
  
        distances[cn]   = mindist; 
        clusters[cn]   = mink; 
        dist_sum_threads[id] += mindist; 
       } 
  
  
       // bad parallel reduction 
       //#pragma omp parallel for reduction(+:dist_sum) 
       //for (cp=0; cp<p; cp++){    
       // dist_sum += dist_sum_threads[cp]; 
       //} 
  
      } 
  
  
      // -- end time -- 
      wtime_end = omp_get_wtime(); 
  
      // -- total wall time -- 
      wtime_spent = wtime_end - wtime_start; 
  
      // sequential reduction 
      for (cp=0; cp<p; cp++)  
       dist_sum += dist_sum_threads[cp]; 
  
  
      // stats 
      *(total_ops) = 0; 
      for (cp=0; cp<p; cp++) 
       *(total_ops) += total_ops_threads[cp]; 
  
      return centroids; 
  } 
  

Answer 1

您的代碼不是一個mcve我只能在這裏發出的假說。然而，這裏就是我的想法（可能）發生（按重要性沒有特定的順序）：

• 你從錯誤共享苦，當你更新dist_sum_threads和total_ops_threads。只需在parallel區域內聲明reduction(+: dist_sum)並直接使用dist_sum即可完全避免前者。您也可以使用同樣使用本地total_ops宣佈爲reduction(+)，並且您最後將其累積到*total_ops中。 （順便說一句，dist_sum是計算，但從來沒有用過......）

• 該代碼看起來無論如何都是內存綁定，因爲你有大量的內存訪問幾乎沒有計算。因此，預期的加速將主要取決於您的內存帶寬和可以在並行代碼時訪問的內存控制器的數量。有關更多詳細信息，請參閱this epic answer。

• 鑑於上述可能遇到的問題的內存綁定特性，請嘗試使用內存放置（numactl可能和/或與proc_bind的線程關聯）。您也可以嘗試使用線程調度策略，或嘗試查看是否有一些循環平鋪無法應用於將數據阻塞到緩存中的問題。

• 您沒有詳細說明測量時間的方式，但請注意，加速時間僅適用於掛鐘時間而非CPU時間。請使用omp_get_wtime()進行任何此類測量。

試着解決這些問題，並根據您的內存架構評估您的實際潛在加速。如果你仍然覺得自己沒有達到你應該做的，那麼就更新你的問題。

---

編輯：

既然你提供了一個完整的例子，我成功地試驗了一下你的代碼，並實施修改我腦子裏想的（減少僞共享居多）。

這裏是什麼功能沒有的樣子：

double **parallel_compute_distances(double **dataset, int n, int d, 
            int k, long int *total_ops) { 
    // reset before parallel region 
    dist_sum = 0;    

    // -- start time -- 
    wtime_start = omp_get_wtime(); 

    long int tot_ops = 0; 

    // parallel loop 
    # pragma omp parallel for reduction(+: dist_sum, tot_ops) 
    for (int cn = 0; cn < n; cn++) { 
     double mindist = DMAX; 
     int mink = 0; 
     for (int ck = 0; ck < k; ck++) { 
      double dist = 0; 
      for (int cd = 0; cd < d; cd++) { 
       double error = dataset[cn][cd] - centroids[ck][cd]; 
       dist += error * error; 
       tot_ops++; 
      } 
      if (dist < mindist) { 
       mindist = dist; 
       mink = ck; 
      } 
     } 
     distances[cn] = mindist; 
     clusters[cn] = mink; 
     dist_sum += mindist; 
    } 

    // -- end time -- 
    wtime_end = omp_get_wtime(); 

    // -- total wall time -- 
    wtime_spent = wtime_end - wtime_start; 

    // stats 
    *(total_ops) = tot_ops; 

    return centroids; 
} 

所以，談幾點看法：

• 如前所述，dist_sum和操作總數的局部變量（tot_ops ）現在宣佈爲reduction(+:)。這避免了每個索引使用一個線程訪問相同的數組，這觸發了false sharing（這幾乎是觸發它的最佳情況）。我使用了局部變量而不是total_ops作爲指針，它不能直接在reduction子句中使用。但是，最後用tot_ops更新它可以完成這項工作。

• 我儘可能延遲了所有的變量聲明。這是一個很好的做法，因爲它可以避免大部分private聲明，這些聲明通常是OpenMP程序員的主要缺陷。現在你只需要考慮兩個reduction變量和兩個數組，它們顯然是shared，因此不需要任何額外的聲明。這簡化了很多parallel指令，並有助於集中於哪些事項

• 現在不需要線程ID的越多，parallel和for指令可以合併爲更好的可讀性（以及可能的表現太）。

• 我刪除了schedule子句，讓編譯器和/或運行時庫使用它們的默認值。如果我有充分的理由，我只會選擇不同的調度策略。

就這樣，在我的雙核筆記本電腦唱GCC 5.3.0，並與-std=c99 -O3 -fopenmp -mtune=native -march=native編譯，我得到不同的線程數和2倍的2個線程加速一致的結果。

在一個10核的機器，使用英特爾編譯器和-std=c99 -O3 -xhost -qopenmp，我從1到10個線程的線性加速...

即使是在至強融核KNC我能得到接近線性的速度 - 從1到60個線程（然後使用更多的硬件線程仍然提供一些加速，但不是相同的擴展）。

觀察到的加速度使我意識到，與我所設想的不同，代碼不受內存限制，因爲您訪問的數組實際上被很好地緩存。原因在於您只能訪問dataset[cn][cd]和centroids[ck][cd]，其中第二個維度非常小（40和16），因此很適合緩存，而要加載下一個cn索引的大塊dataset可以有效地預取。

您仍然遇到此版本的代碼的可擴展性問題？