- #pragma omp parallel for
- for (unsigned i = 0; i < options.num_threads; ++i)
- {
- //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
+ 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);
+ }