基于OpenMP的矩阵乘法全循环并行组合验证问题求助

背景说明 #

现有串行矩阵乘法C代码，需对三层循环（记为Z循环：i遍历行、X循环：j遍历列、Y循环：k累加计算）的所有7种并行组合（[Z]、[X]、[Y]、[Z,X]、[Z,Y]、[X,Y]、[Z,X,Y]）验证，目前仅成功实现Y循环并行，其他组合出现段错误或计算异常，以下是各组合的正确实现及错误原因分析。

串行核心代码 #

void matrix_multiply(int n, int m, int p, int **A, int **B, int **C) {
  int i,j,k,v;
  for (i = 0; i < n; i++) {
    for (j = 0; j < p; j++) {
      v = 0;
      for (k = 0; k < m; k++) {
        v += A[i][k] * B[k][j];
      }
      C[i][j] = v;
    }
  }
}

各并行组合正确实现 #

1. 仅并行Z循环（i循环） #

错误原因（原尝试代码）：

使用#pragma omp parallel而非parallel for，导致所有线程重复执行完整的i循环，多线程同时写入同一C[i][j]引发数据竞争；且错误对v使用reduction，v在循环内每次初始化，不符合规约逻辑。

正确实现：

void matrix_multiply(int n, int m, int p, int **A, int **B, int **C) {
  int i,j,k,v;
  #pragma omp parallel for num_threads(8) private(j,k,v) shared(A,B,C)
  for (i = 0; i < n; i++) {
    for (j = 0; j < p; j++) {
      v = 0;
      for (k = 0; k < m; k++) {
        v += A[i][k] * B[k][j];
      }
      C[i][j] = v;
    }
  }
}

说明：

• 每个线程负责独立的i行，j、k循环在单线程内串行执行
• j、k、v声明为线程私有，避免多线程间变量冲突
• 每个C[i][j]仅由对应线程写入，无数据竞争，无需同步

2. 仅并行X循环（j循环） #

正确实现：

void matrix_multiply(int n, int m, int p, int **A, int **B, int **C) {
  int i,j,k,v;
  for (i = 0; i < n; i++) {
    #pragma omp parallel for num_threads(8) private(k,v) shared(A,B,C,i)
    for (j = 0; j < p; j++) {
      v = 0;
      for (k = 0; k < m; k++) {
        v += A[i][k] * B[k][j];
      }
      C[i][j] = v;
    }
  }
}

说明：

• 外层i循环串行，每个线程负责同一行i下的不同j列
• k、v为线程私有，每个线程独立计算对应j的C[i][j]，无数据竞争

3. 仅并行Y循环（k循环）（已成功实现） #

错误原因（原尝试代码）：

未使用reduction子句，多线程同时写入共享变量v，引发数据竞争，导致计算结果异常。

正确实现（已验证）：

void matrix_multiply(int n, int m, int p, int **A, int **B, int **C) {
  int i,j,k,v;
  for (i = 0; i < n; i++) {
    for (j = 0; j < p; j++) {
      v = 0;
      #pragma omp parallel for reduction(+: v) schedule(static)
      for (k = 0; k < m; k++) {
        v += A[i][k] * B[k][j];
      }
      C[i][j] = v;
    }
  }
}

说明：

• reduction(+: v)让每个线程拥有独立的v副本，累加完成后合并所有副本的值，避免数据竞争

4. 并行Z+X循环（i+j循环） #

正确实现：

void matrix_multiply(int n, int m, int p, int **A, int **B, int **C) {
  int i,j,k,v;
  #pragma omp parallel for collapse(2) num_threads(8) private(k,v) shared(A,B,C)
  for (i = 0; i < n; i++) {
    for (j = 0; j < p; j++) {
      v = 0;
      for (k = 0; k < m; k++) {
        v += A[i][k] * B[k][j];
      }
      C[i][j] = v;
    }
  }
}

说明：

• 使用collapse(2)将i、j两层循环合并为一个逻辑循环，线程分配不同的(i,j)对
• 每个线程独立计算对应(i,j)的C[i][j]，无数据竞争，计算效率优于单独并行某一层

5. 并行Z+Y循环（i+k循环） #

正确实现：

void matrix_multiply(int n, int m, int p, int **A, int **B, int **C) {
  int i,j,k,v;
  #pragma omp parallel for num_threads(8) private(j,k,v) shared(A,B,C)
  for (i = 0; i < n; i++) {
    for (j = 0; j < p; j++) {
      v = 0;
      #pragma omp parallel for reduction(+: v)
      for (k = 0; k < m; k++) {
        v += A[i][k] * B[k][j];
      }
      C[i][j] = v;
    }
  }
}

说明：

• 外层并行i循环，每个线程负责一行数据
• 每行内的j循环串行，每个j对应的k循环并行，通过reduction保证v的累加正确性
• 注意：该组合嵌套并行，需确保编译器支持嵌套OpenMP并行（可通过OMP_NESTED=true环境变量开启）

6. 并行X+Y循环（j+k循环） #

正确实现：

void matrix_multiply(int n, int m, int p, int **A, int **B, int **C) {
  // 先初始化C所有元素为0
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < p; j++) {
      C[i][j] = 0;
    }
  }

  int i,j,k;
  for (i = 0; i < n; i++) {
    #pragma omp parallel for collapse(2) reduction(+: C[i][*]) num_threads(8) shared(A,B,C)
    for (j = 0; j < p; j++) {
      for (k = 0; k < m; k++) {
        C[i][j] += A[i][k] * B[k][j];
      }
    }
  }
}

说明：

• 先初始化C为全0，避免累加初始值错误
• 使用collapse(2)合并j、k循环，reduction(+: C[i][*])保证多线程对C[i][j]的累加操作安全，避免数据竞争

7. 并行Z+X+Y循环（三层全并行） #

正确实现：

void matrix_multiply(int n, int m, int p, int **A, int **B, int **C) {
  // 先初始化C所有元素为0
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < p; j++) {
      C[i][j] = 0;
    }
  }

  int i,j,k;
  #pragma omp parallel for collapse(3) reduction(+: C[:n][:p]) num_threads(8) shared(A,B,C)
  for (i = 0; i < n; i++) {
    for (j = 0; j < p; j++) {
      for (k = 0; k < m; k++) {
        C[i][j] += A[i][k] * B[k][j];
      }
    }
  }
}

说明：

• 三层循环通过collapse(3)合并，线程分配不同的(i,j,k)三元组
• 利用reduction保证C[i][j]的累加安全，需确保编译器支持二维数组的规约操作
• 若编译器不支持二维数组规约，可改为先并行i+j循环，内层k循环串行，等价于Z+X并行方案

内容的提问来源于stack exchange，提问作者user14054021