6
\$\begingroup\$

I am learning OpenMP+MPI hybrid programming. As an example I have chosen Gauss-Seidel+SOR. My implementation uses MPI_THREAD_FUNNELED style hybrid programming, overlapping of communication and computation, neighbourhood collectives, MPI IO etc. Based on my analysis and running many executions it seems to be correct (cannot guarantee however, since I am relatively new to this area).

I am looking for general comments and feedback for improving my OpenMP, MPI and hybrid programming skills. Additionally, I am particularly interested in learning the root cause of the following scalability issue:

Unfortunately MPI outperforms the OpenMP scalability on a single compute node:

enter image description here

  • Hardware: 2x Xeon E5-2698 v3, 16 cores each
  • Compiler: Cray C : Version 8.3.4
  • Compiler Flags: -h pragma=omp -h ipa5
  • MPI: cray-mpich/7.2.2
  • Input functions: f(x) = sqrt(x) and h(x) = 10 * sin(10 * x)

My hypothesis is that this could be due to the ccNUMA architecture of the Haswell-EP CPU. Wikipedia says:

ccNUMA may perform poorly when multiple processors attempt to access the same memory area in rapid succession.

This is the case in my code. Additionally: If I have 8 OpenMP threads, then they are bound by aprun to cores which are all connected by the same ring bus. When I increase to 16 cores the overhead for maintaining cache coherence increases and threads are bound to the cores of the first CPU and communication therefore has to go through the buffered switches. In the case of 32 threads even through QPI. The MPI version will require a lot less cache synchronisation, since the data is spread among several processes. Note that in the OpenMP case the complete data fits into one memory page.

#include "functions.h"
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/param.h>
#include <float.h>
#include <stdbool.h>
#include <omp.h>

void iteration_loop(double* u, const int right_most, const int right_ghost, const int left_most, const int left_ghost);
void redblack(double* u, int right_most, int right_ghost, int left_most, int left_ghost);
void cacheFunctions(const int start, const int end, const int N_local);
void writeresult(const int N_local, const double* u);
double step(double* u, const int n);
extern double f(double x);
extern double r(double x);

static double * ff;
static double * rr;

#define N 6000 /* number of discretization points */
static const double h = 1. / (N - 1); /* discretization step */
static const int MAX_ITER = 1000000; /* iteration limit */
static const int TERMINATION_DETECTION_INTERVALL = 8000;
static const double EPSILON = 0.0000001; /* convergence criterion */
static const double THETA = 0.9; /* relaxation factor for SOR */

static int world_rank, world_size; /* MPI */
#define ndims 1 /* one dimensional Cartesian coordinate system */
static const int peridos[ndims] = {0}; /* grid is not periodic  */
static int dims[ndims] = {0}; /* division of processors in a Cartesian
                                 0 allows MPI_Dims_create to override element */
static MPI_Comm COMM_GRID; /* Cartesian communicator */
static int my_coords[ndims]; /* Cartesian of process */

/* MPI_Request objects for non-blocking communication */
static MPI_Request neighbor = MPI_REQUEST_NULL;
static MPI_Request reduce = MPI_REQUEST_NULL;

static double start, end; /* used to calculate runtime */


////////////////////////////////////////////////////////////////////////////////

int main(int argc, char** argv) {

    /* This code requires MPI_THREAD_FUNNELED. Quit if not met. */
    int provided_thread_support;
    MPI_Init_thread(&argc, &argv, MPI_THREAD_FUNNELED, &provided_thread_support);
    if (provided_thread_support == MPI_THREAD_SINGLE) {
        fprintf(stderr, "Provided only MPI_THREAD_SINGLE\n");
        exit(EXIT_FAILURE);
    }

    MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &world_size);

    /* Divides processors in a Cartesian grid. */
    MPI_Dims_create(world_size, ndims, dims);
    /* Create communicator with Cartesian topology */
    MPI_Cart_create(MPI_COMM_WORLD, ndims, dims, peridos, 0, &COMM_GRID);
    MPI_Cart_coords(COMM_GRID, world_rank, ndims, my_coords);

    /* data distribution, number of cells assigned to processor */
    int N_local = ((N - 2 + dims[0] - my_coords[0] - 1) / world_size);
    /* two additional memory cells are required for ghost cells / boundary condition*/
    int MEM_local = N_local + 2;

    /* Indices */
    int left_ghost = 0; /* Index of left ghost cell */
    int left_most = 1; /* Index of left most local cell */
    int right_ghost = N_local + 1; /* Index of right ghost cell */
    int right_most = N_local; /* Index of right most local cell */

    /* Time measurement */
    MPI_Barrier(MPI_COMM_WORLD);
    start = MPI_Wtime();

    /* Cached r and f */
    rr = malloc(N_local * sizeof (double));
    ff = malloc(N_local * sizeof (double));
    cacheFunctions(1, right_most, N_local);

    /* Allocate initial guess (all zero) */
    double* u = calloc(MEM_local, sizeof (double));
    if (u == NULL) {
        perror("calloc");
        exit(EXIT_FAILURE);
    }

    iteration_loop(u, right_most, right_ghost, left_most, left_ghost);

    /* Time measurement */
    MPI_Barrier(MPI_COMM_WORLD);
    end = MPI_Wtime();

    writeresult(N_local, u);

    if (world_rank == 0) {
        printf("\n\nExecution time: %f\n", end - start);
    }

    /* Cleanup */
    free(u);
    free(rr);
    free(ff);

    MPI_Finalize();
    return (EXIT_SUCCESS);
}


////////////////////////////////////////////////////////////////////////////////

/* Caches results of f(x) and r(x) in ff[] and rr[] respectively */
void cacheFunctions(
        const int start,
        const int end,
        const int N_local) {
    for (int n = start; n <= end; ++n) {
        const double x_n = (my_coords[0] * N_local + n) * h;
        ff[n] = f(x_n);
        rr[n] = r(x_n);
    }
}

////////////////////////////////////////////////////////////////////////////////

/* Perform SOR until termination condition is met.  */
void iteration_loop(
        double* u,
        const int right_most,
        const int right_ghost,
        const int left_most,
        const int left_ghost) {

    /* Maximum deviation from current and previous iteration */
    double maxdelta = 0.0;
    double oldmaxdelta = DBL_MAX;

    /* For iteration loop control */
    bool terminate = false;
    int iter;

    /* Iteration loop */
#    pragma omp parallel default(shared) private(iter) firstprivate(terminate)
    for (iter = 0; iter < MAX_ITER; ++iter) {
        /* Breaking Out of iteration loop  */
        if (!terminate) {
            /* Communicate ghost cells */
#            pragma omp master
            {
                redblack(u, right_most, right_ghost, left_most, left_ghost);
            }
            /* Update cells 2 to right_most - 1 which are independent
               of ghost cells  */
#            pragma omp for schedule(static) reduction(max:maxdelta)
            for (int n = 2; n <= right_most - 1; ++n) {
                double delta = step(u, n);
                maxdelta = MAX(maxdelta, delta);
            }

#            pragma omp master
            {
                /* Wait until ghost cells are completely communicated */
                MPI_Wait(&neighbor, MPI_STATUSES_IGNORE);
                /* Update cells dependent on ghost cells */
                double delta;
                delta = step(u, left_most);
                maxdelta = MAX(maxdelta, delta); /* only visible to master, acceptable */
                delta = step(u, right_most);
                maxdelta = MAX(maxdelta, delta); /* dito */
            }

            /* If reduction of termination condition was initiated in previous
             * iteration. */
            if ((iter - 1) % TERMINATION_DETECTION_INTERVALL == 0) {
                /* Wait until reduction of global maximum deviation has completed */
#                pragma omp master
                {
                    MPI_Wait(&reduce, MPI_STATUSES_IGNORE);
                }
                /* All threads need to wait for master to complete MPI reduction
                 * and they need to see the result in oldmaxdelta in order to
                 * know if they should break the loop */
#                pragma omp barrier
                /* If maxdelta of previous iteration is smaller than EPSILON, terminate */
                if (oldmaxdelta < EPSILON) {
                    if (world_rank == 0) {
                        printf("[%d] breaking after %d of %d iterations\n", omp_get_thread_num(), iter, MAX_ITER);
                    }
                    terminate = true;
                }
            }
#            pragma omp master
            {
                /* Decides if termination detection should take place in this iteration */
                if (iter % TERMINATION_DETECTION_INTERVALL == 0) {
                    //printf("2) iter=%d, oldmaxdelta=%.9f, maxdelta=%.9f \n", iter, oldmaxdelta, maxdelta);
                    /* Save maxdelta to oldmaxdelta in order to overlap reduction with computation */
                    oldmaxdelta = maxdelta;
                    /* Perform reduction to find global maximum deviation in current iteration */
                    MPI_Iallreduce(MPI_IN_PLACE, &oldmaxdelta, 1, MPI_DOUBLE,
                            MPI_MAX, COMM_GRID, &reduce);
                }
            }
            /* master resetts maxdelta, after initiating Iallreduce */
#            pragma omp master
            {
                maxdelta = 0;
            }
            /* all threads need to wait until maxdelta has been updated */
#            pragma omp barrier
        }
    }
}

////////////////////////////////////////////////////////////////////////////////

/* Perform SOR on cell  */
double step(double* u, const int n) {
    const double last_u = u[n];
    u[n] = (u[n - 1] + u[n + 1] - h * h * ff[n]) / (2. - h * h * rr[n]);
    u[n] = u[n] + THETA * (u[n] - last_u);
    return fabs(last_u - u[n]);
}

////////////////////////////////////////////////////////////////////////////////

/* Perform redback communication. */
void redblack(
        double* u,
        const int right_most,
        const int right_ghost,
        const int left_most,
        const int left_ghost) {


    int counts[] = {1, 1};
    int sdispls[] = {left_most, right_most};
    int rdispls[] = {left_ghost, right_ghost};
    MPI_Ineighbor_alltoallv(&u[0], counts, sdispls, MPI_DOUBLE,
            &u[0], counts, rdispls, MPI_DOUBLE, COMM_GRID, &neighbor);
}

////////////////////////////////////////////////////////////////////////////////

/* Write result to file. */
void writeresult(const int N_local, const double* u) {
    /* Global index of processor's first cell */
    const int global_offset = world_rank * ((N - 2) / world_size) + MIN(world_rank, (N - 2) % world_size) + 1;
    /* Length of output line */
    const int LINE_SIZE = 19;
    /* Allocate memory for every output line and \0 terminator */
    char * buf = malloc(N_local * LINE_SIZE + 1);
    /* Write output to buffer */
    for (int i = 1; i <= N_local; ++i) {
        const double x_n = (global_offset + i) * h;
        sprintf(buf + (i - 1) * LINE_SIZE, "%f %+f\n", x_n, u[i]);
    }

    /* MPI IO */
    MPI_File fh;
    MPI_File_open(MPI_COMM_WORLD, "pde.out", MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
    if (world_rank == 0) {
        /* Master writes boundary condition, since this is not contained in
         any local cell */
        MPI_File_seek(fh, 0, MPI_SEEK_SET);
        char * boundery0 = "0.000000 +0.000000\n";
        MPI_File_write(fh, boundery0, LINE_SIZE, MPI_CHAR, MPI_STATUS_IGNORE);
    }
    MPI_File_seek(fh, global_offset * LINE_SIZE, MPI_SEEK_SET);
    MPI_File_write(fh, buf, (N_local * LINE_SIZE), MPI_CHAR, MPI_STATUS_IGNORE);
    MPI_File_close(&fh);
    free(buf);
}
\$\endgroup\$
  • 2
    \$\begingroup\$ I would agree with your suspicion on the NUMA. I think it is pretty rare to see an OpenMP benefit beyond a single NUMA domain. Also, the advantage of a hybrid parallelization often only gets visible for large node counts, so a single node measurement might be a little misleading. Still, 4 OMP threads per MPI Rank and 8 ranks per node, seem to be a good option on your system. So I'd conclude your assessment is correct . \$\endgroup\$ – haraldkl Sep 24 '15 at 16:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.