#include "../single-thread/nbody.c" // Include original code
#include "graphics.h" // For Graphics_Run function only
+
+#ifdef OVER_ENGINEERED
+System * sub_system = NULL;
+#endif //OVER_ENGINEERED
+
/**
* @function Simulation_Run
* @purpose Start the simulation
*/
void Simulation_Run(int argc, char ** argv)
{
- omp_set_nested(1);
- if (omp_get_nested() != 1)
- {
- fprintf(stderr, "Couldn't set nested parallelism, I kind of need that\n");
- exit(EXIT_FAILURE);
- }
-
-
+
if (options.num_threads == 0)
options.num_threads = omp_get_max_threads();
if (options.draw_graphics)
- omp_set_num_threads(2);
- else
- omp_set_num_threads(1);
-
-
-
-
- #pragma omp parallel sections
- {
-
- #pragma omp section
+ {
+ Graphics_Run(argc, argv);
+ if (options.pedantic_graphics == false)
{
- omp_set_num_threads(options.num_threads);
- while (true)
+ options.num_threads += 1;
+ omp_set_nested(1);
+ if (omp_get_nested() == 0)
{
-
- if (runstate != RUN) break; //Check whether the main thread wants to exit
-
- if (ExitCondition())
- {
- break;
- }
+ fprintf(stderr, "Couldn't enable nested parallelism\n.");
+ exit(EXIT_FAILURE);
+ }
- if (options.draw_graphics == false && options.verbosity != 0
- && universe.steps % options.verbosity == 1)
+ #pragma omp parallel sections num_threads(2)
+ {
+ #pragma omp section
{
- DisplayStatistics();
+ Compute();
+ }
+ #pragma omp section
+ {
+ glutMainLoop();
}
-
-
- Compute();
}
+ return;
}
+ }
+ Compute();
+
- #pragma omp section
- {
- if (options.draw_graphics)
- Graphics_Run(argc, argv);
- printf("Got to bit after Graphics_Run()\n");
- }
- }
}
-/**
- * @function Compute
- * @purpose Compute a single step, multithreaded using OpenMP
- */
-void Compute(void)
+
+void Compute(void)
{
- //printf("There are %d threads working and %d available\n", omp_get_num_threads(), omp_get_max_threads());
- //return;
- unsigned n_per_thread = (universe.N) / options.num_threads;
- unsigned remainder = (universe.N) % options.num_threads;
- #pragma omp parallel for
- for (unsigned i = 0; i < options.num_threads; ++i)
+ omp_set_num_threads(options.num_threads);
+ if (omp_get_max_threads() != options.num_threads)
{
- //printf("Thread %d executes iteration %u\n", omp_get_thread_num(), i);
- System s;
- s.body = universe.body + (i * n_per_thread);
- s.N = (i == options.num_threads - 1) ? n_per_thread : n_per_thread + remainder;
- System_Forces(&s, &universe);
+ fprintf(stderr, "Couldn't use %d threads!\n", options.num_threads);
+ exit(EXIT_FAILURE);
}
- #pragma omp parallel for
- for (unsigned i = 0; i < options.num_threads; ++i)
+ #ifdef OVER_ENGINEERED
+ sub_system = Split_System(&universe, options.num_threads); // I could use omp for loops instead, but I like this way
+ #endif //OVER_ENGINEERED
+
+ #pragma omp parallel
{
- //printf("Thread %d executes iteration %u\n", omp_get_thread_num(), i);
- System s;
- s.body = universe.body + (i * n_per_thread);
- s.N = (i == options.num_threads - 1) ? n_per_thread : n_per_thread + remainder;
- System_Positions(&s);
- }
+ #ifdef OVER_ENGINEERED
+ int index = omp_get_thread_num();
+ #endif //OVER_ENGINEERED
+ while (!ExitCondition())
+ {
+
+ #ifdef OVER_ENGINEERED
+ System_Forces(sub_system+index, &universe);
+ #pragma omp barrier // Look! An explicit barrier!
+ System_Positions(sub_system+index);
+ #pragma omp barrier
+ #else
+ #pragma omp for
+ for (unsigned a = 0; a < universe.N; ++a)
+ {
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ (universe.body+a)->F[i] = 0;
+ //printf("Thread %d computes force for body %d\n", omp_get_thread_num(), a);
+ omp_set_num_threads(options.nested_threads);
+ #pragma omp parallel for
+ for (unsigned b = 0; b < universe.N; ++b)
+ {
+ if (b != a)
+ Body_Force(universe.body+a, universe.body+b);
+ }
+ }
+
+ #pragma omp for
+ for (unsigned a = 0; a < universe.N; ++a)
+ {
+ //printf("Thread %d computes position for body %d\n", omp_get_thread_num(), a);
+ Body_Velocity(universe.body+a);
+ Body_Position(universe.body+a);
+ }
- universe.steps += 1;
+ #endif //OVER_ENGINEERED
+ #pragma omp master
+ {
+ // The if statement has a cost, but if can be removed easily in an optimised version.
+ // Only the main thread is allowed to mess with glut, because glut is evil (not thread safe).
+ // Although this function may be nested, that will only happen if the if would be false.
+ // ("if the if would be false"; I am awesome at explaining code)
+ if (options.draw_graphics && options.pedantic_graphics)
+ {
+ glutMainLoopEvent();
+ Graphics_Display();
+ }
+ }
+ #pragma omp single
+ {
+ // Only one thread should execute this, and we don't care which one.
+ // So much easier than the pthread version!
+ universe.steps += 1;
+ if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
+ DisplayStatistics();
+ }
+
+ }
+ }
}
+