OpenMP实现拉普拉斯矩阵行列式计算遇段错误（core dumped）求助

我尝试用OpenMP编写基于拉普拉斯法的矩阵行列式并行计算代码，但运行时出现Segmentation fault (core dumped)错误，怀疑问题出在#pragma omp parallel private中用于递归赋值的变量上，求帮忙排查错误。

以下是我的代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

long double **nowaMacierz(int stopien);

void usunMacierz(long double **macierz);

long double metodaLaplace(long double **macierz, int stopien, int count);

int main() {
    long double **macierz;
    int stopien;
    int count = 1;
    int nr_wiersza, nr_kolumny;

    printf("Proszę podać stopień macierzy n=");
    fscanf(stdin, "%d", &stopien);
    macierz = nowaMacierz(stopien);
    for (nr_wiersza = 0; nr_wiersza < stopien; nr_wiersza++)
        for (nr_kolumny = 0; nr_kolumny < stopien; nr_kolumny++) {
            printf("A[%d,%d]=", nr_wiersza, nr_kolumny);
            fscanf(stdin, "%Lf", &macierz[nr_wiersza][nr_kolumny]);
        }
    printf("det(A) = %Lf\n", metodaLaplace(macierz, stopien, count));
    usunMacierz(macierz);
    getchar();
    return 0;
}

long double **nowaMacierz(int stopien) {
    long double **macierz;
    int nr_wiersza;

    macierz = (long double **) calloc(stopien, sizeof(long double *));
    *macierz = (long double *) calloc(stopien * stopien, sizeof(long double));

    // Przydzielone adresy komórek pamięci są segragowane, tak by utworzyły tablicę dwuwymiarową
    for (nr_wiersza = 1; nr_wiersza < stopien; nr_wiersza++)
        *(macierz + nr_wiersza) = *(macierz + nr_wiersza - 1) + stopien;
    return macierz;
}

void usunMacierz(long double **macierz) {
    free(*macierz), free(macierz);
    return;
}

long double metodaLaplace(long double **macierz, int stopien, int count) {
    long double **dopelnienie;   
    int nr_wiersza, nr_kolumny;
    int nr_kolumny_dop, nr_kolumny_mac;
    long double det = 0.00;   

    printf("test");

    if (count == 0) {
        printf("sekwenc");
        if (stopien <= 2)
            return macierz[0][0] * (stopien > 1 ? macierz[1][1] : 1) -
                   (stopien > 1 ? macierz[0][1] * macierz[1][0] : 0);

        dopelnienie = nowaMacierz(stopien - 1);

        for (nr_kolumny = 0; nr_kolumny < stopien; nr_kolumny++) {
            for (nr_kolumny_dop = 0, nr_kolumny_mac = 0;
                 nr_kolumny_dop < stopien - 1; nr_kolumny_dop++, nr_kolumny_mac++) {
                nr_kolumny_mac += (nr_kolumny_mac == nr_kolumny ? 1
                                                                : 0);
                for (nr_wiersza = 0; nr_wiersza < stopien - 1; nr_wiersza++)
                    dopelnienie[nr_wiersza][nr_kolumny_dop] = macierz[nr_wiersza + 1][nr_kolumny_mac];
            }
            // det = ?(aij * (-1)^i+j * Aij)
            det += (macierz[0][nr_kolumny] * (long double) pow(-1.0, 1.0 + nr_kolumny + 1.0) *
                    metodaLaplace(dopelnienie, stopien - 1, count - 1));
        }
    } else {
        printf("parallel %d", omp_get_num_threads());
        // Wyznacznik stopnia co najwyżej 2. zostanie obliczony z definicji
        if (stopien <= 2)
            return macierz[0][0] * (stopien > 1 ? macierz[1][1] : 1) -
                   (stopien > 1 ? macierz[0][1] * macierz[1][0] : 0);

#pragma omp parallel private(nr_wiersza, nr_kolumny, nr_kolumny_dop, nr_kolumny_mac, dopelnienie,stopien, macierz)
        dopelnienie = nowaMacierz(stopien - 1);

#pragma omp parallel for
        for (nr_kolumny = 0; nr_kolumny < stopien; nr_kolumny++) {
            for (nr_kolumny_dop = 0, nr_kolumny_mac = 0;
                 nr_kolumny_dop < stopien - 1; nr_kolumny_dop++, nr_kolumny_mac++) {
                nr_kolumny_mac += (nr_kolumny_mac == nr_kolumny ? 1
                                                                : 0); 
                for (nr_wiersza = 0; nr_wiersza < stopien - 1; nr_wiersza++)
                    dopelnienie[nr_wiersza][nr_kolumny_dop] = macierz[nr_wiersza + 1][nr_kolumny_mac];
            }
            // det = ?(aij * (-1)^i+j * Aij)
            det += (macierz[0][nr_kolumny] * (long double) pow(-1.0, 1.0 + nr_kolumny + 1.0) *
                    metodaLaplace(dopelnienie, stopien - 1, count - 1));
        }
    }

    usunMacierz(dopelnienie);
    return det;
}

核心错误点分析 #

1. 错误的private变量声明
   你将stopien和macierz声明为private变量，但这两个是函数的核心输入参数，所有线程都需要访问同一个矩阵和矩阵大小。设为private后，每个线程的macierz会被初始化为未定义的垃圾值，访问时直接触发段错误。

2. 并行区域的内存冲突
   你在单独的#pragma omp parallel区域内分配dopelnienie，导致每个线程都创建一个子矩阵副本，但后续的并行for循环会让所有线程同时写入同一个主线程的dopelnienie，引发严重的数据竞争和内存访问错误。

3. 共享变量的竞争条件
   多个线程同时执行det += ...，而+=并非原子操作，不仅会导致计算结果错误，还可能引发内存访问异常。

4. 递归并行的过度嵌套
   每次递归调用都会创建新的并行区域，导致线程数量爆炸，超出系统资源限制，同时加剧内存冲突风险。

修正后的代码 #

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

long double **nowaMacierz(int stopien);
void usunMacierz(long double **macierz);
long double metodaLaplace(long double **macierz, int stopien);

int main() {
    long double **macierz;
    int stopien;
    int nr_wiersza, nr_kolumny;

    printf("Proszę podać stopień macierzy n=");
    fscanf(stdin, "%d", &stopien);
    macierz = nowaMacierz(stopien);
    for (nr_wiersza = 0; nr_wiersza < stopien; nr_wiersza++)
        for (nr_kolumny = 0; nr_kolumny < stopien; nr_kolumny++) {
            printf("A[%d,%d]=", nr_wiersza, nr_kolumny);
            fscanf(stdin, "%Lf", &macierz[nr_wiersza][nr_kolumny]);
        }
    printf("det(A) = %Lf\n", metodaLaplace(macierz, stopien));
    usunMacierz(macierz);
    getchar();
    return 0;
}

long double **nowaMacierz(int stopien) {
    long double **macierz;
    int nr_wiersza;

    macierz = (long double **)calloc(stopien, sizeof(long double *));
    *macierz = (long double *)calloc(stopien * stopien, sizeof(long double));

    for (nr_wiersza = 1; nr_wiersza < stopien; nr_wiersza++)
        *(macierz + nr_wiersza) = *(macierz + nr_wiersza - 1) + stopien;
    return macierz;
}

void usunMacierz(long double **macierz) {
    free(*macierz);
    free(macierz);
}

long double metodaLaplace(long double **macierz, int stopien) {
    // 设置并行阈值，小矩阵用串行计算
    const int PARALLEL_THRESHOLD = 3;
    long double det = 0.0L;

    if (stopien <= 2) {
        return macierz[0][0] * (stopien > 1 ? macierz[1][1] : 1.0L) -
               (stopien > 1 ? macierz[0][1] * macierz[1][0] : 0.0L);
    }

    if (stopien > PARALLEL_THRESHOLD) {
        // 用reduction处理累加，避免竞争条件
        #pragma omp parallel for reduction(+:det) \
            private(nr_wiersza, nr_kolumny, nr_kolumny_dop, nr_kolumny_mac, dopelnienie)
        for (int nr_kolumny = 0; nr_kolumny < stopien; nr_kolumny++) {
            long double **dopelnienie = nowaMacierz(stopien - 1);
            int nr_wiersza, nr_kolumny_dop, nr_kolumny_mac;

            for (nr_kolumny_dop = 0, nr_kolumny_mac = 0; nr_kolumny_dop < stopien - 1; nr_kolumny_dop++, nr_kolumny_mac++) {
                if (nr_kolumny_mac == nr_kolumny) {
                    nr_kolumny_mac++;
                }
                for (nr_wiersza = 0; nr_wiersza < stopien - 1; nr_wiersza++) {
                    dopelnienie[nr_wiersza][nr_kolumny_dop] = macierz[nr_wiersza + 1][nr_kolumny_mac];
                }
            }

            // 简化符号计算，避免浮点pow的误差
            long double sign = ((nr_kolumny % 2) == 0) ? 1.0L : -1.0L;
            det += macierz[0][nr_kolumny] * sign * metodaLaplace(dopelnienie, stopien - 1);

            usunMacierz(dopelnienie);
        }
    } else {
        // 小矩阵串行计算
        long double **dopelnienie = nowaMacierz(stopien - 1);
        for (int nr_kolumny = 0; nr_kolumny < stopien; nr_kolumny++) {
            int nr_wiersza, nr_kolumny_dop, nr_kolumny_mac;

            for (nr_kolumny_dop = 0, nr_kolumny_mac = 0; nr_kolumny_dop < stopien - 1; nr_kolumny_dop++, nr_kolumny_mac++) {
                if (nr_kolumny_mac == nr_kolumny) {
                    nr_kolumny_mac++;
                }
                for (nr_wiersza = 0; nr_wiersza < stopien - 1; nr_wiersza++) {
                    dopelnienie[nr_wiersza][nr_kolumny_dop] = macierz[nr_wiersza + 1][nr_kolumny_mac];
                }
            }

            long double sign = ((nr_kolumny % 2) == 0) ? 1.0L : -1.0L;
            det += macierz[0][nr_kolumny] * sign * metodaLaplace(dopelnienie, stopien - 1);
        }
        usunMacierz(dopelnienie);
    }

    return det;
}

修正说明 #

• 移除了错误的stopien和macierz的private声明，保留核心参数的共享性。
• 使用reduction(+:det)处理行列式累加，避免线程竞争。
• 每个并行迭代独立创建和释放子矩阵，消除内存冲突。
• 设置并行阈值，小矩阵用串行计算，避免递归过度并行的资源浪费。
• 简化符号计算逻辑，替换浮点pow为整数取模，提升效率并避免误差。

内容的提问来源于stack exchange，提问作者K' CHU