C OMP for 并行区域中的循环。非工作共享

Question

我有一个想要并行化的函数。这是连续版。

void parallelCSC_SpMV(float *x, float *b)
{
    int i, j;
    for(i = 0; i < numcols; i++)
    {
        for(j = colptrs[i] - 1; j < colptrs[i+1] - 1; j++)
        {
            b[irem[j] - 1] += xrem[j]*x[i];
        }
    }
}

我想一个不错的方法是让每个线程写入 b 数组的私有副本（它不需要是受保护的关键部分，因为它是私有副本），在线程完成后，它然后将其结果复制到实际的 b 数组。这是我的代码。

void parallelCSC_SpMV(float *x, float *b)
{
    int i, j, k;
    #pragma omp parallel private(i, j, k)
    {
        float* b_local = (float*)malloc(sizeof(b));       
     
        #pragma omp for nowait
        for(i = 0; i < numcols; i++)
        {
            for(j = colptrs[i] - 1; j < colptrs[i+1] - 1; j++)
            {
                float current_add = xrem[j]*x[i];
                int index = irem[j] - 1;
                b_local[index] += current_add;
            }
        }
        
        for (k = 0; k < sizeof(b) / sizeof(b[0]); k++)
        {
            // Separate question: Is this if statement allowed?
            //if (b_local[k] == 0) { continue; }
            #pragma omp atomic
            b[k] += b_local[k];
        }
    }
}

但是，由于第二个 for 循环，我得到了分段错误。我不需要在那个循环上使用"#pragma omp for"，因为我希望每个线程都能完全执行它。如果我注释掉 for 循环内的内容，则不会出现分段错误。我不确定会是什么问题。

Answer 1

首先，正如Jim Cownie所指出的：

在所有这些答案中，b_local 未初始化，但您正在添加 给它。你需要使用 calloc 而不是 malloc

只是添加到接受的答案，我认为您可以尝试以下方法以避免并行调用 malloc，以及调用 #pragma omp atomic 的开销。

void parallelCSC_SpMV(float *x, float *b, int b_size, int num_threads) {
    
    float* b_local[num_threads];
    for(int i = 0; i < num_threads; i++) 
       b_local[i] = calloc(b_size, sizeof(float));
    
    #pragma omp parallel num_threads(num_threads)
    { 
        int tid = omp_get_thread_num();
        #pragma omp for
        for(int i = 0; i < numcols; i++){ 
            for(int j = colptrs[i] - 1; j < colptrs[i+1] - 1; j++){
               float current_add = xrem[j]*x[i];
               int index = irem[j] - 1;
               b_local[tid][index] += current_add;
           }
       }
    }   
    for(int id = 0; id < num_threads; id++)
    {   
        #pragma omp for simd
        for (int k = 0; k < b_size; k++)
        {    
             b[k] += b_local[id][k];
        }
        free(b_local[id]);
    }
}  

我还没有测试过它的性能，所以请随时这样做并提供反馈。

您可以进一步优化，而不是为主线程创建一个local_b，只是复用了原来的b，如下：

void parallelCSC_SpMV(float *x, float *b, int b_size, int num_threads) {
    
    float* b_local[num_threads-1];
    for(int i = 0; i < num_threads-1; i++) 
       b_local[i] = calloc(b_size, sizeof(float));
    
    #pragma omp parallel num_threads(num_threads)
    { 
        int tid = omp_get_thread_num();
        float *thread_b = (tid == 0) ? b : b_local[tid-1];
        #pragma omp for
        for(int i = 0; i < numcols; i++){ 
            for(int j = colptrs[i] - 1; j < colptrs[i+1] - 1; j++){
               float current_add = xrem[j]*x[i];
               int index = irem[j] - 1;
               thread_b[index] += current_add;
           }
       }
    }   
    
    for(int id = 0; id < num_threads-1; id++)
    {   
        #pragma omp for simd
        for (int k = 0; k < b_size; k++)
        {    
             b[k] += b_local[id][k];
        }
        free(b_local[id]);
    }
}  

Answer 2

您可能正在尝试访问动态数组 b_local 中超出范围的位置。

看到sizeof(b) 将返回float* 的字节大小（浮点指针的大小）。

如果您想知道要传递给函数的数组的大小，我建议您将其添加到函数的参数中。

void parallelCSC_SpMV(float *x, float *b, int b_size){
...
    float* b_local = (float*) malloc(sizeof(float)*b_size); 
...
}

而且，如果colptrs 的大小是numcols，我会小心colptrs[i+1]，因为当i=numcols-1 会出现另一个超出范围的问题。