Wave Equation Implementation (Version 3)
Wave Equation Implementation (Version 3)
Take a look at the code below and note the MPI functions.
Click on the name of each MPI routine to read a
short description of that routine's purpose and syntax.
/*---------------------------------------------------------------------
*
* This program implements the concurrent wave equation described
* in Chapter 5 of Fox et al., 1988, Solving Problems on Concurrent
* Processors, vol 1.
*
* A vibrating string is decomposed into points. Each process is
* responsible for updating the amplitude of a number of points over
* time.
*
* At each iteration, each process exchanges boundary points with
* nearest neighbors. This version uses low level sends and receives
* to exchange boundary points.
*
* This version illustrates use of persistent communication.
*
* Revised in 11/97, Saleh Elmohamed.
*---------------------------------------------------------------------*/
#include
#include
#include
#include
#include
#include "mpi.h"
#define MASTER 0
#define MAXPROC 8
#define MAXPOINTS 1000
#define MAXSTEPS 1000000
int E_RtoL = 10;
int E_LtoR = 20;
int E_OUT1 = 30;
int E_OUT2 = 40;
void get_data(void);
void init_line(void);
void pcomm_create (int right, int left);
void update (int left, int right);
void output_master(void);
void output_workers(void);
MPI_Request request;
MPI_Status status;
MPI_Request Prequest[4];
MPI_Status Pstatus[4];
int rank, /* rank */
ntask, /* number of tasks */
nsteps, /* number of time steps */
tpoints, /* total points along string */
npoints, /* number of points handled by this task */
first, /* index of 1st point handled by this task */
rcode; /* generic return code */
double values[MAXPOINTS+2],/* values at time t */
oldval[MAXPOINTS+2], /* values at time (t-dt) */
newval[MAXPOINTS+2]; /* values at time (t+dt) */
/*------------------------------------------------------------------------
* Obtain input values from user
*------------------------------------------------------------------------ */
void get_data(void)
{
char tchar[8];
int buffer[2];
if (rank == MASTER) {
/* get number of points, number of iterations from user */
tpoints = 1000;
nsteps = 10000;
#if 0
while ((tpoints < ntask) || (tpoints > MAXPOINTS))
{
printf("Enter number of points along vibrating string\n");
scanf("%s", tchar);
tpoints = atoi(tchar);
if ((tpoints < ntask) || (tpoints > MAXPOINTS))
printf("enter value between %d and %d\n", ntask, MAXPOINTS);
}
while ((nsteps < 1) || (nsteps > MAXSTEPS))
{
printf("Enter number of time steps\n");
scanf("%s", tchar);
nsteps = atoi(tchar);
if ((nsteps < 1) || (nsteps > MAXSTEPS))
printf("enter value between 1 and %d\n", MAXSTEPS);
}
#endif
printf("%d: points = %d, steps = %d\n", rank, tpoints, nsteps);
buffer[0] = tpoints;
buffer[1] = nsteps;
}
/* then broadcast number of points and number of iterations */
MPI_Bcast(buffer, 2,
MPI_INT, MASTER, MPI_COMM_WORLD);
tpoints = buffer[0];
nsteps = buffer[1];
}
/*------------------------------------------------------------------------
* Initialize points on line
*--------------------------------------------------------------------- */
void init_line(void)
{
int nmin, nleft, npts, i, j, k;
double x, fac;
double PI;
/* calculate initial values based on sine curve */
nmin = tpoints/ntask;
nleft = tpoints%ntask;
PI = 4.0 * atan(1.0);
fac = 2.0 * PI;
for (i = 0, k = 0; i < ntask; i++)
{
npts = (i < nleft) ? nmin + 1 : nmin;
if (rank == i)
{
first = k + 1;
npoints = npts;
printf ("%d: first = %d, npoints = %d\n", rank, first, npts);
for (j = 1; j <= npts; j++, k++)
{
x = (double)k/(double)(tpoints - 1);
values[j] = sin (fac * x);
}
}
else k += npts;
}
for (i = 1; i <= npoints; i++)
oldval[i] = values[i];
}
/*-------------------------------------------------------------------------
* Create persistent communication requests
*-------------------------------------------------------------------------*/
void pcomm_create (int right, int left)
{
int i;
/* Initialize all requests to NULL
* All tasks wait for all four requests, even though some tasks
* do not use all four requests */
for (i=0;i<4;i++)
Prequest[i] = MPI_REQUEST_NULL;
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* Step 1: Create persistent communication requests for: *
* sending to left-hand neighbor Prequest[0] *
* receiving from left-hand neighbor Prequest[1] *
* sending to right-hand neighbor Prequest[2] *
* receiving from right-hand neighbor Prequest[3] *
* Remember: the arguments to the initialization calls are the *
* same as for the non-blocking communications you are replacing *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/* create requests for communication with left-hand neighbor */
if (left != MPI_PROC_NULL)
{
MPI_Send_init(&values[1], 1,
MPI_DOUBLE, left, E_RtoL, MPI_COMM_WORLD, &Prequest[0]);
MPI_Recv_init(&values[0], 1,
MPI_DOUBLE, left, E_LtoR, MPI_COMM_WORLD, &Prequest[1]);
}
/* create requests for communication with right-hand neighbor */
if (right != MPI_PROC_NULL)
{
MPI_Send_init(&values[npoints], 1,
MPI_DOUBLE, right,
E_LtoR, MPI_COMM_WORLD, &Prequest[2]);
MPI_Recv_init(&values[npoints+1], 1,
MPI_DOUBLE, right, E_RtoL, MPI_COMM_WORLD, &Prequest[3]);
}
}
/*---------------------------------------------------------------------
* Calculate new values using wave equation
*---------------------------------------------------------------------*/
void do_math(int i)
{
const double dtime = 0.3;
const double c = 1.0;
const double dx = 1.0;
double tau, sqtau;
tau = (c * dtime / dx);
sqtau = tau * tau;
newval[i] = (2.0 * values[i]) - oldval[i]
+ (sqtau * (values[i-1] - (2.0 * values[i]) + values[i+1]));
}
/*---------------------------------------------------------------------
* All tasks update their points a specified number of times
*---------------------------------------------------------------------*/
void update(int left, int right)
{
int i, j;
struct timeval tv1, tv2;
int dt=0;
/* update values for each point along string */
for (i = 1; i <= nsteps; i++)
{
/* * * * * * * * * * * * * * * * * * * * * *
* Step 2: *
* Replace the four non-blocking communication calls with *
* calls that "start" the four persistent communications. *
* * * * * * * * * * * * * * * * * * * * * */
/* time the communication */
gettimeofday(&tv1, (struct timeval*)0);
/* exchange data with "left-hand" neighbor */
if (left != MPI_PROC_NULL)
{
MPI_Start(&Prequest[0]);
MPI_Start(&Prequest[1]);
}
/* exchange data with "right-hand" neighbor */
if (right != MPI_PROC_NULL)
{
MPI_Start(&Prequest[2]);
MPI_Start(&Prequest[3]);
}
MPI_Waitall(4, Prequest, Pstatus);
/* complete the timing */
gettimeofday(&tv2, (struct timeval*)0);
dt = dt + ((tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec
- tv1.tv_usec);
/* update points along line */
for (j = 1; j <= npoints; j++)
{
/* global endpoints */
if ((first + j - 1 == 1) || (first + j - 1 == tpoints))
newval[j] = 0.0;
else
do_math(j);
}
for (j = 1; j <= npoints; j++)
{
oldval[j] = values[j];
values[j] = newval[j];
}
}
printf("time for communication in update loop uSec = %8d\n", dt);
}
/*------------------------------------------------------------------
* Master receives results from workers and prints
*------------------------------------------------------------------ */
void output_master(void)
{
int i, start, npts, buffer[2], istep, source;
double results[MAXPOINTS];
/* store worker's results in results array */
for (i = 1; i < ntask; i++)
{
/* receive number of points and first point */
MPI_Recv(buffer, 2, MPI_INT,
MPI_ANY_SOURCE, E_OUT1, MPI_COMM_WORLD, &status);
start = buffer[0];
npts = buffer[1];
/* receive results */
source = status.MPI_SOURCE;
MPI_Recv(&results[start-1],
npts, MPI_DOUBLE, source, E_OUT2, MPI_COMM_WORLD, &status);
}
/* store master's results in results array */
for (i = first; i < first + npoints; i++)
results[i-1] = values[i];
istep = (tpoints <= 10) ? 1: tpoints/10;
printf ("Points for validation:\n");
for (i = 0; i < tpoints; i+=istep)
printf ("%d:%4.2f ", i, results[i]);
if (i-istep != tpoints - 1)
printf ("%d:%4.2f ", tpoints-1, results[tpoints-1]);
printf("\n");
}
/*----------------------------------------------------------------
* Workers send the updated values to the master
*----------------------------------------------------------------*/
void output_workers(void)
{
int buffer[2];
/* send first point and number of points handled to master */
buffer[0] = first;
buffer[1] = npoints;
MPI_Isend(buffer, 2, MPI_INT,
0, E_OUT1, MPI_COMM_WORLD, &request);
MPI_Wait(&request, &status);
/* send results to master */
MPI_Isend(&values[1], npoints,
MPI_DOUBLE, 0, E_OUT2, MPI_COMM_WORLD, &request);
MPI_Wait(&request, &status);
}
/*------------------------------------------------------------------
* Main program
*------------------------------------------------------------------*/
int main(int argc, char **argv)
{
int left, right, i;
/* learn number of tasks and rank in MPI_COMM_WORLD */
rcode = MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD,
&rank);
MPI_Comm_size(MPI_COMM_WORLD,
&ntask);
if (rcode >= 0) printf ("%d: Wave Program running\n", rank);
/* determine left and right neighbors */
if (rank == ntask-1) right = MPI_PROC_NULL;
else right = rank + 1;
if (rank == 0) left = MPI_PROC_NULL;
else left = rank - 1;
/* get program parameters and initialize wave values */
get_data();
init_line();
/* create persistent communication requests */
pcomm_create(right, left);
/* update values along the line for nstep time steps */
update(left, right);
/* * * * * * * * * * * * * * * * * * * * *
* Step 3: *
*--- deallocate persistent communication requests here */
for (i=0; i<4; i++)
if (Prequest[i] != MPI_REQUEST_NULL) MPI_Request_free (&Prequest[i]);
/*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*/
/* collect results and print */
if (rank == MASTER) output_master();
else output_workers();
MPI_Finalize();
return(0);
}
A 4-processor run:
1: Wave Program running
0: Wave Program running
2: Wave Program running
3: Wave Program running
0: points = 1000, steps = 10000
0: first = 1, npoints = 250
1: first = 251, npoints = 250
2: first = 501, npoints = 250
3: first = 751, npoints = 250
time for communication in update loop uSec = 27450719
time for communication in update loop uSec = 28843668
time for communication in update loop uSec = 23671080
time for communication in update loop uSec = 23318864
Points for validation:
0:0.00 100:0.59 200:0.95 300:0.95 400:0.59 500:-0.00 600:-0.59 700:-0.95 800:-0.95 900:-0.58 999:0.00