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 omp_lock_t runstate_lock;
  17
  18 /**
  19  * @function Simulation_Run
  20  * @purpose Start the simulation
  21  * @param argc, argv - Passed to Graphics_Run if needed, otherwise unused
  22  */
  23 void Simulation_Run(int argc, char ** argv)
  24 {
  25         omp_init_lock(&runstate_lock);
  26         if (options.num_threads == 0)
  27                 options.num_threads = omp_get_max_threads();
  28
  29         if (options.draw_graphics)
  30         {
  31                 Graphics_Run(argc, argv);
  32                 if (options.pedantic_graphics == false)
  33                 {
  34                         options.num_threads += 1;
  35                         omp_set_nested(1);
  36                         if (omp_get_nested() == 0)
  37                         {
  38                                 fprintf(stderr, "Couldn't enable nested parallelism\n.");
  39                                 exit(EXIT_FAILURE);
  40                         }
  41
  42
  43                         #pragma omp parallel num_threads(2)
  44                         {
  45                                 // This can't be done with sections!
  46                                 // Because glut is useless, and can only be dealt with in the main thread
  47                                 // I just hope that thread "0" is always the main thread!
  48                                 if (omp_get_thread_num() == 0) // #pragma omp section
  49                                 {
  50                                         //printf("Thread %d gets graphics\n", omp_get_thread_num());
  51                                         glutMainLoop();
  52                                         QuitProgram(false); // Program gets to here if the window is quit
  53                                 }
  54                                 else // #pragma omp section
  55                                 {
  56                                         //printf("Thread %d gets computations\n", omp_get_thread_num());
  57                                         Compute();
  58                                 }
  59                         }
  60                         return;
  61                 }
  62         }
  63         Compute();
  64
  65
  66         omp_destroy_lock(&runstate_lock);
  67 }
  68
  69
  70 void Compute(void)
  71 {
  72         omp_set_num_threads(options.num_threads);
  73         if (omp_get_max_threads() != options.num_threads)
  74         {
  75                 fprintf(stderr, "Couldn't use %d threads!\n", options.num_threads);
  76                 exit(EXIT_FAILURE);
  77         }
  78
  79         #ifdef OVER_ENGINEERED
  80         sub_system = Split_System(&universe, options.num_threads); // I could use omp for loops instead, but I like this way
  81         #endif //OVER_ENGINEERED
  82
  83         #pragma omp parallel
  84         {
  85                 #ifdef OVER_ENGINEERED
  86                 int index = omp_get_thread_num();
  87                 #endif //OVER_ENGINEERED
  88                 while (!ExitCondition())
  89                 {
  90
  91                         #ifdef OVER_ENGINEERED
  92                                 System_Forces(sub_system+index, &universe);
  93                                 #pragma omp barrier // Look! An explicit barrier!
  94                                 System_Positions(sub_system+index);
  95                                 #pragma omp barrier
  96                         #else
  97                                 #pragma omp for
  98                                 for (unsigned a = 0; a < universe.N; ++a)
  99                                 {
 100                                         for (unsigned i = 0; i < DIMENSIONS; ++i)
 101                                                 (universe.body+a)->F[i] = 0;
 102                                         //printf("Thread %d computes force for body %d\n", omp_get_thread_num(), a);
 103                                         omp_set_num_threads(options.nested_threads);
 104                                         #pragma omp parallel for
 105                                         for (unsigned b = 0; b < universe.N; ++b)
 106                                         {
 107                                                 if (b != a)
 108                                                         Body_Force(universe.body+a, universe.body+b);
 109                                         }
 110                                 }
 111
 112                                 #pragma omp for
 113                                 for (unsigned a = 0; a < universe.N; ++a)
 114                                 {
 115                                         //printf("Thread %d computes position for body %d\n", omp_get_thread_num(), a);
 116                                         Body_Velocity(universe.body+a);
 117                                         Body_Position(universe.body+a);
 118                                 }
 119
 120                         #endif //OVER_ENGINEERED
 121
 122                         #pragma omp master
 123                         {
 124                                 // The if statement has a cost, but if can be removed easily in an optimised version.
 125                                 // Only the main thread is allowed to mess with glut, because glut is evil (not thread safe).
 126                                 // Although this function may be nested, that will only happen if the if would be false.
 127                                 // ("if the if would be false"; I am awesome at explaining code)
 128                                 if (options.draw_graphics && options.pedantic_graphics)
 129                                 {
 130                                         glutMainLoopEvent();
 131                                         Graphics_Display();
 132                                 }
 133                         }
 134                         #pragma omp single
 135                         {
 136                                 // Only one thread should execute this, and we don't care which one.
 137                                 // So much easier than the pthread version!
 138                                 universe.steps += 1;
 139                                 if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
 140                                         DisplayStatistics();
 141                         }
 142
 143                 }
 144         }
 145         QuitProgram(false);
 146 }
 147
 148 void QuitProgram(bool error)
 149 {
 150         if (runstate == QUIT || runstate == QUIT_ERROR)
 151                 return;
 152         omp_set_lock(&runstate_lock);
 153         runstate = (error) ? QUIT_ERROR : QUIT;
 154         omp_unset_lock(&runstate_lock);
 155 }
 156