main.c

#include <limits.h>
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include "configuration.h"
#include "mpicomm.h"
#include "mesh.h"
#include "field.h"

#ifdef DEBUG_ATTACH
#include <unistd.h>
void wait_for_debugger();
#endif

void print_decomposition(const conf_t *configuraiton, const box_t *decomposition, FILE *f);
void print_adjacency_matrix(const mesh_t *mesh, const conf_t *configuraiton, FILE *f);

void halo_exchange_demo(const mesh_t *mesh, const conf_t *configuration);

int main(int argc, char *argv[])
{
   int mpi_rank, nprocs;
   const char* error;
   mesh_t mesh;

   MPI_Init(&argc, &argv);

   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
   MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);

   // initialize random number generator
   srand((mpi_rank + 17) * (int)time(NULL));

   conf_t configuration;
   conf_init(&configuration);

   if(mpi_rank == 0) {
      // parse comand line args on root process only
      error = conf_set_from_args(&configuration, argc, argv);
      if(error) {
         fprintf(stderr, "%s\n", error);
         MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
      }
   }

   // Broadcast configuration (nonblocking)
   MPI_Request req_bcast_conf;
   mpi_broadcast_configuration_start(&configuration, &req_bcast_conf);

   // Create and broadcast decomposition
   box_t *decomposition = malloc(nprocs * sizeof(*decomposition));
   if(mpi_rank == 0) {
      const box_t *cube = &configuration.global_domain;
      box_decompose(cube, nprocs, decomposition);
   }
   mpi_broadcast_decomposition(decomposition);

   // Finish nonblocking configuration broadcast
   mpi_broadcast_configuration_finish(&req_bcast_conf);

   if(info_enabled(&configuration)) {
      conf_print(&configuration, stdout);
      print_decomposition(&configuration, decomposition, stdout);
   }

#ifdef DEBUG_ATTACH
   if(mpi_rank == configuration.debug_rank) {
      wait_for_debugger();
   }
#endif

   mesh_init(&mesh, &configuration, decomposition);
   print_adjacency_matrix(&mesh, &configuration, stdout);

   halo_exchange_demo(&mesh, &configuration);

   mesh_free(&mesh);

   free(decomposition);
   MPI_Finalize();

   return EXIT_SUCCESS;
}

// --------------------------------------------------------------------------

void halo_exchange_demo(const mesh_t *mesh, const conf_t *configuration)
{
   int mpi_rank, i, n;
   int_field_t field;
   const neighbor_info *nb;

   int_field_init(&field, configuration, mesh);

   MPI_Comm_rank(mesh->communicator, &mpi_rank);

   // Set field data and exchange halo information
   n = configuration->n_iterations;
   if(info_enabled(configuration)) {
      printf("Exchanging halo information (%d iterations)...\n\n", n);
   }

   for(i = 0; i < n; i++) {
      if(trace_enabled(configuration)) {
         printf("Iteration %d\n", i);
      }
      int_field_set_constant(&field, mpi_rank);
      int_field_halo_exchange(&field);
   }

   if(info_enabled(configuration)) {
      printf("Validating...\n\n");
   }

   // Validate local field data
   if(!int_field_check_value(&field, &mesh->local_domain, mpi_rank)) {
      fprintf(stderr, "ERROR (rank %d): Local cell data corrupted.\n", mpi_rank);
      MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
   }

   // Validate halo data received from each neighbor
   for(i = 0; i < mesh->n_neighbors; i++) {
      nb = &mesh->neighbors[i];
      if(!int_field_check_value(&field, &nb->halo_incoming.domain, nb->mpi_rank)) {
         fprintf(stderr, "ERROR (rank %d): Halo values for neighbor rank %d not correct.\n",
                 mpi_rank, nb->mpi_rank);
         MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
      }
   }

   if(info_enabled(configuration)) {
      printf("Validation successful.\n\n");
   }

   int_field_free(&field);
}

void print_decomposition(const conf_t *configuration, const box_t *boxes, FILE *f)
{
   int nprocs;
   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);

   long volume, min_volume = LONG_MAX, max_volume = 0;
   int i, j;
   for(i = 0; i < nprocs; i++) {
      volume = boxes[i].volume;
      if(min_volume > volume) min_volume = boxes[i].volume;
      if(max_volume < volume) max_volume = boxes[i].volume;
   }
   fprintf(f, "Cells per processor (min-max): %ld - %ld\n\n", min_volume, max_volume);

   if(trace_enabled(configuration)) {
      char buf[128];
      box_t intersection;
      fprintf(f, "Decomposition:\n");
      for(i = 0; i < nprocs; i++) {
         fprintf(f, " * box %5d: %s\n", i, box_to_string(&boxes[i], buf, sizeof(buf))); 

         for(j = 0; j < i; j++) {
            box_intersect(&boxes[i], &boxes[j], &intersection);
            if(!box_is_empty(&intersection)){
               fprintf(f, "   * intersects with box %d\n", j);
            }
         }
      }
      fprintf(f, "\n");
   }
}

void print_adjacency_matrix(const mesh_t *mesh, const conf_t *configuration, FILE *f)
{
   if(configuration->verbosity_level >= DEBUG)
   {
      int mpi_rank, nprocs;
      int *sendbuf, *recvbuf, *row;
      neighbor_info *nb;
      int i, j, min, max, n;

      MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
      MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);

      sendbuf = calloc(nprocs, sizeof(*sendbuf));
      recvbuf = mpi_rank == 0 ? calloc(nprocs * nprocs, sizeof(*recvbuf)) : NULL;

      for(i = 0; i < mesh->n_neighbors; i++) {
         nb = &mesh->neighbors[i];
         sendbuf[nb->mpi_rank] = nb->halo_incoming.domain.volume;
      }

      MPI_Gather(
            sendbuf, nprocs, MPI_INT,
            recvbuf, nprocs, MPI_INT,
            0, MPI_COMM_WORLD);

      if(mpi_rank == 0) {
         fprintf(f, "Adjacency matrix:\n");
         fprintf(f, "%2s ", "");
         for(i = 0; i < nprocs; i++) {
            fprintf(f, "%2d ", i);
         }
         fprintf(f, "\n");

         for(j = 0; j < nprocs; j++) {
            row = recvbuf + j * nprocs;

            fprintf(f, "%2d ", j);
            for(i = 0; i < nprocs; i++) {
               fprintf(f, "%2s ", row[i] ? "X" : "");
            }
            fprintf(f, "\n");
         }
         fprintf(f, "\n");

         min = INT_MAX; max = 0;
         for(i = 0; i < nprocs * nprocs; i++) {
            n = recvbuf[i];
            if(n > 0 && min > n) min = recvbuf[i];
            if(max < n) max = n;
         }
         fprintf(f, "Number of adjacent cells per neighbor (min-max): %d - %d\n\n", min, max);
         free(recvbuf);
      }
      free(sendbuf);
   }
}

#ifdef DEBUG_ATTACH
void wait_for_debugger()
{
   int rank;
   volatile int resume = 0;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);

   printf("Rank %d (pid %d) waiting for debugger to attach...\n", rank, getpid());
   while(!resume) {
      sleep(1);
   }
}
#endif