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

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;
        }
    }
}

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

我正在使用-O3 -fno-alias -parallel -par-report3 进行编译，并获得了很多依赖项ala：assumed FLOW dependence between matrix line 75 and matrix line 73. 或assumed ANTI dependence between matrix line 73 and matrix line 75.，仅第75 行也是如此。编译器有什么问题？显然，我可以准确地告诉它如何处理一些 pragma，但我想了解编译器可以单独获得什么。

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; 
         } 
    }
} 

底线只是给你正在写的那些暂时不同的名字来欺骗编译器。我知道这有点脏而且我一般不会推荐这种编程方式。但是如果你确定你没有数据依赖，那很好。

事实上，我会在它周围放置一些 cmets，以便未来看到此代码的人非常清楚这是一种解决方法，以及您这样做的原因。

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

Answer 2

icc 12.1 上存在相同的自动并行化问题。所以我用这个较新的版本进行实验。

将输出矩阵添加到函数的参数列表并将第三个循环的主体更改为此

out[row][col] = matrix[row][col] - matrix[pivot][col] * tmp;

修复了“FLOW 依赖”问题。这意味着，“-fno-alias”只影响函数参数，而单个参数的内容仍然被怀疑是别名。我不知道为什么这个选项不会影响一切。由于矩阵的不同部分并没有真正相互别名，您可以将这个附加参数留给函数并通过这个参数传递相同的矩阵。

有趣的是，在抱怨 'matrix' 时，编译器对 'vector' 只字未提，这确实存在别名问题：这条线 vector[row] -= vector[pivot] * tmp; 可能会导致错误别名（在一个线程中写入 vector[row] 可能会触及缓存线，存储vector[pivot]，供每个线程使用）。

“FLOW 依赖”不是这段代码中唯一的问题。修复后，编译器仍然拒绝并行化第二个和第三个循环，因为“计算工作不足”。所以我试着给它一些额外的工作：

float tmp = matrix[row][pivot] * pivotVal;
...
out[row][col] = matrix[row][col] - matrix[pivot][col] *tmp /pivotVal /pivotVal;

在这一切之后，第二个循环终于并行化了，虽然我不确定它是否获得了任何速度提升。

---

更新：我找到了一个更好的替代方案，可以让计算机“做一些额外的工作”。选项-par-threshold50 可以解决问题。

Answer 3

我无法访问 icc 来测试我的想法，但我怀疑编译器担心别名：矩阵定义为浮点**：指向浮点数组的指针数组。所有这些指针都可以指向同一个浮点数组，因此并行化这将非常危险。这没有任何意义，但编译器无法知道。