小贝子编程

C-使OpenMP程序与Pthreads一起使用，Segfault错误



我有一个书面程序，该程序在C中执行高斯消除并返回矩阵的L2 Norm。该程序被称为./exec n k，其中n是N矩阵的A N大小，K是用于执行程序的线程数（Max 4）。我为n 1矩阵分配空间，因为具有增强矩阵是高斯消除的一部分。

它在OpenMP中完美工作。如下面的代码所示，我只有1个平行。我现在的目标是使循环使用PTHreads而不是OpenMP同时运行。我制作了将循环并行的循环，并创建pthreads来处理它。我的猜测是，pthreads并非每个人都在做循环的不同部分（基本上是j的迭代不同），而是4个pthreads只是在运行整个循环。我运行了像./gauss 30 4这样的程序，有时可以工作，有时是segfaults，尽管当它起作用时，L2标准不是0（如果程序完美地工作，则L2将返回0），因此显然，线程部分有些事情。当我通过GDB运行它时，由于某种原因，它会在循环中以循环播放，但是相同的循环在OpenMP中完美运行...有人可以帮助我

GDB

https://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 threadNumn" );
            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: %fn", l2Norm);
        gettimeofday(&end, NULL);
        delta = ((end.tv_sec  - start.tv_sec) * 1000000u + 
             end.tv_usec - start.tv_usec) / 1.e6;
        printf("end time: %fn", 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 threadNumn" );
        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 problemn");
            exit(1);
            }
        }
    while(threadcount != 0){
        if(pthread_join (workerThreads[threadcount-1], &retval ) != 0){
        fprintf(stderr, "Error: consumer create problemn");
        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: %fn", l2Norm);
    gettimeofday(&end, NULL);
    delta = ((end.tv_sec  - start.tv_sec) * 1000000u + 
         end.tv_usec - start.tv_usec) / 1.e6;
    printf("end time: %fn", delta);
}
        }

请注意，每个线程中应将j , z声明为本地（私有）变量。
在OpenMP代码中，您在第100行中关闭了For Loop的支撑：

    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 threadNumn" );
        exit(1);
    }
}
void *parallelstuff(void *args){
    threader temp = *((threader *)args);
    int i, matrixSize;
    i = temp.counter;
    matrixSize = temp.matrixl;
    //printf("matrixSize=%d counter=%dn" , 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=%dn" , 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 problemn");
            exit(1);
            }
        }
    while(threadcount != 0){
        if(pthread_join (workerThreads[threadcount-1], &retval ) != 0){
        fprintf(stderr, "Error: consumer create problemn");
        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: %fn", l2Norm);
    gettimeofday(&end, NULL);
    delta = ((end.tv_sec  - start.tv_sec) * 1000000u + 
         end.tv_usec - start.tv_usec) / 1.e6;
    printf("end time: %fn", delta);
}

书面的两个代码不是等效的。观察OpenMP代码：

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

OpenMP中的parallel for组合构建体是一个工作共享结构，即它在团队中的线程之间在以下循环上分发迭代。给定schedule(static)子句，迭代空间分为#threads块，每个块分配给另一个线程。

您的pthreads代码没有共享工作：

i = temp.counter;
matrixSize = temp.matrixl;
...
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];
    }
}

鉴于相同的stuff对象传递给所有线程，它们都接收到i和matrixSize的相同值，并在整个迭代空间上循环，因此结果错误。

您要做的就是模拟#pragma omp for schedule(static)的作用，即使每个线程仅执行某些matrixSize - (i+1) + 1迭代。您应该传递每个线程一个唯一的数据对象，其中包含启动和结束迭代：

typedef struct threader {
    int start;
    int end;
    int i;
    int matrixSize;
} threader;
...
void *parallelstuff(void *args){
    threader *temp = (threader *)args;
    int start, end, i, matrixSize;
    start = temp->start;
    end = temp->end;
    i = temp->i;
    matrixSize = temp->matrixSize;
    double temp2;
    int j , z; //houssam
    for(j = start + 1; j<end; 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];
        }
    }
}
...
threader stuff[threadNum];
//MAKE EACH THREAD DO SOMETHING DIFF
// parallelstuff(int i, int matrixSize, double **a){
for(threadcount = 0; threadcount < threadNum; threadcount++){
    stuff[threadcount].start = i + threadcount*(matrixSize / threadNum);
    stuff[threadcount].end = i + (threadcount+1)*(matrixSize / threadNum);
    stuff[threadcount].i = i;
    stuff[threadcount].matrixSize = matrixSize;
    if(pthread_create (&workerThreads[threadcount], NULL, parallelstuff, (void *) &stuff ) != 0){
        fprintf(stderr, "Error: consumer create problemn");
        exit(1);
    }
}

从理论上讲，您还可以让每个线程知道那里有多少个线程，并让其计算迭代范围本身，但是这对于Pthreads API缺乏相当于omp_get_thread_num()的事实而变得复杂。有一个高级技巧，它采用了对齐的内存分配，并且在传递的指针的最后位中编码的数字线程ID。

相关内容

最新更新


热门标签：

javascript python java c# php android html jquery c++ css ios sql mysql arrays asp.net json python-3.x ruby-on-rails .net sql-server django objective-c excel regex ruby linux ajax iphone xml vba spring asp.net-mvc database wordpress string postgresql wpf windows xcode bash git oracle list vb.net multithreading eclipse algorithm macos powershell visual-studio image forms numpy scala function api selenium