将 OpenMP 添加到程序中以计算 n x n 矩阵 n x n 的行列式

Add OpenMP to program to calculate the determinant of an n x n matrix n x n

这是查找矩阵 n x n 的行列式的代码。

#include <iostream>

using namespace std;

int determinant(int *matrix[], int size);
void ijMinor(int *matrix[], int *minorMatrix[], int size, int row, int column);

int main()
{
    int size;
    cout << "What is the size of the matrix for which you want to find the determinant?:\t";
    cin >> size;

    int **matrix;
    matrix = new int*[size];
    for (int i = 0 ; i < size ; i++)
        matrix[i] = new int[size];

    cout << "\nEnter the values of the matrix seperated by spaces:\n\n";
    for(int i = 0; i < size; i++)
        for(int j = 0; j < size; j++)
            cin >> matrix[i][j];

    cout << "\nThe determinant of the matrix is:\t" << determinant(matrix, size) << endl;

    return 0;
}

int determinant(int *matrix[], int size){
    if(size==1)return matrix[0][0];
    else{
        int result=0, sign=-1;
        for(int j = 0; j < size; j++){

            int **minorMatrix;
            minorMatrix = new int*[size-1];
            for (int k = 0 ; k < size-1 ; k++)
                minorMatrix[k] = new int[size-1];

            ijMinor(matrix, minorMatrix, size, 0, j);

            sign*=-1;
            result+=sign*matrix[0][j]*determinant(minorMatrix, size-1);
            for(int i = 0; i < size-1; i++){
                delete minorMatrix[i];
            }
        }

        return result;
    }
}

void ijMinor(int *matrix[], int *minorMatrix[], int size, int row, int column){
    for(int i = 0; i < size; i++){
        for(int j = 0; j < size; j++){
            if(i < row){
                if(j < column)minorMatrix[i][j] = matrix[i][j];
                else if(j == column)continue;
                else minorMatrix[i][j-1] = matrix[i][j];
            }
            else if(i == row)continue;
            else{
                if(j < column)minorMatrix[i-1][j] = matrix[i][j];
                else if(j == column)continue;
                else minorMatrix[i-1][j-1] = matrix[i][j];
            }
        }
    }
}

添加 OpenMP 编译指示后,我更改了行列式函数,现在它看起来像这样:

int determinant(int *matrix[], int size){
    if(size==1)return matrix[0][0];
    else{
        int result=0, sign=-1;
        #pragma omp parallel for default(none) shared(size,matrix,sign) private(j,k)  reduction(+ : result)

        for(int j = 0; j < size; j++){

            int **minorMatrix;

            minorMatrix = new int*[size-1];

            for (int k = 0 ; k < size-1 ; k++)
                minorMatrix[k] = new int[size-1];

            ijMinor(matrix, minorMatrix, size, 0, j);

            sign*=-1;
            result+=sign*matrix[0][j]*determinant(minorMatrix, size-1);
            for(int i = 0; i < size-1; i++){
                delete minorMatrix[i];
            }
        }

        return result;    
        delete [] matrix;
    }
}

我的问题是每次的结果都不一样。有时它给出正确的值,但大多数时候它是错误的。我认为这是因为 sign 变量。我遵循以下公式:

如您所见,在我的 for 循环的每次迭代中都应该有不同的 sign 但是当我使用 OpenMP 时,出现了问题。如何使用 OpenMP 将此程序制作成 运行?

最后,我的第二个问题是使用 OpenMP 并不会使程序 运行 比不使用 OpenMP 更快。我也尝试制作一个 100,000 x 100,000 的矩阵,但是我的程序报告分配内存的错误。我如何 运行 这个程序有非常大的矩阵?

我看到你的问题如下:

1) 正如 Hristo 所指出的,关于 sign 变量,您的线程正在互相踩踏对方的数据。它应该对每个线程都是私有的,这样它们就可以完全 read/write 访问它而不必担心竞争条件。然后,您只需要一个算法来计算 sign 是正负 1,具体取决于迭代 j 独立于其他迭代 。稍微思考一下,您就会发现 Hristo 的建议是正确的:sign = (j % 2) ? -1 : 1; 应该可以解决问题。

2) 您的 determinant() 函数是递归的。照原样,这意味着循环的每次迭代,在形成你的未成年人之后,你然后在那个未成年人上再次调用你的函数。因此,单个线程将执行它的迭代,进入递归函数,然后尝试将自己拆分为 nthreads 个更多线程。 你现在可以看到你是如何通过启动比您实际拥有的内核更多的线程来超额订阅您的系统。两个简单的解决方案:

  • omp parallel 代码中调用您的原始串行函数。这是最快的方法,因为这可以避免任何 OpenMP-startup 开销。
  • 在首次调用 determinant() 之前调用 omp_set_nested(0); 关闭嵌套并行性。
  • 在并行指令中添加一个 if 子句:if(omp_in_parallel())

3) 你的内存问题是因为递归的每次迭代,你都在分配更多的内存。如果您解决了问题 #2,那么您应该在串行情况下使用与并行情况相当的内存量。话虽如此,在进入算法之前分配所有你想要的内存会好得多。分配大块内存(然后释放它!),尤其是并行分配,是代码中的一个可怕瓶颈。

在进入第一个循环之前(在纸上)计算您需要的内存量,并一次分配所有内存。我还强烈建议您考虑连续分配内存(又名一维)以更好地利用缓存。请记住,每个线程都应该有自己独立的工作区域。然后,将您的功能更改为:

int determinant(int *matrix, int *startOfMyWorkspace, int size)

无需在循环内分配新的 (size-1)x(size-1) 矩阵,您只需利用工作区的下一个 (size-1)*(size-1) 整数,更新下一个递归的 startOfMyWorkspace打电话,继续。