我编写了一个程序，用 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);
}
}
最佳答案
请注意，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);
}
关于c - 使 OpenMP 程序与 Pthreads 一起工作，出现段错误，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/30017462/