c++ - 三角矩阵转换和自动并行化

Question

我在 ICC 中玩了一下自动并行化(11.1；旧的，但对此无能为力)，我想知道为什么编译器不能并行化内部循环以进行简单的高斯消除:

void makeTriangular(float **matrix, float *vector, int n) {
    for (int pivot = 0; pivot < n - 1; pivot++) {
        // swap row so that the row with the largest value is
        // at pivot position for numerical stability
        int swapPos = findPivot(matrix, pivot, n);
        std::swap(matrix[pivot], matrix[swapPos]);
        std::swap(vector[pivot], vector[swapPos]);
        float pivotVal = matrix[pivot][pivot];
        for (int row = pivot + 1; row < n; row++) { // line 72; should be parallelized
            float tmp = matrix[row][pivot] / pivotVal;  
            for (int col = pivot + 1; col < n; col++) { // line 74
                matrix[row][col] -= matrix[pivot][col] * tmp;
            }
            vector[row] -= vector[pivot] * tmp;
        }
    }
}

我们只写入依赖于私有(private)行(和列)变量的数组，并且行保证大于枢轴，因此对于编译器来说应该很明显我们没有覆盖任何内容。

我正在使用-O3 -fno-alias -parallel -par-report3进行编译，并获得大量依赖项ala:假设矩阵行75和矩阵行73之间存在FLOW依赖关系。 或假定矩阵行 73 和矩阵行 75 之间的 ANTI 依赖性。 并且仅对第 75 行相同。编译器有什么问题？显然我可以准确地告诉它如何处理一些编译指示，但我想了解编译器可以单独得到什么。

Answer 1

基本上，编译器无法确定不存在依赖关系，因为名称 matrix 和名称 vector 都被读取和写入(即使使用不同的地区)。您也许可以通过以下方式解决这个问题(虽然有点脏):

void makeTriangular(float **matrix, float *vector, int n)
{     
    for (int pivot = 0; pivot < n - 1; pivot++) 
    {         
         // swap row so that the row with the largest value is    
         // at pivot position for numerical stability       
         int swapPos = findPivot(matrix, pivot, n);    
         std::swap(matrix[pivot], matrix[swapPos]);   
         std::swap(vector[pivot], vector[swapPos]);     
         float pivotVal = matrix[pivot][pivot];     
         float **matrixForWriting = matrix;  // COPY THE POINTER
         float *vectorForWriting = vector;   // COPY THE POINTER
         // (then parallelize this next for loop as you were)
         for (int row = pivot + 1; row < n; row++)  { 
              float tmp = matrix[row][pivot] / pivotVal;               
              for (int col = pivot + 1; col < n; col++) {
                  // WRITE TO THE matrixForWriting VERSION
                  matrixForWriting[row][col] = matrix[row][col] - matrix[pivot][col] * tmp; 
              } 
              // WRITE TO THE vectorForWriting VERSION
              vectorForWriting[row] = vector[row] - vector[pivot] * tmp; 
         } 
    }
} 

底线只是给你正在写入的内容一个临时不同的名称来欺骗编译器。我知道这有点肮脏并且我一般不会推荐这种编程。但如果您确定没有数据依赖性，那就完全没问题了。

事实上，我会在它周围添加一些注释，对于 future 看到此代码的人来说，这是一种解决方法，以及为什么要这样做。

编辑:我认为这个答案基本上是由@FPK 触及的，并且答案是由@Evgeny Kluev 发布的。然而，在 @Evgeny Kluev 的回答中，他建议将此作为输入参数，这可能会并行化，但不会给出正确的值，因为矩阵中的条目不会更新。我认为我上面发布的代码也会给出正确的答案。