c - 使 OpenMP 程序与 Pthreads 一起工作，出现段错误

Question

我编写了一个程序，用 C 语言执行高斯消元法并返回矩阵的 L2 范数。该程序的调用方式类似于 ./exec n k，其中 n 是 n × n 矩阵的大小，k 是将用于执行该程序的线程数(最大4).我为 n x n+1 矩阵分配空间，因为增广矩阵是高斯消除的一部分。

它在 OpenMP 中完美运行。正如下面的代码所示，我只有 1 个并行。我现在的目标是使用 Pthreads 而不是 OpenMP 使并行 for 循环同时运行。我将 for 循环并行化为一个单独的函数，并创建 pthread 来处理它。我的猜测是，pthreads 并不是各自执行循环的不同部分(基本上是 j 的不同迭代)，而是 4 个 pthreads 只是运行整个循环。我运行像 ./gauss 30 4 这样的程序，它有时可以工作，有时会出现段错误，尽管当它确实工作时，L2 范数不是 0(如果程序完美工作，L2 将返回 0)，所以有些东西是显然与线程部分有关。当我通过 GDB 运行它时，由于某种原因，它在循环中出现段错误，但同一个循环在 OpenMP 中运行完美......有人可以帮助我吗

GDB

http://i.stack.imgur.com/V99yt.png

OpenMP 代码:

 #include <stdlib.h>
    #include <stdio.h>
    #include <math.h>
    #include <omp.h>
    #include <time.h>
    #include <sys/time.h>
    //globals
    double **a, *vect, *bvect, scalar, ratio, sum, delta, *temp;
    int i,j,k,ptr, z;
    int y,z;
    int bvectcount = 0;
    struct timeval start, end;
    // a is matrix, b is vector, x is the solution vector, and n is the size 
    double L2(double **a, double *bvect, double *vect, int matrixSize) {
        double sum;
        double res[matrixSize];
        int i, j;
        for (i=0; i < matrixSize; i++) {
            sum = (double) 0;
            for (j=0; j < matrixSize; j++) {
                sum += a[i][j] * vect[j];
            }
            res[i] = sum;
        }
        for (i=0; i < matrixSize; i++) {
            res[i] -= vect[i];
        }
        double sum_squares = (double) 0;
        for (i=0; i < matrixSize; i++) {
            sum_squares += res[i] * res[i];
        }
        return sqrt(sum_squares);
    }
    int checkargs(int argc, char* argv[]){
        if(argc != 3){
            fprintf(stderr, "Error: Usage is size threadNum\n" );
            exit(1);
        }
    }
    int main(int argc, char* argv[]){
        //check for args
        checkargs(argc, argv);
        int matrixSize = atoi(argv[1]);
        int threadNum = atoi(argv[2]);
        int chunk = matrixSize/threadNum;
        //memory allocation
        a = (double**)malloc(matrixSize*sizeof(double*));
        for(i = 0; i < matrixSize ; i++)
            a[i] = (double*)malloc(matrixSize*sizeof(double) * matrixSize);
        vect = (double*)malloc(matrixSize*sizeof(double));
        bvect = (double*)malloc(matrixSize*sizeof(double));
        temp = (double*)malloc(matrixSize*sizeof(double));
        for(i = 0; i < matrixSize; ++i){
            for(j = 0; j < matrixSize + 1; ++j){
                a[i][j] = drand48(); 
            }
        }   
        j = 0;
        j += matrixSize;
        for(i = 0; i < matrixSize; ++i){
            bvect[i] = a[i][j];
        }
        //generation of scalar matrix (diagonal vector)
        gettimeofday(&start, NULL);
        for(i=0; i<matrixSize; i++){
            scalar = a[i][i];
            //initialization of p to travel throughout matrix
            ptr = i;
            //find largest number in column and row number of it
            for(k = i+1; k < matrixSize; k++){
            if(fabs(scalar) < fabs(a[k][i])){
                //k is row of scalar, while
               scalar = a[k][i];
               ptr = k;
            }
        }
        //swaping the elements of diagonal row and row containing largest no
        for(j = 0; j <= matrixSize; j++){
            temp[0] = a[i][j];
            a[i][j]= a[ptr][j];
            a[ptr][j] = temp[0];
        }
        //calculating triangular matrix
        //threading needs to be done HERE
        ratio = a[i][i];
        for(k = 0; k < matrixSize + 1; k++){
            a[i][k] = a[i][k] / ratio;
        }
        double temp2;
        #pragma omp parallel default(none) num_threads(threadNum)  shared(a,i,matrixSize,vect) private(j,z,ratio,temp2)
        {
        #pragma omp for schedule(static)
        for(j = i + 1; j<matrixSize; j++){
            temp2 = a[j][i]/a[i][i];
            for(z = 0; z<matrixSize + 1; z++){
                a[j][z] = a[j][z] - temp2 * a[i][z];
                }
            }
        }
    }
        //backward substitution method
        for(i=matrixSize-1; i >=0; i--){
            for(k = i; k > 0; k--){
                a[k-1][matrixSize] -= a[k-1][i] * a[i][matrixSize];
                a[k-1][i] -= a[k-1][i] * a[i][i];
            }     
        }
        for(i = 0; i < matrixSize; ++i){
            vect[i] = a[i][matrixSize];
        }
        double l2Norm;
        l2Norm = L2(a, bvect, vect, matrixSize);
        printf("THIS IS L2 NORM: %f\n", l2Norm);
        gettimeofday(&end, NULL);
        delta = ((end.tv_sec  - start.tv_sec) * 1000000u + 
             end.tv_usec - start.tv_usec) / 1.e6;
        printf("end time: %f\n", delta);
    }

Pthreads 代码:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <omp.h>
#include <time.h>
#include <sys/time.h>
#include <pthread.h>

//globals
double **a, *vect, *bvect, scalar, ratio, sum, delta, *temp;
int i,j,k,ptr, z;
int y,z;
int bvectcount = 0;
int threadcount;
pthread_t workerThreads[4];
typedef struct threader {
    int counter;
    int matrixl;
} threader;
struct timeval start, end;
    void *retval;

int checkargs(int argc, char* argv[]);
// a is matrix, b is vector, x is the solution vector, and n is the size 
double L2(double **a, double *bvect, double *vect, int matrixSize) {
    double sum;
    double res[matrixSize];
    int i, j;
    for (i=0; i < matrixSize; i++) {
        sum = (double) 0;
        for (j=0; j < matrixSize; j++) {
            sum += a[i][j] * vect[j];
        }
        res[i] = sum;
    }
    for (i=0; i < matrixSize; i++) {
        res[i] -= vect[i];
    }
    double squaresum = (double) 0;
    for (i=0; i < matrixSize; i++) {
        squaresum += res[i] * res[i];
    }
    return sqrt(squaresum);
}
int checkargs(int argc, char* argv[]){
    if(argc != 3){
        fprintf(stderr, "Error: Usage is size threadNum\n" );
        exit(1);
    }
}

void *parallelstuff(void *args){
    threader temp = *((threader *)args);
    int i, matrixSize;
    i = temp.counter;
    matrixSize = temp.matrixl;
    double temp2;
    for(j = i + 1; j<matrixSize; j++){
        temp2 = a[j][i]/a[i][i];
        for(z = 0; z<matrixSize + 1; z++){
            a[j][z] = a[j][z] - temp2 * a[i][z];
        }
    }
}
int main(int argc, char* argv[]){
    //check for args
    checkargs(argc, argv);
    int matrixSize = atoi(argv[1]);
    int threadNum = atoi(argv[2]);
    //memory allocation
    a = (double**)malloc(matrixSize*sizeof(double*));
    for(i = 0; i < matrixSize ; i++)
        a[i] = (double*)malloc(matrixSize*sizeof(double) * matrixSize);
    vect = (double*)malloc(matrixSize*sizeof(double));
    bvect = (double*)malloc(matrixSize*sizeof(double));
    temp = (double*)malloc(matrixSize*sizeof(double));
    for(i = 0; i < matrixSize; ++i){
        for(j = 0; j < matrixSize + 1; ++j){
            a[i][j] = drand48(); 
        }
    }
    j = 0;
    j += matrixSize;
    for(i = 0; i < matrixSize; ++i){
        bvect[i] = a[i][j];
    }
//generation of scalar matrix (diagonal vector)
    gettimeofday(&start, NULL);
    for(i=0; i<matrixSize; i++){
        scalar = a[i][i];
    //initialization of p to travel throughout matrix
        ptr = i;
    //find largest number in column and row number of it
        for(k = i+1; k < matrixSize; k++){
        if(fabs(scalar) < fabs(a[k][i])){
            //k is row of scalar, while
           scalar = a[k][i];
           ptr = k;
        }
    }
    //swaping the elements of diagonal row and row containing largest no
    for(j = 0; j <= matrixSize; j++)
    {
        temp[0] = a[i][j];
        a[i][j]= a[ptr][j];
        a[ptr][j] = temp[0];
    }
    ratio = a[i][i];
    for(k = 0; k < matrixSize + 1; k++){
        a[i][k] = a[i][k] / ratio;
    }   
    threader stuff;
    stuff.counter = i;
    stuff.matrixl = matrixSize;
    //MAKE EACH THREAD DO SOMETHING DIFF
    // parallelstuff(int i, int matrixSize, double **a){
    for(threadcount = 0; threadcount < threadNum; threadcount++){
        if(pthread_create (&workerThreads[threadcount], NULL, parallelstuff, (void *) &stuff ) != 0){
            fprintf(stderr, "Error: consumer create problem\n");
            exit(1);
            }
        }
    while(threadcount != 0){
        if(pthread_join (workerThreads[threadcount-1], &retval ) != 0){
        fprintf(stderr, "Error: consumer create problem\n");
        exit(1);
        }
        threadcount--;
    }

    //create matrix of n size 
    //backward substitution method
    for(i=matrixSize-1; i >=0; i--){
        for(k = i; k > 0; k--){
            a[k-1][matrixSize] -= a[k-1][i] * a[i][matrixSize];
            a[k-1][i] -= a[k-1][i] * a[i][i];
        }  
    }
    for(i = 0; i < matrixSize; ++i){
        vect[i] = a[i][matrixSize];
        }       
    double l2Norm;
    l2Norm = L2(a, bvect, vect, matrixSize);
    printf("THIS IS L2 NORM: %f\n", l2Norm);
    gettimeofday(&end, NULL);
    delta = ((end.tv_sec  - start.tv_sec) * 1000000u + 
         end.tv_usec - start.tv_usec) / 1.e6;
    printf("end time: %f\n", delta);
}
        }

Answer 1

请注意，j , z 应在每个线程中声明为本地(私有(private))变量。
在 OpenMP 代码中，您在第 100 行关闭了 for 循环的大括号:

    gettimeofday(&start, NULL);
    for(i=0; i<matrixSize; i++){
        scalar = a[i][i];
        //initialization of p to travel throughout matrix
        .......
        ......
        .....

} //line 100

但是在 pthreads 代码中，您在第 149 行关闭了它，因此完整代码:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <omp.h>
#include <time.h>
#include <sys/time.h>
#include <pthread.h>

//globals
double **a, *vect, *bvect, scalar, ratio, sum, delta, *temp;
int i,j,k,ptr, z;
int y; //z?
int bvectcount = 0;
int threadcount;
pthread_t workerThreads[4];
typedef struct threader {
    int counter;
    int matrixl;
} threader;
struct timeval start, end;
    void *retval;

int checkargs(int argc, char* argv[]);
// a is matrix, b is vector, x is the solution vector, and n is the size 
double L2(double **a, double *bvect, double *vect, int matrixSize) {
    double sum;
    double res[matrixSize];
    int i, j;
    for (i=0; i < matrixSize; i++) {
        sum = (double) 0;
        for (j=0; j < matrixSize; j++) {
            sum += a[i][j] * vect[j];
        }
        res[i] = sum;
    }
    for (i=0; i < matrixSize; i++) {
        res[i] -= vect[i];
    }
    double squaresum = (double) 0;
    for (i=0; i < matrixSize; i++) {
        squaresum += res[i] * res[i];
    }
    return sqrt(squaresum);
}
int checkargs(int argc, char* argv[]){
    if(argc != 3){
        fprintf(stderr, "Error: Usage is size threadNum\n" );
        exit(1);
    }
}

void *parallelstuff(void *args){
    threader temp = *((threader *)args);
    int i, matrixSize;
    i = temp.counter;
    matrixSize = temp.matrixl;
    //printf("matrixSize=%d counter=%d\n" , matrixSize ,temp.counter );
    double temp2;
    int j , z; //houssam
    for(j = i + 1; j<matrixSize; j++){
        temp2 = a[j][i]/a[i][i];
        for(z = 0; z<matrixSize + 1; z++){
            a[j][z] = a[j][z] - temp2 * a[i][z];
        }
    }


}
int main(int argc, char* argv[]){
    //check for args
    checkargs(argc, argv);
    int matrixSize = atoi(argv[1]);
    int threadNum = atoi(argv[2]);
    //memory allocation
    a = (double**)malloc(matrixSize*sizeof(double*));
    for(i = 0; i < matrixSize ; i++)
        a[i] = (double*)malloc(matrixSize*sizeof(double) * matrixSize);
    vect = (double*)malloc(matrixSize*sizeof(double));
    bvect = (double*)malloc(matrixSize*sizeof(double));
    temp = (double*)malloc(matrixSize*sizeof(double));
    for(i = 0; i < matrixSize; ++i){
        for(j = 0; j < matrixSize + 1; ++j){
            a[i][j] = drand48();
        }
    }
    j = 0;
    j += matrixSize;
    for(i = 0; i < matrixSize; ++i){
        bvect[i] = a[i][j];
    }
//generation of scalar matrix (diagonal vector)
    gettimeofday(&start, NULL);
    for(i=0; i<matrixSize; i++){
        scalar = a[i][i];
    //initialization of p to travel throughout matrix
        ptr = i;
    //find largest number in column and row number of it
        for(k = i+1; k < matrixSize; k++){
        if(fabs(scalar) < fabs(a[k][i])){
            //k is row of scalar, while
           scalar = a[k][i];
           ptr = k;
        }
    }
    //swaping the elements of diagonal row and row containing largest no
    for(j = 0; j <= matrixSize; j++)
    {
        temp[0] = a[i][j];
        a[i][j]= a[ptr][j];
        a[ptr][j] = temp[0];
    }
    ratio = a[i][i];
    for(k = 0; k < matrixSize + 1; k++){
        a[i][k] = a[i][k] / ratio;
    }   

    threader  stuff;
    stuff.counter = i;
    stuff.matrixl = matrixSize;
    //printf("i=%d\n" , i);
    //MAKE EACH THREAD DO SOMETHING DIFF
    // parallelstuff(int i, int matrixSize, double **a){
    for(threadcount = 0; threadcount < threadNum; threadcount++){
        if(pthread_create (&workerThreads[threadcount], NULL, parallelstuff, (void *) &stuff ) != 0){
            fprintf(stderr, "Error: consumer create problem\n");
            exit(1);
            }
        }

    while(threadcount != 0){
        if(pthread_join (workerThreads[threadcount-1], &retval ) != 0){
        fprintf(stderr, "Error: consumer create problem\n");
        exit(1);
        }
        threadcount--;
    }
}
    //create matrix of n size 
    //backward substitution method
    for(i=matrixSize-1; i >=0; i--){
        for(k = i; k > 0; k--){
            a[k-1][matrixSize] -= a[k-1][i] * a[i][matrixSize];
            a[k-1][i] -= a[k-1][i] * a[i][i];
        }  
    }
    for(i = 0; i < matrixSize; ++i){
        vect[i] = a[i][matrixSize];
        }       
    double l2Norm;
    l2Norm = L2(a, bvect, vect, matrixSize);
    printf("THIS IS L2 NORM: %f\n", l2Norm);
    gettimeofday(&end, NULL);
    delta = ((end.tv_sec  - start.tv_sec) * 1000000u + 
         end.tv_usec - start.tv_usec) / 1.e6;
    printf("end time: %f\n", delta);

}