course/semester2/pprog/assignment1/openmp/nbody.c

   1 /**
   2  * @file nbody.c
   3  * @purpose OpenMP version of N-Body simulator, implementation
   4  * @author Sam Moore (20503628) - 2012
   5  */
   6
   7 #include "nbody.h"
   8 #include "../single-thread/nbody.c" // Include original code
   9 #include "graphics.h" // For Graphics_Run function only
  10
  11
  12 #ifdef OVER_ENGINEERED
  13 System * sub_system = NULL;
  14 #endif //OVER_ENGINEERED
  15
  16 /**
  17  * @function Simulation_Run
  18  * @purpose Start the simulation
  19  * @param argc, argv - Passed to Graphics_Run if needed, otherwise unused
  20  */
  21 void Simulation_Run(int argc, char ** argv)
  22 {
  23
  24         if (options.num_threads == 0)
  25                 options.num_threads = omp_get_max_threads();
  26
  27         if (options.draw_graphics)
  28         {
  29                 Graphics_Run(argc, argv);
  30                 if (options.pedantic_graphics == false)
  31                 {
  32                         options.num_threads += 1;
  33                         omp_set_nested(1);
  34                         if (omp_get_nested() == 0)
  35                         {
  36                                 fprintf(stderr, "Couldn't enable nested parallelism\n.");
  37                                 exit(EXIT_FAILURE);
  38                         }
  39
  40
  41                         #pragma omp parallel sections num_threads(2)
  42                         {
  43                                 #pragma omp section
  44                                 {
  45                                         Compute();
  46                                 }
  47                                 #pragma omp section
  48                                 {
  49                                         glutMainLoop();
  50                                 }
  51                         }
  52                         return;
  53                 }
  54         }
  55         Compute();
  56
  57
  58
  59 }
  60
  61
  62 void Compute(void)
  63 {
  64         omp_set_num_threads(options.num_threads);
  65         if (omp_get_max_threads() != options.num_threads)
  66         {
  67                 fprintf(stderr, "Couldn't use %d threads!\n", options.num_threads);
  68                 exit(EXIT_FAILURE);
  69         }
  70
  71         #ifdef OVER_ENGINEERED
  72         sub_system = Split_System(&universe, options.num_threads); // I could use omp for loops instead, but I like this way
  73         #endif //OVER_ENGINEERED
  74
  75         #pragma omp parallel
  76         {
  77                 #ifdef OVER_ENGINEERED
  78                 int index = omp_get_thread_num();
  79                 #endif //OVER_ENGINEERED
  80                 while (!ExitCondition())
  81                 {
  82
  83                         #ifdef OVER_ENGINEERED
  84                                 System_Forces(sub_system+index, &universe);
  85                                 #pragma omp barrier // Look! An explicit barrier!
  86                                 System_Positions(sub_system+index);
  87                                 #pragma omp barrier
  88                         #else
  89                                 #pragma omp for
  90                                 for (unsigned a = 0; a < universe.N; ++a)
  91                                 {
  92                                         for (unsigned i = 0; i < DIMENSIONS; ++i)
  93                                                 (universe.body+a)->F[i] = 0;
  94                                         //printf("Thread %d computes force for body %d\n", omp_get_thread_num(), a);
  95                                         omp_set_num_threads(options.nested_threads);
  96                                         #pragma omp parallel for
  97                                         for (unsigned b = 0; b < universe.N; ++b)
  98                                         {
  99                                                 if (b != a)
 100                                                         Body_Force(universe.body+a, universe.body+b);
 101                                         }
 102                                 }
 103
 104                                 #pragma omp for
 105                                 for (unsigned a = 0; a < universe.N; ++a)
 106                                 {
 107                                         //printf("Thread %d computes position for body %d\n", omp_get_thread_num(), a);
 108                                         Body_Velocity(universe.body+a);
 109                                         Body_Position(universe.body+a);
 110                                 }
 111
 112                         #endif //OVER_ENGINEERED
 113
 114                         #pragma omp master
 115                         {
 116                                 // The if statement has a cost, but if can be removed easily in an optimised version.
 117                                 // Only the main thread is allowed to mess with glut, because glut is evil (not thread safe).
 118                                 // Although this function may be nested, that will only happen if the if would be false.
 119                                 // ("if the if would be false"; I am awesome at explaining code)
 120                                 if (options.draw_graphics && options.pedantic_graphics)
 121                                 {
 122                                         glutMainLoopEvent();
 123                                         Graphics_Display();
 124                                 }
 125                         }
 126                         #pragma omp single
 127                         {
 128                                 // Only one thread should execute this, and we don't care which one.
 129                                 // So much easier than the pthread version!
 130                                 universe.steps += 1;
 131                                 if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
 132                                         DisplayStatistics();
 133                         }
 134
 135                 }
 136         }
 137 }
 138