/**************************************************************** * * * gjmpi.c * * Gauss-Jordan method for solving systems of linear equations * * MPI-Version (row-oriented) * * Joerg Meyer University of Nebraska at Omaha * * 07/15/94 Computer Science Department * * 07/02/95 Bug fix * * * ****************************************************************/ #include #include #include #include #include "mpi.h" #define N 50 /* size of linear system */ #define EPSILON 1e-20 #define FALSE 0 #define TRUE 1 #ifndef MICROSEC #define MICROSEC 1000000 #endif typedef double mrow_t[N+1]; /* used to reference dynamically created */ typedef mrow_t *p_mrow_t; /* two-dimensional array */ /**************************************************************** * * * initialize_system (a) - randomly generates system of * * linear equations * * parameter: a - pointer to matrix * * * ****************************************************************/ void initialize_system (p_mrow_t a) { int i, j, seed; /* time (&seed); srand ((unsigned int)seed);*/ srand (0); for (i = 0; i < N; i++) { a[i][N] = 0.0; for (j = 0; j < N; j++) { a[i][j] = rand (); a[i][N] += (j - 50) * a[i][j]; } } } /**************************************************************** * * * print_results (a) - prints solution of system of * * linear equations * * parameter: a - pointer to matrix * * * ****************************************************************/ void print_results (p_mrow_t a) { int i; char s[80]; for (i = 0; i < N; i++) { printf ("x[%ld] = %10.6lf\n", i, a[i][N]); } } /**************************************************************** * * * gauss_jordan_elimination (m) - * * row-oriented Gauss-Jordan algorithm implementation * * solves system of linear equations in parallel * * using MPI calls, all tasks have to call this * * function, only task 0 actually delivers pointer * * to matrix (matrix overwritten) * * parameter: m - pointer to matrix * * result : solution vector in m[i][N] * * * ****************************************************************/ int gauss_jordan_elimination (p_mrow_t m) { int myrank, ranksize; MPI_Status status; int startrow, lastrow, nrow, sr, lr; p_mrow_t mp; double *pivotrow; int *pivotp, *marked; int i, j, k, picked; double tmp; struct { double val; int rank; } in, out; /* assume MPI_Init has already been called */ MPI_Comm_rank (MPI_COMM_WORLD, &myrank); MPI_Comm_size (MPI_COMM_WORLD, &ranksize); /* rows of matrix I have to process */ startrow = (myrank * N) / ranksize; lastrow = ((myrank + 1) * N) / ranksize; nrow = lastrow - startrow; /* dynamically allocate data structures */ /* copy of my portion of matrix */ mp = (p_mrow_t) malloc (nrow * (N+1) * sizeof (double)); /* row to be used next for reduction of matrix */ pivotrow = (double *) malloc ((N+1) * sizeof (double)); /* stores which column was reduced using this row */ pivotp = (int *) malloc (N * sizeof (int)); /* stores if row was used for reduction of column */ marked = (int *) malloc (nrow * sizeof (int)); if (myrank == 0) /* process 0 */ { /* distribute data among processes */ for (i = 1; i < ranksize; i++) { sr = (i * N) / ranksize; lr = ((i + 1) * N) / ranksize; MPI_Send (m[sr], (lr-sr) * (N+1), MPI_DOUBLE, i, 1234, MPI_COMM_WORLD); } for (i = startrow; i < lastrow; i++) for (j = 0; j < N + 1; j++) mp[i][j] = m[i][j]; } else { /* receive portion of matrix */ MPI_Recv (mp, nrow * (N+1), MPI_DOUBLE, 0, 1234, MPI_COMM_WORLD, &status); } for (i = 0; i < nrow; i++) marked[i] = 0; for (i = 0; i < N; i++) /* for each column */ { /* WinMPI extension MPI_Win_yield ();*/ /* find local maximum in column */ tmp = 0.0; for (j = 0; j < nrow; j++) { if (!marked[j] && (fabs (mp[j][i]) > tmp)) { tmp = fabs (mp[j][i]); picked = j; } } /* find global maximum in column and which process stores it */ in.val = tmp; in.rank = myrank; MPI_Allreduce (&in, &out, 1, MPI_DOUBLE_INT, MPI_MAXLOC, MPI_COMM_WORLD); if (out.rank == myrank) /* I am process with the global maximum */ { marked[picked] = 1; pivotp[picked] = i; for (j = 0; j < N + 1; j++) pivotrow[j] = mp[picked][j]; } /* process with global maximum broadcasts entire row (pivot row) */ MPI_Bcast (pivotrow, N + 1, MPI_DOUBLE, out.rank, MPI_COMM_WORLD); if (fabs (pivotrow[i]) < EPSILON) /* no solution */ { if (myrank == 0) printf ("Exits on iteration %d\n", i); return (FALSE); } /* reduce all rows using broadcasted pivot row */ for (j = 0; j < nrow; j++) if (!(marked[j] && pivotp[j] == i)) { tmp = mp[j][i] / pivotrow[i]; for (k = i; k < N + 1; k++) mp[j][k] -= pivotrow[k] * tmp; } } if (myrank == 0) /* process 0 */ /* compute local part of solution and collect results of other processes */ { for (i = 0; i < nrow; i++) m[i][N] = mp[i][N] / mp[i][pivotp[i]]; for (i = 1; i < ranksize; i++) { sr = (i * N) / ranksize; lr = ((i + 1) * N) / ranksize; MPI_Recv (pivotrow, lr - sr, MPI_DOUBLE, i, 1235, MPI_COMM_WORLD, &status); for (j = 0; j < lr - sr; j++) m[sr + j][N] = pivotrow[j]; MPI_Recv (pivotp + sr, lr - sr, MPI_INT, i, 1236, MPI_COMM_WORLD, &status); } /* sort result vector - solution stored in m[i][N] */ for (i = 0; i < N; i++) m[pivotp[i]][0] = m[i][N]; for (i = 0; i < N; i++) m[i][N] = m[i][0]; } else /* process != 0 */ { /* compute and send local part of solution */ for (i = 0; i < nrow; i++) pivotrow[i] = mp[i][N] / mp[i][pivotp[i]]; MPI_Send (pivotrow, nrow, MPI_DOUBLE, 0, 1235, MPI_COMM_WORLD); MPI_Send (pivotp, nrow, MPI_INT, 0, 1236, MPI_COMM_WORLD); } return (TRUE); /* solution found */ } int main (int argc, char **argv) { p_mrow_t a; int solution; int myrank, ranksize; struct timeval t1, t2; struct timezone tz; long sec, usec; MPI_Init (&argc, &argv); MPI_Comm_rank (MPI_COMM_WORLD, &myrank); MPI_Comm_size (MPI_COMM_WORLD, &ranksize); if (myrank == 0) /* process 0 */ { /* process 0 provides matrix to be solved */ a = (p_mrow_t) malloc (N * sizeof(mrow_t)); initialize_system (a); } else /* process != 0 */ a = NULL; /* wait until all processes are running (measuring time) */ MPI_Barrier (MPI_COMM_WORLD); if (myrank == 0) if (gettimeofday(&t1, &tz) == -1) printf("error1\n"); solution = gauss_jordan_elimination (a); MPI_Barrier (MPI_COMM_WORLD); if (myrank == 0) /* process 0 */ { if (gettimeofday(&t2, &tz) == -1) printf("error1\n"); sec = t2.tv_sec - t1.tv_sec; if ((usec = t2.tv_usec - t1.tv_usec) <0) { sec--; usec += MICROSEC; } if (solution == TRUE) print_results (a); else printf ("No solution\n"); printf("%8d sec and %8d u-sec elapsed\n", sec, usec); } MPI_Finalize (); return (0); }