A small (more or less) example of programming with MPI in C. This is a copy of my solution to the second practical for Parallel Programming Languages and Systems, which our course organiser Dr Murray Cole has kindly allowed me to post in this blog.
The algorithm calculates adaptive quadrature using MPI and the "bag of tasks" approach.
"Adaptive Quadrature is a recursive algorithm that computes an approximation of the integral of a function F(x), using static quadrature rules on adaptively reļ¬ned sub-intervals of the integration domain."
To compile and run it you need some kind of MPI libraries installed on your system (I used OpenMPI) and the following commands:
mpicc -o
to compile and
mpirun -c 5
to run. Here's the source code (a bit lengthy, need to find a way to minimise code listings...):
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#define EPSILON 1e-3
#define F(arg) cosh(arg)*cosh(arg)*cosh(arg)*cosh(arg)
#define A 0.0
#define B 5.0
#define SLEEPTIME 1
/***************************************************************************
Notes on implementation:
Tested on my own machine with 2 CPU cores. MPI interface: OpenMPI 1.3.2-3ubuntu1.1
Description:
The implementation is based on a standard "bag of tasks" technique. The farmer
(aka controller) and worker functions initiate two loops:
Farmer:
After initial declarations, farmer loop is initiated with exit conditions of:
- All workers have finished computing area (maintained by the int idle variable
which increments every time a worker returns results and decrements each time
some worker is passed task to perform)
- Stack is empety (additional method isEmpty was implemented in stack.c to
check that condition)
The loop starts with wild-card synchronous MPI receive function. When data is
received, controller checks whether it is a new task or computed area (indicated by a tag).
If it is a new task, it is pushed into the stack, otherwise, adds the received value to
the total area. Then idle process counter is incremented and the value of the proc_waiting
array slot, corresponding to the worker's ID, changes to 1, indicating that this
worker is idle. Then, if stack is not empty, the controller iterates over proc_waiting,
starting from the point where it last finished (which ensures that farmer would
not be hijacked by any worker), and if the current value in the array is 1, sends new task
to the worker and continues the main loop. When stack is empty and all workers are idle it means
that the area has been computed, farmer breaks the main loop and sends exit signal to all workers.
Worker:
Worker loop starts with synchronous MPI receive comand waiting for input from
controller. The received input is processed according to the algorithm provided
and the results are sent back to the controller. When exit signal (message with
specific tag) is received, worker breaks the loop and terminates.
MPI primitives:
Synchronous blocking wild-card MPI receive is used on the controller in order
to avoid useless iterating over all workers with asynchronous receive scanning
for results. Blocking receive ensures synchronisation with all workers. MPI gather
is not used since workers can finish their tasks with different speed, and waiting
for input from all workers together on every iteration of the controller loop would
have decreased the overall performance, leaving alone the fact that implementation
would have not been so straight forward in case of using MPI gather.
Semi-synchronous (but blocking) send is used since neither the farmer nor the workers need to wait
till data is received by other side as long as they know it is going to be received eventually,
which is ensured by the blocking property of the mechanism.
****************************************************************************/
/*
Results:
mpirun -c 10 ./aquadPartA
Area=7583461.801486
Tasks Per Process
0 1 2 3 4 5 6 7 8 9
0 733 728 709 744 746 748 730 722 707
*/
typedef struct stack_node_tag stack_node;
typedef struct stack_tag stack;
struct stack_node_tag {
double data[2];
stack_node *next;
};
struct stack_tag {
stack_node *top;
};
stack *new_stack();
void free_stack(stack *);
void push(double *, stack *);
double *pop (stack *);
int isEmpty(stack*);
int *tasks_per_process;
double farmer(int);
void worker(int);
int main(int argc, char **argv ) {
int i, myid, numprocs;
double area, a, b;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
MPI_Comm_rank(MPI_COMM_WORLD,&myid);
if(numprocs < 2) {
fprintf(stderr, "ERROR: Must have at least 2 processes to run\n");
MPI_Finalize();
exit(1);
}
if (myid == 0) { // Farmer
// init counters
tasks_per_process = (int *) malloc(sizeof(int)*(numprocs));
for (i=0; i<numprocs; i++) {
tasks_per_process[i]=0;
}
}
if (myid == 0) { // Farmer
area = farmer(numprocs);
} else { //Workers
worker(myid);
}
if(myid == 0) {
fprintf(stdout, "Area=%lf\n", area);
fprintf(stdout, "\nTasks Per Process\n");
for (i=0; i<numprocs; i++) {
fprintf(stdout, "%d\t", i);
}
fprintf(stdout, "\n");
for (i=0; i<numprocs; i++) {
fprintf(stdout, "%d\t", tasks_per_process[i]);
}
fprintf(stdout, "\n");
free(tasks_per_process);
}
MPI_Finalize();
return 0;
}
double farmer(int numprocs) {
int n_workers = numprocs - 1;
//total number of idle workers
int idle=0;
//iterator over the list of workers
int iter=0;
// list of workers. values: 1 if waiting for input, 0 if computing
int* proc_waiting = (int*) malloc(sizeof(int)*(n_workers));
double result = 0;
MPI_Status status;
double* temp = (double*) malloc(sizeof(double)*2);
stack* bag;
bag = new_stack();
temp[0] = A;
temp[1] = B;
push(temp,bag);
int i=0;
for (i;i<n_workers;i++){
proc_waiting[i] = 0;
}
//Controller loop
do{
// Receiving data from workers.
MPI_Recv(temp, 2, MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
idle++;
proc_waiting[status.MPI_SOURCE - 1] = 1;
if (status.MPI_TAG == 1){
result += temp[0];
}else{
push(temp,bag);
int tag = status.MPI_TAG;
int source = status.MPI_SOURCE;
MPI_Recv(temp, 2, MPI_DOUBLE, source, tag, MPI_COMM_WORLD, &status);
push(temp,bag);
}
//Iterating over proc_waiting list to find and idle worker and send a task to it
while(!isEmpty(bag) && idle>0){
if (proc_waiting[iter]){
MPI_Send(pop(bag),2,MPI_DOUBLE,iter+1,0,MPI_COMM_WORLD);
proc_waiting[iter]=0;
idle--;
tasks_per_process[iter+1]++;
}
iter = (iter+1)%n_workers;
}
}while(!(isEmpty(bag) && idle==n_workers));
i = 0;
// Sending exit signal to all workers
for (i;i<n_workers;i++){
temp[0] = 0;
temp[1] = 0;
MPI_Send(temp,2,MPI_DOUBLE,i+1,1,MPI_COMM_WORLD);
}
return result;
}
// mypid argument is not used in this implementation
void worker(int mypid) {
MPI_Status status;
double buf[]={0,0};
// Sending 0 to the farmer to indicate that the worker has initialized and
// ready to receive a task.
MPI_Send(buf,2,MPI_DOUBLE,0,1,MPI_COMM_WORLD);
// Worker loop
while (1){
MPI_Recv(buf,2,MPI_DOUBLE,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
int tag = status.MPI_TAG;
// exit signal
if (tag == 1)
break;
double left = buf[0];
double right = buf[1];
double lrarea = (F(left) + F(right)) * (right - left) / 2;
double fleft = F(left);
double fright = F(right);
double mid, fmid, larea, rarea;
mid = (left + right) / 2;
fmid = F(mid);
larea = (fleft + fmid) * (mid - left) / 2;
rarea = (fmid + fright) * (right - mid) / 2;
if( fabs((larea + rarea) - lrarea) > EPSILON ) {
buf[0] = left;
buf[1] = mid;
// Sending the resulted task 1
MPI_Send(buf,2,MPI_DOUBLE,0,0,MPI_COMM_WORLD);
buf[0] = mid;
buf[1] = right;
// Sending the resulted task 2
MPI_Send(buf,2,MPI_DOUBLE,0,0,MPI_COMM_WORLD);
}else{
buf[0] = larea+rarea;
buf[1] = 0;
MPI_Send(buf,2,MPI_DOUBLE,0,1,MPI_COMM_WORLD);
}
}
}
int isEmpty(stack * s){
if (s==NULL || s->top == NULL){
return 1;
}else{
return 0;
}
}
// Simple type for stack of doubles
// creating a new stack
stack * new_stack()
{
stack *n;
n = (stack *) malloc (sizeof(stack));
n->top = NULL;
return n;
}
// cleaning up after use
void free_stack(stack *s)
{
free(s);
}
// Push data to stack s, data has to be an array of 2 doubles
void push (double *data, stack *s)
{
stack_node *n;
n = (stack_node *) malloc (sizeof(stack_node));
n->data[0] = data[0];
n->data[1] = data[1];
if (s->top == NULL) {
n->next = NULL;
s->top = n;
} else {
n->next = s->top;
s->top = n;
}
}
// Pop data from stack s
double * pop (stack * s)
{
stack_node * n;
double *data;
if (s == NULL || s->top == NULL) {
return NULL;
}
n = s->top;
s->top = s->top->next;
data = (double *) malloc(2*(sizeof(double)));
data[0] = n->data[0];
data[1] = n->data[1];
free (n);
return data;
}
No comments:
Post a Comment