Hello,

I have implemented pthreads for the sequential Matrix Inverse Program. The program will do the division step first and then the elimination is done for all rows except one row. I tried implementing pthreads for the elimination part only. But for every row iteration, I create and join (destroying) threads so that the execution time is the same as the sequential program. But I have an idea to reduce the time by creating a certain number of threads first and then using those threads for both the division step and the elimination step for all the rows, and inserting barriers in between the steps. Finally destroying the threads at once, so that the thread creation overhead will be reduced and execution time will be reduced. I don't get any logic to implement this scenario. I will be thankful if any suggestions or help.

The pthread implemented code is as follows:

#include #include #include #include #include

#define MAX_SIZE 4096

typedef double matrix[MAX_SIZE][MAX_SIZE];

int N; /* matrix size */ int maxnum; /* max number of element*/ char* Init; /* matrix init type */ int PRINT; /* print switch */ matrix A; /* matrix A */ matrix I = {0.0}; /* The A inverse matrix, which will be initialized to the identity matrix */

struct intervals* args; int nothreads = 2, load; int* data = NULL; pthread_t* threads;

/* forward declarations */ void find_inverse(void); void Init_Matrix(void); void Print_Matrix(matrix M, char name[]); void Init_Default(void); int Read_Options(int, char**); void* find_inverse_par(void* data);

int main(int argc, char** argv) { printf("Matrix Inverse "); int i, timestart, timeend, iter;

Init_Default(); /* Init default values */ Read_Options(argc, argv); /* Read arguments */ Init_Matrix(); /* Init the matrix */ threads = (pthread_t*)malloc(nothreads * sizeof(pthread_t)); find_inverse();

if (PRINT == 1) { //Print_Matrix(A, "End: Input"); Print_Matrix(I, "Inversed"); } }

void find_inverse(void) { int i, j, k;

int row, col, p; // 'p' stands for pivot (numbered from 0 to N-1) double pivalue; // pivot value

/* Bringing the matrix A to the identity form */ for (p = 0; p < N; p++) { /* Outer loop */ pivalue = A[p][p]; for (col = 0; col < N; col++) { A[p][col] = A[p][col] / pivalue; /* Division step on A */ I[p][col] = I[p][col] / pivalue; /* Division step on I */ } assert(A[p][p] == 1.0);

load = (N/nothreads); //printf("load %d for k %d and n %d ",load,k,N); for (int i =0; i

data[0]=p; //data[1]=i; //if(i==0){ data[1]=load*i; data[2]=load*i+load; /*else{ data[1]=load*i+1; data[2]=load*i+1+load+1;}*/ if(i==nothreads-1){ int diff=N-data[2]; data[2]+=diff; } //printf("The Thread %d load is from %d to %d",i,data[1], data[2]-1); pthread_create(&(threads[i]), NULL, find_inverse_par, (void *)data); } for (unsigned int id = 0; id < nothreads; id++) { pthread_join(threads[id], NULL ); } } }

void* find_inverse_par(void* data) { int* dat=(int*) data; int p=dat[0]; int start=dat[1]; int end=dat[2];

double multiplier; int row,col; for (row = start; row < end; row++) { multiplier = A[row][p]; if (row != p) // Perform elimination on all except the current pivot row { for (col = 0; col < N; col++) { A[row][col] = A[row][col] - A[p][col] * multiplier; /* Elimination step on A */ I[row][col] = I[row][col] - I[p][col] * multiplier; /* Elimination step on I */ } assert(A[row][p] == 0.0); } }

}

void Init_Matrix() { int row, col;

// Set the diagonal elements of the inverse matrix to 1.0 // So that you get an identity matrix to begin with for (row = 0; row < N; row++) { for (col = 0; col < N; col++) { if (row == col) I[row][col] = 1.0; } }

printf(" size = %dx%d ", N, N); printf(" maxnum = %d ", maxnum); printf("Init = %s ", Init); printf("Initializing matrix...");

if (strcmp(Init, "rand") == 0) { for (row = 0; row < N; row++) { for (col = 0; col < N; col++) { if (row == col) /* diagonal dominance */ A[row][col] = (double)(rand() % maxnum) + 5.0; else A[row][col] = (double)(rand() % maxnum) + 1.0; } } } if (strcmp(Init, "fast") == 0) { for (row = 0; row < N; row++) { for (col = 0; col < N; col++) { if (row == col) /* diagonal dominance */ A[row][col] = 5.0; else A[row][col] = 2.0; } } }

printf("done "); if (PRINT == 1) { //Print_Matrix(A, "Begin: Input"); //Print_Matrix(I, "Begin: Inverse"); } }

void Print_Matrix(matrix M, char name[]) { int row, col;

printf("%s Matrix: ", name); for (row = 0; row < N; row++) { for (col = 0; col < N; col++) printf(" %5.2f", M[row][col]); printf(" "); } printf(" "); }

void Init_Default() { N = 5; Init = "fast"; maxnum = 15.0; PRINT = 0; }

int Read_Options(int argc, char** argv) { char* prog;

prog = *argv; while (++argv, --argc > 0) if (**argv == '-') switch (*++ * argv) { case 'n': --argc; N = atoi(*++argv); break; case 'h': printf(" HELP: try matinv -u "); exit(0); break; case 'u': printf(" Usage: matinv [-n problemsize] "); printf(" [-D] show default values "); printf(" [-h] help "); printf(" [-I init_type] fast/rand "); printf(" [-m maxnum] max random no "); printf(" [-P print_switch] 0/1 "); exit(0); break; case 'D': printf(" Default: n = %d ", N); printf(" Init = rand"); printf(" maxnum = 5 "); printf(" P = 0 "); exit(0); break; case 'I': --argc; Init = *++argv; break; case 'm': --argc; maxnum = atoi(*++argv); break; case 'P': --argc; PRINT = atoi(*++argv); break; default: printf("%s: ignored option: -%s ", prog, *argv); printf("HELP: try %s -u ", prog); break; } }