|
ML Version of the Day
|
|
ML Version of the Day
|
#include <math.h> #include "ml_include.h" #ifdef ML_MPI #include "mpi.h" extern int Poisson_getrow(ML_Operator *mat_in, int N_requested_rows, int requested_rows[], int allocated_space, int columns[], double values[], int row_lengths[]); extern int Poisson_matvec(ML_Operator *mat_in, int in_length, double p[], int out_length, double ap[]); extern int Poisson_comm(double x[], void *A_data); extern int send_msg(char *send_buffer, int length, int neighbor); extern int recv_msg(char *recv_buffer, int length, int neighbor, USR_REQ *request); extern int post_msg(char *recv_buffer, int length, int neighbor, USR_REQ *request); /* Simple example corresponding to Poisson on a line with a total of 5 grid points */ /* Thus, the global matrix is */ /* */ /* | 2 -1 | */ /* | -1 2 -1 | */ /* | -1 2 -1 | */ /* | -1 2 -1 | */ /* | -1 2 | */ /* */ /* Processor 0 is assigned global rows 0 and 4 which are stored locally as 1 and */ /* 0. Processor 1 is assigned global rows 1-3 which are stored locally as 0-2 */ /* Processor 0's local matrix is */ /* */ /* | 2 -1 | */ /* | 2 1 -1 | */ /* */ /* and processor 1's local matrix is */ /* */ /* | 2 -1 -1 | */ /* | -1 2 -1 | */ /* | -1 2 -1 | */ /* */ /* See the ML guide for more details. */ int main(int argc, char *argv[]){ ML *ml_object; int i, N_grids = 3, N_levels; double sol[5], rhs[5]; ML_Aggregate *agg_object; int proc, nlocal, nlocal_allcolumns; MPI_Init(&argc,&argv); ML_Set_PrintLevel(15); for (i = 0; i < 5; i++) sol[i] = 0.; for (i = 0; i < 5; i++) rhs[i] = 2.; ML_Create (&ml_object, N_grids); proc = ml_object->comm->ML_mypid; if (ml_object->comm->ML_nprocs != 2) { if (proc == 0) printf("Must be run on two processors\n"); ML_Destroy(&ml_object); MPI_Finalize(); exit(1); } if (proc == 0) {nlocal = 2; nlocal_allcolumns = 4;} else if (proc == 1){nlocal = 3; nlocal_allcolumns = 5;} else {nlocal = 0; nlocal_allcolumns = 0;} ML_Init_Amatrix (ml_object, 0, nlocal, nlocal, &proc); ML_Set_Amatrix_Getrow(ml_object, 0, Poisson_getrow, Poisson_comm, nlocal_allcolumns); ML_Set_Amatrix_Matvec(ml_object, 0, Poisson_matvec); ML_Aggregate_Create(&agg_object); ML_Aggregate_Set_MaxCoarseSize(agg_object,1); N_levels = ML_Gen_MGHierarchy_UsingAggregation(ml_object, 0, ML_INCREASING, agg_object); ML_Gen_Smoother_Jacobi(ml_object, ML_ALL_LEVELS, ML_PRESMOOTHER, 1, ML_DEFAULT); ML_Gen_Solver (ml_object, ML_MGV, 0, N_levels-1); ML_Iterate(ml_object, sol, rhs); if (proc == 0) { printf("sol(0) = %e\n",sol[1]); fflush(stdout); } ML_Comm_GsumInt(ml_object->comm,1); /* just used for synchronization */ if (proc == 1) { printf("sol(1) = %e\n",sol[0]); printf("sol(2) = %e\n",sol[1]); printf("sol(3) = %e\n",sol[2]); fflush(stdout); } ML_Comm_GsumInt(ml_object->comm,1); /* just used for synchronization */ if (proc == 0) { printf("sol(4) = %e\n",sol[0]); fflush(stdout); } ML_Aggregate_Destroy(&agg_object); ML_Destroy(&ml_object); MPI_Finalize(); return 0; } /* Application specific getrow. */ int Poisson_getrow(ML_Operator *mat_in, int N_requested_rows, int requested_rows[], int allocated_space, int cols[], double values[], int row_lengths[]) { int m = 0, i, row, proc, *itemp, start; itemp = (int *) ML_Get_MyGetrowData(mat_in); proc = *itemp; for (i = 0; i < N_requested_rows; i++) { row = requested_rows[i]; if (allocated_space < m+3) return(0); values[m] = 2; values[m+1] = -1; values[m+2] = -1; start = m; if (proc == 0) { if (row == 0) {cols[m++] = 0; cols[m++] = 2; } if (row == 1) {cols[m++] = 1; cols[m++] = 3;} } if (proc == 1) { if (row == 0) {cols[m++] = 0; cols[m++] = 1; cols[m++] = 4;} if (row == 1) {cols[m++] = 1; cols[m++] = 0; cols[m++] = 2;} if (row == 2) {cols[m++] = 2; cols[m++] = 1; cols[m++] = 3;} } row_lengths[i] = m - start; } return(1); } /* Application specific matrix-vector product. */ int Poisson_matvec(ML_Operator *mat_in, int in_length, double p[], int out_length, double ap[]) { int i, proc, *itemp; double new_p[5]; itemp = (int *) ML_Get_MyMatvecData(mat_in); proc = *itemp; for (i = 0; i < in_length; i++) new_p[i] = p[i]; Poisson_comm(new_p, &proc); for (i = 0; i < out_length; i++) ap[i] = 2.*new_p[i]; if (proc == 0) { ap[0] -= new_p[2]; ap[1] -= new_p[3]; } if (proc == 1) { ap[0] -= new_p[1]; ap[0] -= new_p[4]; ap[1] -= new_p[2]; ap[1] -= new_p[0]; ap[2] -= new_p[3]; ap[2] -= new_p[1]; } return 0; } /* Communication routine that should be performed before doing a matvec */ /* See ML Guide. */ int Poisson_comm(double x[], void *A_data) { int proc, neighbor, length, *itemp; double send_buffer[2], recv_buffer[2]; MPI_Request request; itemp = (int *) A_data; proc = *itemp; length = 2; if (proc == 0) { neighbor = 1; send_buffer[0] = x[0]; send_buffer[1] = x[1]; post_msg((char *) recv_buffer, length, neighbor, &request); send_msg((char *) send_buffer, length, neighbor); recv_msg((char *) recv_buffer, length, neighbor, &request); x[2] = recv_buffer[1]; x[3] = recv_buffer[0]; } else { neighbor = 0; send_buffer[0] = x[0]; send_buffer[1] = x[2]; post_msg((char *) recv_buffer, length, neighbor, &request); send_msg((char *) send_buffer, length, neighbor); recv_msg((char *) recv_buffer, length, neighbor, &request); x[3] = recv_buffer[1]; x[4] = recv_buffer[0]; } return 0; } /* Simple communication wrappers for use with MPI */ int send_msg(char *send_buffer, int length, int neighbor) { ML_Comm_Send(send_buffer, length*sizeof(double), neighbor, 123, MPI_COMM_WORLD); return 0; } int recv_msg(char *recv_buffer, int length, int neighbor, USR_REQ *request) { MPI_Status status; MPI_Wait(request, &status); return 0; } int post_msg(char *recv_buffer, int length, int neighbor, USR_REQ *request) { int type = 123; ML_Comm_Irecv(recv_buffer, length*sizeof(double), &neighbor, &type, MPI_COMM_WORLD, request); return 0; } #else int main(int argc, char *argv[]) { printf("In order to use this example, you must configure with the option\n--with-mpi-compilers=full_path_to_your_mpi_compilers and recompile.\n"); /* returns ok not to break the test harness */ return(EXIT_SUCCESS); } #endif /* ifdef ML_MPI */
1.7.4
