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

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

看到 sizeof(b) 将 return 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 会出现另一个超出范围的问题。

Answer 2

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

In all of these answers, b_local is uninitialised, yet you are adding to it. You need to use calloc instead of 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]);
    }
}

C OMP for 并行区域中的循环。不分工

C OMP for loop in parallel region. Not work-shared

c

parallel-processing

performance

multithreading

openmp