/* parallel_jacobi.c -- parallel implementation of Jacobi's method * for solving the linear system Ax = b. Uses block distribution * of vectors and block-row distribution of A. * * Input: * n: order of system * tol: convergence tolerance * max_iter: maximum number of iterations * A: coefficient matrix * b: right-hand side of system * * Output: * x: the solution if the method converges * max_iter: if the method fails to converge * * Notes: * 1. A should be strongly diagonally dominant in * order to insure convergence. * 2. A, x, and b are statically allocated. * * See Chap 10, pp. 220 & ff in PPMPI. */ #include #include "mpi.h" #include #define Swap(x,y) {float* temp; temp = x; x = y; y = temp;} #define MAX_DIM 12 typedef float MATRIX_T[MAX_DIM][MAX_DIM]; int Parallel_jacobi( MATRIX_T A_local /* in */, float x_local[] /* out */, float b_local[] /* in */, int n /* in */, float tol /* in */, int max_iter /* in */, int p /* in */, int my_rank /* in */); void Read_matrix(char* prompt, MATRIX_T A_local, int n, int my_rank, int p); void Read_vector(char* prompt, float x_local[], int n, int my_rank, int p); void Print_matrix(char* title, MATRIX_T A_local, int n, int my_rank, int p); void Print_vector(char* title, float x_local[], int n, int my_rank, int p); main(int argc, char* argv[]) { int p; int my_rank; MATRIX_T A_local; float x_local[MAX_DIM]; float b_local[MAX_DIM]; int n; float tol; int max_iter; int converged; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &p); MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); if (my_rank == 0) { printf("Enter n, tolerance, and max number of iterations\n"); scanf("%d %f %d", &n, &tol, &max_iter); } MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&tol, 1, MPI_FLOAT, 0, MPI_COMM_WORLD); MPI_Bcast(&max_iter, 1, MPI_INT, 0, MPI_COMM_WORLD); Read_matrix("Enter the matrix", A_local, n, my_rank, p); Read_vector("Enter the right-hand side", b_local, n, my_rank, p); converged = Parallel_jacobi(A_local, x_local, b_local, n, tol, max_iter, p, my_rank); if (converged) Print_vector("The solution is", x_local, n, my_rank, p); else if (my_rank == 0) printf("Failed to converge in %d iterations\n", max_iter); MPI_Finalize(); } /* main */ /*********************************************************************/ /* Return 1 if iteration converged, 0 otherwise */ /* MATRIX_T is a 2-dimensional array */ int Parallel_jacobi( MATRIX_T A_local /* in */, float x_local[] /* out */, float b_local[] /* in */, int n /* in */, float tol /* in */, int max_iter /* in */, int p /* in */, int my_rank /* in */) { int i_local, i_global, j; int n_bar; int iter_num; float x_temp1[MAX_DIM]; float x_temp2[MAX_DIM]; float* x_old; float* x_new; float Distance(float x[], float y[], int n); n_bar = n/p; /* Initialize x */ MPI_Allgather(b_local, n_bar, MPI_FLOAT, x_temp1, n_bar, MPI_FLOAT, MPI_COMM_WORLD); x_new = x_temp1; x_old = x_temp2; iter_num = 0; do { iter_num++; /* Interchange x_old and x_new */ Swap(x_old, x_new); for (i_local = 0; i_local < n_bar; i_local++){ i_global = i_local + my_rank*n_bar; x_local[i_local] = b_local[i_local]; for (j = 0; j < i_global; j++) x_local[i_local] = x_local[i_local] - A_local[i_local][j]*x_old[j]; for (j = i_global+1; j < n; j++) x_local[i_local] = x_local[i_local] - A_local[i_local][j]*x_old[j]; x_local[i_local] = x_local[i_local]/ A_local[i_local][i_global]; } MPI_Allgather(x_local, n_bar, MPI_FLOAT, x_new, n_bar, MPI_FLOAT, MPI_COMM_WORLD); } while ((iter_num < max_iter) && (Distance(x_new,x_old,n) >= tol)); if (Distance(x_new,x_old,n) < tol) return 1; else return 0; } /* Jacobi */ /*********************************************************************/ float Distance(float x[], float y[], int n) { int i; float sum = 0.0; for (i = 0; i < n; i++) { sum = sum + (x[i] - y[i])*(x[i] - y[i]); } return sqrt(sum); } /* Distance */ /*********************************************************************/ void Read_matrix( char* prompt /* in */, MATRIX_T A_local /* out */, int n /* in */, int my_rank /* in */, int p /* in */) { int i, j; MATRIX_T temp; int n_bar; n_bar = n/p; /* Fill dummy entries in temp with zeroes */ for (i = 0; i < n; i++) for (j = n; j < MAX_DIM; j++) temp[i][j] = 0.0; if (my_rank == 0) { printf("%s\n", prompt); for (i = 0; i < n; i++) for (j = 0; j < n; j++) scanf("%f",&temp[i][j]); } MPI_Scatter(temp, n_bar*MAX_DIM, MPI_FLOAT, A_local, n_bar*MAX_DIM, MPI_FLOAT, 0, MPI_COMM_WORLD); } /* Read_matrix */ /*********************************************************************/ void Read_vector( char* prompt /* in */, float x_local[] /* out */, int n /* in */, int my_rank /* in */, int p /* in */) { int i; float temp[MAX_DIM]; int n_bar; n_bar = n/p; if (my_rank == 0) { printf("%s\n", prompt); for (i = 0; i < n; i++) scanf("%f", &temp[i]); } MPI_Scatter(temp, n_bar, MPI_FLOAT, x_local, n_bar, MPI_FLOAT, 0, MPI_COMM_WORLD); } /* Read_vector */ /*********************************************************************/ void Print_matrix(char* title, MATRIX_T A_local, int n, int my_rank, int p); void Print_matrix( char* title /* in */, MATRIX_T A_local /* in */, int n /* in */, int my_rank /* in */, int p /* in */) { int i, j; MATRIX_T temp; int n_bar; n_bar = n/p; MPI_Gather(A_local, n_bar*MAX_DIM, MPI_FLOAT, temp, n_bar*MAX_DIM, MPI_FLOAT, 0, MPI_COMM_WORLD); if (my_rank == 0) { printf("%s\n", title); for (i = 0; i < n; i++) { for (j = 0; j < n; j++) printf("%4.1f ", temp[i][j]); printf("\n"); } } } /* Print_matrix */ /*********************************************************************/ void Print_vector(char* title, float x_local[], int n, int my_rank, int p); void Print_vector( char* title /* in */, float x_local[] /* in */, int n /* in */, int my_rank /* in */, int p /* in */) { int i; float temp[MAX_DIM]; int n_bar; n_bar = n/p; MPI_Gather(x_local, n_bar, MPI_FLOAT, temp, n_bar, MPI_FLOAT, 0, MPI_COMM_WORLD); if (my_rank == 0) { printf("%s\n", title); for (i = 0; i < n; i++) printf("%4.1f ", temp[i]); printf("\n"); } } /* Print_vector */