From bb7fa31ea517a1fba064e723b37d5b8d8bd7dd72 Mon Sep 17 00:00:00 2001 From: Sam Moore Date: Sun, 30 Sep 2012 23:56:02 +0800 Subject: [PATCH] Parallel Programming - Finished OpenMP I think --- .../pprog/assignment1/mthread/nbody.c | 33 +- .../pprog/assignment1/mthread/nbody.c~ | 505 ++++++++++++++++++ .../pprog/assignment1/mthread/nbody.h | 2 +- .../pprog/assignment1/mthread/nbody.h~ | 50 ++ .../pprog/assignment1/openmp/graphics.c~ | 388 ++++++++++++++ .../pprog/assignment1/openmp/nbody.c | 199 ++++--- .../pprog/assignment1/openmp/nbody.c~ | 142 +++++ .../pprog/assignment1/openmp/nbody.h | 22 +- .../pprog/assignment1/openmp/nbody.h~ | 48 ++ .../assignment1/single-thread/graphics.c | 8 +- .../pprog/assignment1/single-thread/nbody.c | 72 ++- .../pprog/assignment1/single-thread/nbody.c~ | 305 +++++++++++ .../pprog/assignment1/single-thread/nbody.h | 5 +- .../pprog/assignment1/single-thread/nbody.h~ | 140 +++++ 14 files changed, 1755 insertions(+), 164 deletions(-) create mode 100644 course/semester2/pprog/assignment1/mthread/nbody.c~ create mode 100644 course/semester2/pprog/assignment1/mthread/nbody.h~ create mode 100644 course/semester2/pprog/assignment1/openmp/graphics.c~ create mode 100644 course/semester2/pprog/assignment1/openmp/nbody.c~ create mode 100644 course/semester2/pprog/assignment1/openmp/nbody.h~ create mode 100644 course/semester2/pprog/assignment1/single-thread/nbody.c~ create mode 100644 course/semester2/pprog/assignment1/single-thread/nbody.h~ diff --git a/course/semester2/pprog/assignment1/mthread/nbody.c b/course/semester2/pprog/assignment1/mthread/nbody.c index 13831cd1..a551d7d5 100644 --- a/course/semester2/pprog/assignment1/mthread/nbody.c +++ b/course/semester2/pprog/assignment1/mthread/nbody.c @@ -460,38 +460,7 @@ void Simulation_Run(int argc, char ** argv) -/** - * @function Split_System - * @purpose Helper to divide one system into an array of systems - * Each sub system will have N = (s->N / n) bodies in it - * @param s - The original system (typically &universe) - * @param n - The number of sub systems in the array - * - * WARNING: It is the caller's responsibility to free() the returned array - */ -System * Split_System(System * s, unsigned n) -{ - System * result = (System*)(calloc(n, sizeof(System))); - if (result == NULL) - { - perror("Couldn't create array of sub systems"); - QuitProgram(true); - pthread_exit(NULL); - } - unsigned n_per_system = (s->N) / n; - unsigned remainder = (s->N) % n; - - for (unsigned i = 0; i < n; ++i) - { - result[i].N = n_per_system; - if (i == n-1) - result[i].N += remainder; - result[i].body = (s->body) + (n_per_system * i); - result[i].steps = 0; - } - return result; -} /** * @function Allocate_Threads @@ -527,7 +496,7 @@ void * StepFunction(void * arg) System * s = (System*)(arg); s->steps += 1; //Increment number of steps computed - if (options.verbosity != 0 && s->steps % options.verbosity == 1) + if (options.verbosity != 0 && s->steps % options.verbosity == 0) DisplayStatistics(); diff --git a/course/semester2/pprog/assignment1/mthread/nbody.c~ b/course/semester2/pprog/assignment1/mthread/nbody.c~ new file mode 100644 index 00000000..03def4b5 --- /dev/null +++ b/course/semester2/pprog/assignment1/mthread/nbody.c~ @@ -0,0 +1,505 @@ +/** + * @file nbody.c + * @purpose Implentation of multi-threaded N-Body simulator using pthreads + * @author Sam Moore (20503628) - 2012 + */ + +#include "nbody.h" // Declarations +#include "../single-thread/nbody.c" // Include all functions from the single threaded version +#include +#include "graphics.h" // For declaration of Graphics_Run only +#include "barrier.c" + +// --- Variable declarations --- // + +pthread_t compute_thread; // The thread responsible for computations; it spawns worker threads (* terms and conditions apply) + +pthread_t * worker_thread = NULL; //Array of worker threads responsible for Force and Position updates +System * sub_system = NULL; //Array of Systems used to divide up the main "universe" System for worker threads +System * nested_sub_system = NULL; //Array of Systems for division of "universe" for the nested worker threads + +pthread_mutex_t mutex_runstate; // Mutex around the runstate + + + + +Barrier force_barrier; // I laughed at this variable name. A bit sad really. +Barrier position_barrier; + + +Barrier graphics_barrier; + +pthread_mutex_t mutex_threads_running; +int threads_running = 0; + +/** + * @function Compute_Thread + * @purpose Thread - Continuously computes steps for a system of bodies. Seperate from graphics functions. + * Spawns worker threads to divide up computation. + * @param arg - Can be cast to the System for which computations are performed + * + * NOTE: + * This will always be (void*)(&universe) where universe is the global System variable that has every Body in it. + * But I don't like global variables. And since the argument had to be passed, I thought I might as well use it. + * That way, when I change "universe" to "solar_system", I only have to change the argument where this is called, not all through it. + * Find and replace? Who uses that!? + */ + + +void * Compute_Thread(void * arg) +{ + + System * s = (System*)(arg); //cast argument to a System* + + + + // If no number of threads provided, use the default value, unless someone (me) changed that to a stupid value + if (options.num_threads <= 0) + options.num_threads = (DEFAULT_WORKING_THREADS > 1) ? DEFAULT_WORKING_THREADS : 1; // Fear the ternary operator! + + if (options.nested_threads <= 0) + options.nested_threads = 1; + + // Do a sanity check; there is no point spawning more threads than bodies. + if (options.num_threads > s->N) + { + fprintf(stderr, + "Warning: Using %u threads instead of %u specified, because there are only %u bodies to simulate!\n", + s->N, options.num_threads, s->N); + options.num_threads = s->N; + } + + // Initialise the barriers ("shields up!") + Barrier_Init(&force_barrier, options.num_threads); + Barrier_Init(&position_barrier, options.num_threads); + Barrier_Init(&graphics_barrier, 1); + + + + if (options.num_threads > 1) // If we require additional worker threads... + { + // Allocate worker threads and sub systems + + sub_system = Split_System(&universe, options.num_threads); + + if (options.nested_threads > 1) + nested_sub_system = Split_System(&universe, options.nested_threads); + + #ifdef PERSISTENT_THREADS // Code for the smart way of doing it (spawn threads once, keep them running) + worker_thread = Allocate_Threads(options.num_threads-1); + // Spawn a bunch of worker threads, and let them all do their thing. + // Note the "-1". Because this thread can do work too! + for (unsigned i = 0; i < options.num_threads-1; ++i) + { + if (pthread_create(worker_thread+i, NULL, Worker_Thread, (void*)(sub_system+i)) != 0) + { + perror("In compute thread, couldn't create worker thread"); + QuitProgram(true); + pthread_exit(NULL); + } + } + Worker_Thread((void*)(sub_system+options.num_threads-1)); // This thread becomes a worker thread + #else + worker_thread = Allocate_Threads(options.num_threads); + #endif //PERSISTENT_THREADS + + } + #ifdef PERSISTENT_THREADS + else // We only require one worker thread... + { + // So just do all computations in this thread + while (!ExitCondition()) + { + + System_Forces(s, s); + // If required, wait for graphics to finish drawing + if (options.draw_graphics && options.pedantic_graphics) + Barrier_Wait(&graphics_barrier); + System_Positions(s); + StepFunction(s); + + if (options.draw_graphics && options.pedantic_graphics) + Barrier_ForceExit(&position_barrier); //Make the graphics continue + } + QuitProgram(false); + pthread_exit(NULL); + } + + #else // Code for the stupid way of doing it (respawn threads each step) + // (ie: The way I immediately implemented and didn't realise was stupid until someone told me) + + + + + // Run until we can't run anymore + while (!ExitCondition()) + { + + + if (options.num_threads <= 1) // If there is only 1 worker thread... + { + // Just do everything in this thread + System_Forces(s, s); + // If required, wait for graphics to finish drawing + if (options.draw_graphics && options.pedantic_graphics) + Barrier_Join(&graphics_barrier); + + System_Positions(s); + StepFunction(s); + + if (options.draw_graphics && options.pedantic_graphics) + Barrier_Join(&position_barrier); //Make the graphics continue + continue; + } + + + + //Compute forces by spawning threads, each thread gets a sub system + for (unsigned i = 0; i < options.num_threads; ++i) + { + if (pthread_create(worker_thread+i, NULL, Force_Thread, (void*)(sub_system+i)) != 0) + { + perror("In compute thread, couldn't create worker thread (force)"); + QuitProgram(true); + pthread_exit(NULL); + } + + } + + for (unsigned i = 0; i < options.num_threads; ++i) + pthread_join(worker_thread[i], NULL); + + + // If required, wait for graphics to finish drawing + if (options.draw_graphics && options.pedantic_graphics) + Barrier_Wait(&graphics_barrier); + + + Barrier_Enter(&position_barrier); + + //Compute positions by spawning a bunch of threads to do it + for (unsigned i = 0; i < options.num_threads; ++i) + { + if (pthread_create(worker_thread+i, NULL, Position_Thread, (void*)(sub_system+i)) != 0) + { + perror("In compute thread, couldn't create worker thread (position)"); + QuitProgram(true); + pthread_exit(NULL); + } + } + + + for (unsigned i = 0; i < options.num_threads; ++i) + pthread_join(worker_thread[i], NULL); + + + StepFunction(s); // Execute single threaded stuff + + + } + QuitProgram(false); + #endif //PERSISTENT_THREADS + return NULL; +} + +/** + * @function BeforeDraw + * @purpose Called in graphics thread before the draw loop + * When --pedantic-graphics enabled, will wait for position computations to finish before drawing + * Otherwise does nothing + * + * This originally seemed like a good place to put the code now in StepFunction(), since only one thread runs this + * But then I realised that the graphics might be disabled, + * and there was no point having a thread that only existed to call that code. + * + * So I changed it to the horrible solution that I currently have. + */ +void BeforeDraw() +{ + + //printf("BEFORE DRAW\n"); + if (!options.pedantic_graphics) + return; + + //printf("Graphics thread waits on position barrier\n"); + Barrier_Wait(&position_barrier); + //printf("\tGraphics thread wakes up\n"); + Barrier_Enter(&graphics_barrier); +} + +/** + * @function AfterDraw + * @purpose Called in graphics thread after the draw loop + * When --pedantic-graphics is supplied, will signal computation threads that drawing is finished + * So that positions can be safely altered + * Otherwise does nothing + */ +void AfterDraw() +{ + //universe.steps += 1; + if (!options.pedantic_graphics) + return; + Barrier_Join(&graphics_barrier); + +} + +#ifdef PERSISTENT_THREADS + +/** + * @function Worker_Thread + * @purpose Thread - A self contained worker thread to compute a particular sub system of bodies + * + * This is the "smart" way to do it, because threads are only created once, and compute both force and position. + * The greatest difficulty with pthreads is getting a *single* thread from the team to execute certain code + * (ie: The stuff in StepFunction()). + * With the "continuously respawning threads of stupidity" approach, + * because there is one "master" thread (not necessarilly the main thread... don't get confused now) + * to keep respawning the workers, the single threaded code can just be executed in the master thread. + * + * With this approach, I have created a hacky solution so that the *last* thread to leave the position barrier gets to call StepFunction. + * + */ +void * Worker_Thread(void * arg) +{ + System * s = (System*)(arg); // This is mainly to save typing the RHS a lot of times + + // Each thread runs until the whole program is supposed to end + while (!ExitCondition()) + { + + + System_Forces(s, &universe); // Each thread computes the forces for its share of bodies + + // Do not confuse with "Barrier_Wait". + // Barrier_Wait does not affect the barrier; it just waits for it + // Barrier_Join actively updates the state of the barrier, and wakes up sleeping threads if required. + + Barrier_Join(&force_barrier); // All threads must reach here before moving on. + if (ExitCondition()) return NULL; + + + //fprintf(stderr,"Thread %p - force barrier finished\n", arg); + //printf("Computed ALL forces\n"); + + + // If required, wait for the graphics to finish drawing stuff + if (options.draw_graphics && options.pedantic_graphics) + { + //printf("Worker %p waits on graphics barrier\n", arg); + Barrier_Wait(&graphics_barrier); + //printf("\tWorker %p wakes up after graphics barrier\n", arg); + if (ExitCondition()) return NULL; + } + + + + Barrier_Enter(&position_barrier); + System_Positions(s); // Each thread updates the positions for its share of bodies + + + // Barrier_JoinCall behaves in the same way as Barrier_Join, except the *last* thread + // (ie: the one that wakes up the others) also calls the function with arguments given. + Barrier_JoinCall(&position_barrier, StepFunction, (void*)(&universe)); + if (ExitCondition()) return NULL; + //Barrier_Join(&position_barrier); + + // All threads have computed positions, and *one* thread calls StepFunction() + + } + QuitProgram(false); // Set the run state of the program + return NULL; +} + +#endif //PERSISTENT_THREADS + +/** + * @function Force_Thread + * @purpose Thread - Calculates the forces on objects in a System + * @param s - Can be cast to the System* + */ +void * Force_Thread(void * s) +{ + //System_ForcePair * pair = (System_ForcePair*)s; + + + System_Forces(s, &universe); //Simple wrapper + //printf("Force_Thread waits\n"); + + return NULL; +} + + +/** + * @function Position_Thread + * @purpose Thread - Calculates the positions of objects in a System + * @param s - Can be cast to the System* + */ +void * Position_Thread(void * s) +{ + + System_Positions((System*)s); // Simple wrapper + Barrier_Join(&position_barrier); // This needed so that graphics will wait + + return NULL; +} + +/** + * @function QuitProgram + * @purpose This function can be called in any thread to signal all threads to exit + * Repeated calls to this function have no effect + * + * All threads periodically call ExitCondition(), which will return true if the program should exit. + * One (not the only way) to return true is if this function has been called. + * Threads will call this function if they detect ExitCondition() is true. Only the first call has any effect. + */ +inline void QuitProgram(bool error) +{ + //If already quitting, don't do anything + if (runstate == QUIT || runstate == QUIT_ERROR) + return; + + + + // set the runstate (checked in ExitCondition()) + + pthread_mutex_lock(&mutex_runstate); // aquire mutex + if (error) + runstate = QUIT_ERROR; // Program is exiting due to an error + else + runstate = QUIT; // Program is exiting naturally + pthread_mutex_unlock(&mutex_runstate); //release mutex + + +} + +/** + * @function Thread_Cleanup + * @purpose Will be called in the *main* thread when exit() is called + * Ensures working threads will exit, and waits for them to finish. + * Then waits for them to finish. + * Also frees memory associated with the worker threads. + */ +void Thread_Cleanup(void) +{ + + + // Threads recheck the exit condition whenever they leave a barrier. + // These calls will stop any threads waiting forever in a barrier for threads that exited before getting to the barrier. + Barrier_ForceExit(&force_barrier); + Barrier_ForceExit(&position_barrier); + + + if (options.draw_graphics) // If the graphics are enabled... + { + // Then there is a computation thread, since graphics are done in the main thread + pthread_join(compute_thread, NULL); + } + + #ifdef PERSISTENT_THREADS + for (unsigned i = 0; i < options.num_threads-1; ++i) + { + pthread_join(worker_thread[i], NULL); + } + #else + // All other worker threads (if they were spawned) are terminated in Compute_Thread + #endif //PERSISTENT_THREADS + + // Scary memory management here. + if (worker_thread != NULL) + free(worker_thread); + if (sub_system != NULL) + free(sub_system); + worker_thread = NULL; + sub_system = NULL; + +} + + +/** + * @function Simulation_Run + * @purpose Initialise and start the simulation. Will be called in the main thread. + * Replaces the single-threaded macro that does nothing, and sets up the graphics and computation threads + * @param argc - Number of arguments - Passed to Graphics_Run if needed + * @param argv - Argument strings - Passed to Graphics_Run if needed + */ +void Simulation_Run(int argc, char ** argv) +{ + atexit(Thread_Cleanup); + if (options.draw_graphics) // The graphics are enabled + { + // I have chosen to do graphics in the main thread in this case. + // A *single* seperate thread is spawned here to do computations. + // This computation thread will spawn any additional worker threads required. + if (pthread_create(&compute_thread, NULL, Compute_Thread, (void*)&universe) != 0) + { + perror("Error creating compute thread"); + exit(EXIT_FAILURE); + } + + // This is run in the main thread + // It is effectively the graphics initialisation, followed by the glut loop + Graphics_Run(argc, argv); + + // The main thread reaches here after leaving the glut loop when ExitCondition() returns true. + + QuitProgram(false); + + exit(EXIT_SUCCESS); // This is the main thread; use exit() + + } + else //The graphics are disabled + { + // If graphics are disabled, there is no point spawning an extra thread. + // In this case, the *main* thread starts computations. + // Note that it will probably spawn additional worker threads (unless options.num_threads <= 1) + Compute_Thread((void*)(&universe)); + QuitProgram(false); + exit(EXIT_SUCCESS); + } +} + + + + + +/** + * @function Allocate_Threads + * @purpose Helper function to allocate an array of pthread_t objects + * Handles all the pointless, er, "important" error checking that should be done + * @param n - Number of threads in the array + * + * WARNING: Remember to free() the array!!! + */ +pthread_t * Allocate_Threads(unsigned n) +{ + pthread_t * result = (pthread_t*)(calloc(n, sizeof(pthread_t))); + if (result == NULL) + { + perror("Unable to allocate memory for threads"); + QuitProgram(true); + pthread_exit(NULL); + } + return result; +} + +/** + * @function StepFunction + * @purpose Helper to perform stuff in a single thread every step, after position computations are done + * The reason this has void* all over the place is so that I can pass the function pointer (to horrible dragons and fiendish demons). + * @param arg - Can be cast to System* for which steps are to be updated + * Will always be (void*)(&universe). But I have been brainwashed into the "global variables are baaaaad" philosophy. + * @returns arg + */ +void * StepFunction(void * arg) +{ + //fprintf(stderr, "StepFunction called\n"); + System * s = (System*)(arg); + s->steps += 1; //Increment number of steps computed + + if (options.verbosity != 0 && s->steps % options.verbosity == 1) + DisplayStatistics(); + + + return arg; +} + diff --git a/course/semester2/pprog/assignment1/mthread/nbody.h b/course/semester2/pprog/assignment1/mthread/nbody.h index 9948e524..375c2c94 100644 --- a/course/semester2/pprog/assignment1/mthread/nbody.h +++ b/course/semester2/pprog/assignment1/mthread/nbody.h @@ -32,7 +32,7 @@ void * Compute_Thread(void * system); //Thread - Continuously perform computatio -System * Split_System(System * s, unsigned n); // Splits one system into a number of other systems, returns an array of size n + pthread_t * Allocate_Threads(unsigned n); // Allocates space for threads - handles errors #ifdef PERSISTENT_THREADS diff --git a/course/semester2/pprog/assignment1/mthread/nbody.h~ b/course/semester2/pprog/assignment1/mthread/nbody.h~ new file mode 100644 index 00000000..9948e524 --- /dev/null +++ b/course/semester2/pprog/assignment1/mthread/nbody.h~ @@ -0,0 +1,50 @@ +#ifndef _NBODY_MTHREAD_H +#define _NBODY_MTHREAD_H + + +#include "../single-thread/nbody.h" //Use original simulation code +#include +#include "barrier.h" // Barriers + +#undef SINGLE_THREADED +#define PTHREADED + +#define DEFAULT_WORKING_THREADS 2 + +//#define PERSISTENT_THREADS //If defined, threads will not be continually destroyed and then respawned + + + +//Undefine default macros, replace with functions +#undef Simulation_Run +void Simulation_Run(int argc, char ** argv); +#undef QuitProgram +void QuitProgram(bool error); + +#undef BeforeDraw +void BeforeDraw(); +#undef AfterDraw +void AfterDraw(); + + +void * Compute_Thread(void * system); //Thread - Continuously perform computations for a System of bodies. May spawn additional worker threads. + + + + +System * Split_System(System * s, unsigned n); // Splits one system into a number of other systems, returns an array of size n +pthread_t * Allocate_Threads(unsigned n); // Allocates space for threads - handles errors + +#ifdef PERSISTENT_THREADS +void * Worker_Thread(void * arg); +#endif //PERSISTENT_THREADS +void * Force_Thread(void * system); //Thread - Compute forces for all objects in a system +void * Position_Thread(void * system); //Thread - Compute positions for all objects in a system + +void Thread_Cleanup(void); //Called at program exit to safely join computation thread + + + +void * StepFunction(void * s); //Function called at the end of every step + +#endif //_NBODY_MTHREAD_H diff --git a/course/semester2/pprog/assignment1/openmp/graphics.c~ b/course/semester2/pprog/assignment1/openmp/graphics.c~ new file mode 100644 index 00000000..627cea12 --- /dev/null +++ b/course/semester2/pprog/assignment1/openmp/graphics.c~ @@ -0,0 +1,388 @@ +/** + * @file graphics.c + * @author Sam Moore (20503628) 2012 + * - Adapted from template program provided by UWA + * @purpose Definition of graphics related functions + * NOTE: This file is identical for both the single-threaded and multi-threaded versions of the program + */ + + +#include +#include +#include // Needed for clock +#include // Maths functions (sin, cos) + +#include "graphics.h" //Function declarations +#include "nbody.h" //The simulation + + +// --- Variable definitions --- // +double previousTime, eyeTheta, eyePhi, eyeRho; +float look[3]; +int windowWidth, windowHeight, upY; +double scale = 1.0; + +#ifdef FLYING_CAMERA +Camera eye; +#endif //FLYING_CAMERA + +/** + * @function Graphics_Run + * @purpose Setup and start the graphics engine + * @param argc - Number of arguments to main function, passed to glutInit + * @param argv - Argument strings for main function, passed to glutInit + */ +void Graphics_Run(int argc, char ** argv) +{ + if (options.draw_graphics == false) + { + // This message is left here for when I inevitably accidentally call the function + fprintf(stderr, "Graphics_Run should not be called if graphics are disabled!\n"); + exit(EXIT_FAILURE); + } + + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(WIDTH, HEIGHT); + glutInitWindowPosition(POSITION_X, POSITION_Y); + + //Set window title based on version of program + #ifdef SINGLE_THREADED + glutCreateWindow("N-Body : Single-Threaded"); + #elif defined PTHREADED + glutCreateWindow("N-Body : Multi-Threaded (pthread)"); + #elif defined OMP_THREADED + glutCreateWindow("N-Body : Multi-Threaded (OpenMP)"); + #else + glutCreateWindow("N-Body"); + #endif + glutDisplayFunc(Graphics_Display); + glutIdleFunc(Graphics_Display); + glutKeyboardFunc(Graphics_Keyboard); + glutReshapeFunc(Graphics_Reshape); + + + glClearColor(1.0,1.0,1.0,0.0); + glColor3f(0.0f, 0.0f, 0.0f); + glPointSize(POINT_SIZE); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + + /*init lighting */ + + GLfloat mat_specular[] = { 1.0, 1.0, 1.0, 1.0 }; + GLfloat mat_shininess[] = { 50.0 }; + GLfloat light_position[] = { 1.0, 1.0, 0.0, 0.0 }; + glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular); + glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess); + glLightfv(GL_LIGHT0, GL_POSITION, light_position); + + glColorMaterial(GL_FRONT,GL_DIFFUSE); // Set Color Capability + + glEnable(GL_LIGHTING); + glEnable(GL_LIGHT0); + glEnable(GL_DEPTH_TEST); + + glEnable(GL_COLOR_MATERIAL); // Enable color + + + windowWidth = WIDTH; + windowHeight = HEIGHT; + previousTime = clock(); + eyeTheta = 0; + eyePhi = 0.5 * M_PI; + eyeRho = RHO; + upY = 1; + look[0] = 0; + look[1] = 0; + look[2] = 0; + gluPerspective(VIEW_ANGLE, (double)WIDTH/(double)HEIGHT, WORLD_NEAR, WORLD_FAR); + + #ifdef FLYING_CAMERA + eye.x[0] = 0; eye.x[1] = 0; eye.x[2] = 1; + eye.y[0] = 0; eye.y[1] = 1; eye.y[2] = 0; + eye.z[0] = 1; eye.z[1] = 0; eye.z[2] = 0; + eye.p[0] = 0; eye.p[1] = 0; eye.p[2] = -50; + #endif //FLYING_CAMERA + + + //This option must be set, or when glut leaves the main loop. the exit(3) function is called... annoying + glutSetOption(GLUT_ACTION_ON_WINDOW_CLOSE, GLUT_ACTION_CONTINUE_EXECUTION); + + + #ifndef OMP_THREADED + glutMainLoop(); + #endif //OMP_THREADED +} + +/** + * @function Graphics_Display + * @purpose This function redraws the screen after the positions of particles + * have been updated. + * It also calls System_Compute, and checks for exit conditions, in the single-threaded version only + */ +void Graphics_Display() +{ + + if (options.draw_graphics == false) + { + // This message is left here for when I inevitably accidentally call the function + fprintf(stderr, "Graphics_Display should not be called if graphics are disabled!\n"); + exit(EXIT_FAILURE); + } + + //Check whether the graphics thread should exit for any reason + if (ExitCondition()) + { + glutLeaveMainLoop(); // original glut does not have this, which makes "nicely" exiting a bit annoying + return; + } + + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); + + #ifdef FLYING_CAMERA + + gluLookAt(eye.p[0], eye.p[1], eye.p[2], + eye.p[0] + eye.x[0], eye.p[1] + eye.x[1], eye.p[2] + eye.x[2], + eye.z[0], eye.z[1], eye.z[2]); + #else + gluLookAt(eyeRho * sin(eyePhi) * sin(eyeTheta), eyeRho * cos(eyePhi), + eyeRho * sin(eyePhi) * cos(eyeTheta),look[0], look[1], look[2], 0, upY, 0); + #endif //FLYING_CAMERA + + BeforeDraw(); // Stuff to do before graphics is allowed to draw + // Single-thread - perform a computation step, obviously! It's not done anywhere else + // Multiple threads - ensure that all body positions are updated (ie: not halfway through step). + // (We don't care about the force and velocity variables though) + + for (int i = 0; i < universe.N; ++i) + { + + Body * b = universe.body+i; + glColor3f(0.0f, b->mass/1e11*100, 0.0f); + //glColor3f(1.0f, 0.0f, 0.0f); + glPushMatrix(); // to save the current matrix + glTranslated(scale*b->x[0], scale*b->x[1], scale*b->x[2]); + glutSolidSphere (BALL_SIZE, 10, 10); + glPopMatrix(); // restore the previous matrix + } + glFlush(); + + AfterDraw(); // Stuff to do after graphics is done drawing + // Single-thread - Nothing at all + // Multiple threads - signal threads it is safe to change position variables + +} + +/** + * @function Graphics_Keyboard + * @purpose This function is to manipulate with the image + * @param theKey key pressed + * @param mouseX, mouseY position of the mouse in the window + */ +void Graphics_Keyboard(unsigned char theKey, int mouseX, int mouseY) +{ + if (theKey == 'x' || theKey == 'X') + { + + QuitProgram(false); + + return; + } + + #ifdef FLYING_CAMERA + switch (theKey) + { + // Translate in direction of camera + case 'W': + case 'w': + for (unsigned i = 0; i < DIMENSIONS; ++i) + eye.p[i] += eye.x[i]; + break; + // Translate backwards from camera direction + case 'S': + case 's': + for (unsigned i = 0; i < DIMENSIONS; ++i) + eye.p[i] -= eye.x[i]; + break; + // Translate left from camera direction + case 'A': + case 'a': + for (unsigned i = 0; i < DIMENSIONS; ++i) + eye.p[i] -= eye.y[i]; + break; + // Translate right from camera direction + case 'D': + case 'd': + for (unsigned i = 0; i < DIMENSIONS; ++i) + eye.p[i] += eye.y[i]; + break; + // Translate up from camera direction + case 'Q': + case 'q': + for (unsigned i = 0; i < DIMENSIONS; ++i) + eye.p[i] += eye.z[i]; + break; + // Translate down from camera direction + case 'E': + case 'e': + for (unsigned i = 0; i < DIMENSIONS; ++i) + eye.p[i] -= eye.z[i]; + break; + + // Rotate camera direction "down" + // (pitch control) + case 'I': + case 'i': + { + float theta = M_PI/80.0; + Camera old = eye; + + for (unsigned i = 0; i < DIMENSIONS; ++i) + { + eye.z[i] = old.z[i] * cos(theta) + old.x[i] * sin(theta); + eye.x[i] = old.x[i] * cos(theta) - old.z[i] * sin(theta); + } + break; + } + // Rotate camera direction "up" + // (pitch control) + case 'K': + case 'k': + { + float theta = -M_PI/80.0; + Camera old = eye; + + for (unsigned i = 0; i < DIMENSIONS; ++i) + { + eye.z[i] = old.z[i] * cos(theta) + old.x[i] * sin(theta); + eye.x[i] = old.x[i] * cos(theta) - old.z[i] * sin(theta); + } + break; + } + // Rotate camera direction "left" + // (yaw control) + case 'J': + case 'j': + { + float theta = +M_PI/80.0; + Camera old = eye; + + for (unsigned i = 0; i < DIMENSIONS; ++i) + { + eye.y[i] = old.y[i] * cos(theta) + old.x[i] * sin(theta); + eye.x[i] = old.x[i] * cos(theta) - old.y[i] * sin(theta); + } + break; + } + // Rotate camera direction "right" + // (yaw control) + case 'L': + case 'l': + { + float theta = -M_PI/80.0; + Camera old = eye; + + for (unsigned i = 0; i < DIMENSIONS; ++i) + { + eye.y[i] = old.y[i] * cos(theta) + old.x[i] * sin(theta); + eye.x[i] = old.x[i] * cos(theta) - old.y[i] * sin(theta); + } + break; + } + // Rotate camera direction CCW about its axis + // (roll control) + case 'U': + case 'u': + { + float theta = -M_PI/80.0; + Camera old = eye; + + for (unsigned i = 0; i < DIMENSIONS; ++i) + { + eye.z[i] = old.z[i] * cos(theta) + old.y[i] * sin(theta); + eye.y[i] = old.y[i] * cos(theta) - old.z[i] * sin(theta); + } + break; + } + // Rotate camera direction CW about its axis + // (roll control) + case 'O': + case 'o': + { + float theta = +M_PI/80.0; + Camera old = eye; + + for (unsigned i = 0; i < DIMENSIONS; ++i) + { + eye.z[i] = old.z[i] * cos(theta) + old.y[i] * sin(theta); + eye.y[i] = old.y[i] * cos(theta) - old.z[i] * sin(theta); + } + break; + } + } + + /* Code used for debugging the camera + system("clear"); + printf("Camera status:\n"); + printf("Position: %f %f %f\n", eye.p[0], eye.p[1], eye.p[2]); + printf("x: %f %f %f (%f)\n", eye.x[0], eye.x[1], eye.x[2], sqrt(eye.x[0]*eye.x[0] + eye.x[1]*eye.x[1] + eye.x[2]*eye.x[2])); + printf("y: %f %f %f (%f)\n", eye.y[0], eye.y[1], eye.y[2], sqrt(eye.y[0]*eye.y[0] + eye.y[1]*eye.y[1] + eye.y[2]*eye.y[2])); + printf("z: %f %f %f (%f)\n", eye.z[0], eye.z[1], eye.z[2], sqrt(eye.z[0]*eye.z[0] + eye.z[1]*eye.z[1] + eye.z[2]*eye.z[2])); + */ + #else // The original view code + + // I like how 'I' is used three times. + + if (theKey == 'i' || theKey == 'I') { + eyePhi -= M_PI / 20; + if (eyePhi == 0) + eyePhi = 2 * M_PI; + } else if (theKey == 'm' || theKey == 'I') { + eyePhi += M_PI / 20; + } else if (theKey == 'j' || theKey == 'J') { + eyeTheta -= M_PI / 20; + } else if (theKey == 'k' || theKey == 'K') { + eyeTheta += M_PI / 20; + } else if (theKey == ',') { + eyeRho += 0.5; + } else if (theKey == '.' || theKey == 'I') { + eyeRho -= 0.5; + } else if (theKey == 'w' || theKey == 'W') { + look[1] += 0.5; + } else if (theKey == 'z' || theKey == 'Z') { + look[1] -= 0.5; + } else if (theKey == 'a' || theKey == 'A') { + look[0] -= 0.5; + } else if (theKey == 's' || theKey == 'S') { + look[0] += 0.5; + } else if (theKey == '+') { + scale *= 1.1; + } else if (theKey == '-') { + scale *= 0.9; + } + if (sin(eyePhi) > 0) upY = 1; + else upY = 1; + #endif //FLYING_CAMERA +} + +/** + * @function Graphics_Reshape + * @purpose Resize the view window + * @param width, height - New size of the window + */ +void Graphics_Reshape(int width, int height) +{ + double displayRatio = 1.0 * width / height; + windowWidth = width; + windowHeight = height; + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluPerspective(VIEW_ANGLE, displayRatio, WORLD_NEAR, WORLD_FAR); + glViewport(0, 0, windowWidth, windowHeight); +} + + + diff --git a/course/semester2/pprog/assignment1/openmp/nbody.c b/course/semester2/pprog/assignment1/openmp/nbody.c index 1d2f3690..053625bb 100644 --- a/course/semester2/pprog/assignment1/openmp/nbody.c +++ b/course/semester2/pprog/assignment1/openmp/nbody.c @@ -8,8 +8,10 @@ #include "../single-thread/nbody.c" // Include original code #include "graphics.h" // For Graphics_Run function only -omp_lock_t graphics_lock; -bool graphics_busy = false; + +#ifdef OVER_ENGINEERED +System * sub_system = NULL; +#endif //OVER_ENGINEERED /** * @function Simulation_Run @@ -18,138 +20,119 @@ bool graphics_busy = false; */ 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); - + { + Graphics_Run(argc, argv); + if (options.pedantic_graphics == false) + { + options.num_threads += 1; + omp_set_nested(1); + if (omp_get_nested() == 0) + { + fprintf(stderr, "Couldn't enable nested parallelism\n."); + exit(EXIT_FAILURE); + } - omp_init_lock(&graphics_lock); - #pragma omp parallel sections - { - - #pragma omp section - { - omp_set_num_threads(options.num_threads); - while (true) + #pragma omp parallel sections num_threads(2) { - - if (runstate != RUN) break; //Check whether the main thread wants to exit - - if (ExitCondition()) + #pragma omp section { - break; - } - - if (options.verbosity != 0 && universe.steps % options.verbosity == 1) + Compute(); + } + #pragma omp section { - DisplayStatistics(); + glutMainLoop(); } - - - Compute(); } - //printf("Left compute loop\n"); + return; } + } + Compute(); - #pragma omp section - { - if (options.draw_graphics) - Graphics_Run(argc, argv); - //printf("Got to bit after Graphics_Run()\n"); - } - } - omp_destroy_lock(&graphics_lock); } -/** - * @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); } - // Wait for graphics to finish drawing - - if (options.draw_graphics && options.pedantic_graphics) - { - while (graphics_busy); - } + #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 for - for (unsigned i = 0; i < options.num_threads; ++i) + #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); - } - - universe.steps += 1; + #ifdef OVER_ENGINEERED + int index = omp_get_thread_num(); + #endif //OVER_ENGINEERED + while (!ExitCondition()) + { - // Tell graphics to continue - if (options.draw_graphics && options.pedantic_graphics) - { - omp_set_lock(&graphics_lock); - graphics_busy = true; - omp_unset_lock(&graphics_lock); + #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); + } -} + #endif //OVER_ENGINEERED -/** - * @function BeforeDraw - * @purpose Called before the graphics thread draws bodies. - * If pedantic graphics are used, will wait until graphics_busy is set to true before drawing - */ -void BeforeDraw() -{ - if (options.verbosity != 0 && universe.steps % options.verbosity == 0) - DisplayStatistics(); - if (!options.pedantic_graphics) - return; - while (graphics_busy == false); + #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(); + } + + } + } } -/** - * @function AfterDraw - * @purpose Called after the graphics thread draws bodies - * If pedantic graphics used, will unset graphics_busy when done drawing - */ -void AfterDraw() -{ - if (!options.pedantic_graphics) - return; - omp_set_lock(&graphics_lock); - graphics_busy = false; - omp_unset_lock(&graphics_lock); -} diff --git a/course/semester2/pprog/assignment1/openmp/nbody.c~ b/course/semester2/pprog/assignment1/openmp/nbody.c~ new file mode 100644 index 00000000..0b6d160f --- /dev/null +++ b/course/semester2/pprog/assignment1/openmp/nbody.c~ @@ -0,0 +1,142 @@ +/** + * @file nbody.c + * @purpose OpenMP version of N-Body simulator, implementation + * @author Sam Moore (20503628) - 2012 + */ + +#include "nbody.h" +#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 + * @param argc, argv - Passed to Graphics_Run if needed, otherwise unused + */ +void Simulation_Run(int argc, char ** argv) +{ + + if (options.num_threads == 0) + options.num_threads = omp_get_max_threads(); + + if (options.draw_graphics) + { + Graphics_Run(argc, argv); + if (options.pedantic_graphics == false) + { + options.num_threads += 1; + omp_set_nested(1); + if (omp_get_nested() == 0) + { + fprintf(stderr, "Couldn't enable nested parallelism\n."); + exit(EXIT_FAILURE); + } + + + #pragma omp parallel sections num_threads(2) + { + #pragma omp section + { + Compute(); + } + #pragma omp section + { + glutMainLoop(); + } + } + return; + } + } + #ifdef OVER_ENGINEERED + Compute(); + #else + Compute2(); + #endif //OVER_ENGINEERED + + + +} + + +void Compute(void) +{ + omp_set_num_threads(options.num_threads); + if (omp_get_max_threads() != options.num_threads) + { + fprintf(stderr, "Couldn't use %d threads!\n", options.num_threads); + exit(EXIT_FAILURE); + } + + #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 + { + #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); + } + + #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(); + } + + } + } +} + diff --git a/course/semester2/pprog/assignment1/openmp/nbody.h b/course/semester2/pprog/assignment1/openmp/nbody.h index 3e3d2718..aea97bfc 100644 --- a/course/semester2/pprog/assignment1/openmp/nbody.h +++ b/course/semester2/pprog/assignment1/openmp/nbody.h @@ -21,14 +21,26 @@ #undef Simulation_Run void Simulation_Run(int argc, char ** argv); //#undef QuitProgram -//void QuitProgram(bool error); +//#define QuitProgram(x) break #undef BeforeDraw -void BeforeDraw(); -#undef AfterDraw -void AfterDraw(); +#define BeforeDraw() (void)0 // Do nothing (apparently this is how to do nothing with a macro) +//void BeforeDraw(); -void Compute(); //Compute a single step +#undef AfterDraw +#define AfterDraw() (void)0 +//void AfterDraw(); + + + +void Compute(void); +void Compute2(void); + +//#define OVER_ENGINEERED +// define to (my original approach) manually split the System array (same as pthread version) between threads +// Manually splitting the System array is not required in openmp, since "#pragma omp for" automatically does it +// I still like doing it this way, since you can explicitly see what each thread is doing, but it would probably lose me marks. +// However I can't bring myself to delete it. Hence, compromise. #endif //_NBODY_OPENMP_H diff --git a/course/semester2/pprog/assignment1/openmp/nbody.h~ b/course/semester2/pprog/assignment1/openmp/nbody.h~ new file mode 100644 index 00000000..ebab8927 --- /dev/null +++ b/course/semester2/pprog/assignment1/openmp/nbody.h~ @@ -0,0 +1,48 @@ +#ifndef _NBODY_OPENMP_H +#define _NBODY_OPENMP_H + +/** + * @file nbody.h + * @purpose OpenMP version of N-Body simulator, declarations + * @author Sam Moore (20503628) - 2012 + */ + +#include "../single-thread/nbody.h" //Include original code +#include + + +#undef SINGLE_THREADED +#define OMP_THREADED + + +#define CRAPPY_VERSION + +// Replace default macros with thread-safe functions +#undef Simulation_Run +void Simulation_Run(int argc, char ** argv); +//#undef QuitProgram +//#define QuitProgram(x) break + +#undef BeforeDraw +#define BeforeDraw() (void)0 // Do nothing (apparently this is how to do nothing with a macro) +//void BeforeDraw(); + +#undef AfterDraw +#define AfterDraw() (void)0 +//void AfterDraw(); + +omp_lock_t graphics_lock; + + +void Compute(void); +void Compute2(void); + +//#define OVER_ENGINEERED +// define to (my original approach) manually split the System array (same as pthread version) between threads +// Manually splitting the System array is not required in openmp, since "#pragma omp for" automatically does it +// I still like doing it this way, since you can explicitly see what each thread is doing, but it would probably lose me marks. +// However I can't bring myself to delete it. Hence, compromise. + +#endif //_NBODY_OPENMP_H + +//EOF diff --git a/course/semester2/pprog/assignment1/single-thread/graphics.c b/course/semester2/pprog/assignment1/single-thread/graphics.c index c3d55f52..e960ffe6 100644 --- a/course/semester2/pprog/assignment1/single-thread/graphics.c +++ b/course/semester2/pprog/assignment1/single-thread/graphics.c @@ -109,8 +109,10 @@ void Graphics_Run(int argc, char ** argv) //This option must be set, or when glut leaves the main loop. the exit(3) function is called... annoying glutSetOption(GLUT_ACTION_ON_WINDOW_CLOSE, GLUT_ACTION_CONTINUE_EXECUTION); - glutMainLoop(); - + + #ifndef OMP_THREADED + glutMainLoop(); + #endif //OMP_THREADED } /** @@ -184,7 +186,9 @@ void Graphics_Keyboard(unsigned char theKey, int mouseX, int mouseY) { if (theKey == 'x' || theKey == 'X') { + QuitProgram(false); + glutLeaveMainLoop(); return; } diff --git a/course/semester2/pprog/assignment1/single-thread/nbody.c b/course/semester2/pprog/assignment1/single-thread/nbody.c index ee510ed4..30e5c761 100644 --- a/course/semester2/pprog/assignment1/single-thread/nbody.c +++ b/course/semester2/pprog/assignment1/single-thread/nbody.c @@ -37,11 +37,9 @@ inline void Body_Print(Body * a, FILE * out) * @param a - The Body * @param b - The System */ -inline void Body_Force(Body * a, System * s) +inline void Body_Forces(Body * a, System * s) { - double distance; - double con; - double gd; + for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero a->F[i] = 0; @@ -52,18 +50,32 @@ inline void Body_Force(Body * a, System * s) if (b == a) continue; //Otherwise we would have infinite force - //Calculate distance between a and b - distance = 0.0; - for (unsigned i = 0; i < DIMENSIONS; ++i) - distance += square(b->x[i] - a->x[i]); - distance = sqrt(distance); - con = G * a->mass * b->mass / square(distance); - gd = con / distance; - for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force - a->F[i] += gd * (b->x[i] - a->x[i]); + Body_Force(a, b); } } +/** + * @function Body_Force + * @purpose Compute force on body a due to body b + * @param a - Body a + * @param b - Body b + */ +inline void Body_Force(Body * a, Body * b) +{ + double distance; + double con; + double gd; + //Calculate distance between a and b + distance = 0.0; + for (unsigned i = 0; i < DIMENSIONS; ++i) + distance += square(b->x[i] - a->x[i]); + distance = sqrt(distance); + con = G * a->mass * b->mass / square(distance); + gd = con / distance; + for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force + a->F[i] += gd * (b->x[i] - a->x[i]); +} + /** * @function Body_Velocity * @purpose Compute velocity of a body @@ -110,7 +122,7 @@ inline void System_Forces(System * s1, System * s2) //printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N); for (unsigned i = 0; i < s1->N; ++i) { - Body_Force(s1->body+i, s2); + Body_Forces(s1->body+i, s2); Body_Velocity(s1->body+i); } //printf("Compute positions for %p - Done!\n", (void*)s1); @@ -262,3 +274,35 @@ inline bool ExitCondition(void) //fprintf(stderr,"runstate %d\n timeout %d\n steps %d\n", (int)(runstate != RUN), (int)(options.timeout > 0.00 && ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)), (int)(options.num_steps > 0 && universe.steps > options.num_steps)); return result; } + +/** + * @function Split_System + * @purpose Helper to divide one system into an array of systems + * Each sub system will have N = (s->N / n) bodies in it + * @param s - The original system (typically &universe) + * @param n - The number of sub systems in the array + * + * WARNING: It is the caller's responsibility to free() the returned array + */ +System * Split_System(System * s, unsigned n) +{ + System * result = (System*)(calloc(n, sizeof(System))); + if (result == NULL) + { + perror("Couldn't create array of sub systems"); + exit(EXIT_FAILURE); + } + + unsigned n_per_system = (s->N) / n; + unsigned remainder = (s->N) % n; + + for (unsigned i = 0; i < n; ++i) + { + result[i].N = n_per_system; + if (i == n-1) + result[i].N += remainder; + result[i].body = (s->body) + (n_per_system * i); + result[i].steps = 0; + } + return result; +} diff --git a/course/semester2/pprog/assignment1/single-thread/nbody.c~ b/course/semester2/pprog/assignment1/single-thread/nbody.c~ new file mode 100644 index 00000000..df7cd9d0 --- /dev/null +++ b/course/semester2/pprog/assignment1/single-thread/nbody.c~ @@ -0,0 +1,305 @@ +/** + * @file nbody.c + * @author Sam Moore (20503628) 2012 + * @purpose N-Body simulator - Definition of simulation functions; single threaded version + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "nbody.h" + +RUNSTATE runstate = RUN; + +/** + * @function Body_Print + * @purpose Print the body to a file (or stdout) + * @param a - Body to print' + * @param out - FILE* to print to + */ +inline void Body_Print(Body * a, FILE * out) +{ + //fprintf(out,"Body %p M=%f X=%f Y=%f Z=%f Fx=%f Fy=%f Fz=%f Vx=%f Vy=%f Vz=%f\n", + //(void*)a, a->mass, a->x[0], a->x[1], a->x[2], a->F[0], a->F[1], a->F[2], a->v[0], a->v[1], a->v[2]); + fprintf(out, "%f %f %f %f %f %f %f %f %f %f\n", + a->mass, a->x[0], a->x[1], a->x[2], a->v[0], a->v[1], a->v[2], a->F[0], a->F[1], a->F[2]); +} + +/** + * @function Body_Force + * @purpose Calculates the force on a single Body due to all bodies in a System + * @param a - The Body + * @param b - The System + */ +inline void Body_Forces(Body * a, System * s) +{ + + + for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero + a->F[i] = 0; + + for (unsigned index = 0; index < s->N; ++index) // Step through all bodies in the System + { + Body * b = s->body+index; //Current body + if (b == a) + continue; //Otherwise we would have infinite force + + Body_Force(a, b); + } +} + +/** + * @function Body_Force + * @purpose Compute force on body a due to body b + * @param a - Body a + * @param b - Body b + */ +inline void Body_Force(Body * a, Body * b) +{ + //Calculate distance between a and b + distance = 0.0; + for (unsigned i = 0; i < DIMENSIONS; ++i) + distance += square(b->x[i] - a->x[i]); + distance = sqrt(distance); + con = G * a->mass * b->mass / square(distance); + gd = con / distance; + for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force + a->F[i] += gd * (b->x[i] - a->x[i]); +} + +/** + * @function Body_Velocity + * @purpose Compute velocity of a body + * @param a - The Body + */ +inline void Body_Velocity(Body * a) +{ + for (unsigned i = 0; i < DIMENSIONS; ++i) + a->v[i] += a->F[i] / a->mass * DELTA_T; +} + +/** + * @function Body_Position + * @purpose Compute position of a body + * @param a - The Body + */ +inline void Body_Position(Body * a) +{ + for (unsigned i = 0; i < DIMENSIONS; ++i) + a->x[i] += a->v[i] * DELTA_T; +} + +/** + * @function System_Compute + * @purpose Compute forces and then positions for bodies in System + * NOTE: Only used in the single threaded version of the program + */ +inline void System_Compute(System * s) +{ + System_Forces(s, s); + System_Positions(s); +} + +/** + * @function System_Forces + * @purpose Compute forces on each object in System s1 due to all bodies in System s2 + * @param s1, s2 - The two systems + * NOTE: In single threaded version, called with s1 == s2 == &universe + * In multi threaded version, called with s2 == &universe, but s1 is constructed for each thread + * In nested multi-threaded version, s2 is constructed for the nested threads as well + */ +inline void System_Forces(System * s1, System * s2) +{ + //printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N); + for (unsigned i = 0; i < s1->N; ++i) + { + Body_Forces(s1->body+i, s2); + Body_Velocity(s1->body+i); + } + //printf("Compute positions for %p - Done!\n", (void*)s1); +} + +/** + * @function System_Positions + * @purpose Compute positions of all objects in System + * @param s - The system + */ +inline void System_Positions(System * s) +{ + //printf("Compute positions for %p - %d bodies...\n", (void*)s, s->N); + for (unsigned i = 0; i < s->N; ++i) + Body_Position(s->body+i); + //printf("Compute positions for %p - Done!\n", (void*)s); + +} + + + + +/** + * @function System_Init + * @purpose Initialise a System from an input file + * @param s - The System + * @param fileName - The input file + */ +inline void System_Init(System * s, const char *fileName) +{ + char line[LINE_SIZE]; + char * token; + FILE * file; + + file = fopen(fileName, "rt"); + s->N = atoi(fgets(line, LINE_SIZE, file)); + s->body = (Body*) calloc((size_t)s->N, sizeof(Body)); + s->steps = 0; + + //printf("----------------------Initial field-------------------------------\n"); + + for (unsigned i = 0; i < s->N; ++i) + { + if (fgets(line, LINE_SIZE, file) != NULL) + { + Body * a = s->body+i; + token = strtok(line, ",; \t"); + a->mass = atof(token); + + for (unsigned j = 0; j < DIMENSIONS; ++j) + { + token = strtok(NULL, ",; \t"); + a->x[j] = atof(token); + } + for (unsigned j = 0; j < DIMENSIONS; ++j) + { + token = strtok(NULL, ",; \t"); + a->v[j] = atof(token); + } + + // Ignore any additional tokens + while (token != NULL) + token = strtok(NULL, ",; \t"); + //Body_Print(a); + } + } + + //printf("--------------------------------------------------------------\n"); + + fclose(file); +} + +/** + * @function Universe_Cleanup + * @purpose Cleans up the universe by freeing all memory + * Also prints the bodies in the universe to a file specified in "options" if required. + */ +inline void Universe_Cleanup() +{ + //fprintf(stderr, "Called Universe_Cleanup()\n"); + if (options.output != NULL) // An output file was specified... + { + FILE * save = NULL; + if (strcmp("stdout", options.output) == 0) + save = stdout; + else if (strcmp("stderr", options.output) == 0) + save = stderr; + else + save = fopen(options.output, "w"); + if (save == NULL) + perror("Couldn't save universe to file."); + else + { + // Print the output in the same format as the initial field file + fprintf(save, "%u\n", universe.N); + + for (unsigned i = 0; i < universe.N; ++i) // Print all bodies to the file + { + Body_Print(universe.body+i, save); + + + } + fclose(save); + } + + } + free(universe.body); // Cleanup memory used for bodies + +} + +/** + * @function DisplayStatistics + * @purpose Display useful information about the computations done, + * - Total number of steps computed + * - Time in seconds elapsed + * - Clock cycles executed + * - Equivelant time for a single thread to execute same number of clock cycles + */ +inline void DisplayStatistics() +{ + clock_t end = clock(); + struct timeval end_time; + if (gettimeofday(&end_time, NULL) != 0) + { + perror("Couldn't get time of day.\n"); + return; + } + float runtime = (float)(end_time.tv_sec-options.start_time.tv_sec) + (float)(end_time.tv_usec-options.start_time.tv_usec) / 1.0e6; + float clock_time = (float)(end - options.start_clock) / (float)(CLOCKS_PER_SEC); + printf("%u\t%f\t%u\t%f\n", universe.steps, runtime, (unsigned)(end - options.start_clock), clock_time); +} + + +/** + * @function ExitCondition + * @purpose Helper to check whether the program is supposed to exit + * Does not check for user triggered quit + * Checks for either a timeout, or computation of the required number of steps + * + * The reason timeouts are integers is because this function is called a lot in lots of threads + * So using floats (and working out microseconds every time) is expensive + * + */ +inline bool ExitCondition(void) +{ + bool result = (runstate != RUN || + (options.num_steps >= 0 && universe.steps >= options.num_steps) || + (options.timeout >= 0 && ((int)(time(NULL) - options.start_time.tv_sec) >= options.timeout))); + //fprintf(stderr,"runstate %d\n timeout %d\n steps %d\n", (int)(runstate != RUN), (int)(options.timeout > 0.00 && ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)), (int)(options.num_steps > 0 && universe.steps > options.num_steps)); + return result; +} + +/** + * @function Split_System + * @purpose Helper to divide one system into an array of systems + * Each sub system will have N = (s->N / n) bodies in it + * @param s - The original system (typically &universe) + * @param n - The number of sub systems in the array + * + * WARNING: It is the caller's responsibility to free() the returned array + */ +System * Split_System(System * s, unsigned n) +{ + System * result = (System*)(calloc(n, sizeof(System))); + if (result == NULL) + { + perror("Couldn't create array of sub systems"); + exit(EXIT_FAILURE); + } + + unsigned n_per_system = (s->N) / n; + unsigned remainder = (s->N) % n; + + for (unsigned i = 0; i < n; ++i) + { + result[i].N = n_per_system; + if (i == n-1) + result[i].N += remainder; + result[i].body = (s->body) + (n_per_system * i); + result[i].steps = 0; + } + return result; +} diff --git a/course/semester2/pprog/assignment1/single-thread/nbody.h b/course/semester2/pprog/assignment1/single-thread/nbody.h index a9bfec44..2ce28b05 100644 --- a/course/semester2/pprog/assignment1/single-thread/nbody.h +++ b/course/semester2/pprog/assignment1/single-thread/nbody.h @@ -108,7 +108,8 @@ typedef struct } Options; void Body_Print(Body * a, FILE * out); //Print body a -void Body_Force(Body * a, System * s); //Compute force on body a due to system of bodies s +void Body_Forces(Body * a, System * s); //Compute force on body a due to system of bodies s +void Body_Force(Body * a, Body * b); // Compute force on body a due to body b void Body_Velocity(Body * a); //Compute velocity of body a void Body_Position(Body * a); //Compute position of body a @@ -116,7 +117,7 @@ void System_Init(System * s, const char * fileName); //Initialise System (array void System_Compute(System * s); void System_Forces(System * s1, System * s2); //Compute forces for bodies in s1 due to bodies in s2 (also updates velocities) void System_Positions(System * s); //Update positions for bodies in s1 - +System * Split_System(System * s, unsigned n); // Splits one system into a number of other systems, returns an array of size n void Universe_Cleanup(); //Cleanup universe and write bodies to file diff --git a/course/semester2/pprog/assignment1/single-thread/nbody.h~ b/course/semester2/pprog/assignment1/single-thread/nbody.h~ new file mode 100644 index 00000000..f0327f3c --- /dev/null +++ b/course/semester2/pprog/assignment1/single-thread/nbody.h~ @@ -0,0 +1,140 @@ +#ifndef _NBODY_H +#define _NBODY_H + +/** + * @file nbody.h + * @author Sam Moore (205030628) + * @purpose N-Body simulator: declarations of simulation related parameters + */ + +#include +#include +#include +#include +#include //POSIX time + + + +#define SINGLE_THREADED + +// --- The below macros will be undefined by the multithreaded versions, and replaced with functions --- // + +//Sets up the simulation; in multithreaded versions, will spawn threads +#define Simulation_Run(argc, argv) \ +if (options.draw_graphics) \ + Graphics_Run(argc, argv); \ +else \ +{ \ + while (ExitCondition() == false) \ + { BeforeDraw(); AfterDraw(); } \ +} + +#define QuitProgram(error) (runstate = (error == true) ? QUIT : QUIT_ERROR) //Prepares to exit program, is thread safe in multithreaded versions + +//Macro to be overwritten in multithreaded versions, called before the graphics is allowed to draw anything +#define BeforeDraw() \ +System_Compute(&universe); \ +universe.steps += 1; \ +if (options.verbosity != 0 && universe.steps % options.verbosity == 0) \ + DisplayStatistics(); \ + + + + +//Macro to be overwritten in multithreaded versions, called after the graphics has finished drawing. +#define AfterDraw() + +// --- Constants and other Macros --- // + +#define M_PI 3.14159265358979323846264338327950288 /* pi */ +#define G 6.67428E-11 +#define DELTA_T 0.05 +#define DIMENSIONS 3 +#define LINE_SIZE 1000 +#define square(x) ((x)*(x)) + + +/** + * Structure to represent a single Body + * @param mass - Mass of the body + * @param x - Position vector (array) + * @param v - Velocity vector (array) + * @param F - Net force vector (array) + */ +typedef struct +{ + + double mass; + double x[DIMENSIONS]; + double v[DIMENSIONS]; + double F[DIMENSIONS]; + +} Body; + +/** + * Structure to store an array of bodies, along with the size of the array. + * The universe is represented by a single System. + * In the multithreaded program, the universe is subdivided into one system for each working thread to use. + * @param N - Size of the array + * @param body - The array of bodies + */ +typedef struct +{ + unsigned N; // Number of bodies in the System + Body * body; // Array of bodies + + unsigned steps; //Number of steps simulated + +} System; + +/** + * Structure to represent options passed to the program. + */ +typedef struct +{ + const char * input; // initial body field + const char * output; // file to write final positions / velocities of bodies to + const char * program; // program name + int num_threads; // number of worker threads to spawn (must be greater than 1 for any to be spawned) + int nested_threads; // number of threads to nest computations with (must be greater than 1 for any to be spawned) + int num_steps; // number of steps to run before stopping (run indefinately if less than zero) + int timeout; // number of seconds to run before stopping (run indefinately if less than zero) + bool draw_graphics; // whether or not to actually draw graphics + bool pedantic_graphics; // whether the graphics thread will synchronise with the computation thread (true) or just draw as fast as possible (false) + bool print_positions; // print positions of bodies to stdout on every step + int verbosity; // print statistics every number of steps indicated by this variable + clock_t start_clock; // clock cycles done when simulation starts + struct timeval start_time; // time at which simulation starts +} Options; + +void Body_Print(Body * a, FILE * out); //Print body a +void Body_Force(Body * a, System * s); //Compute force on body a due to system of bodies s +void Body_Force(Body * a, Body * b); // Compute force on body a due to body b +void Body_Velocity(Body * a); //Compute velocity of body a +void Body_Position(Body * a); //Compute position of body a + +void System_Init(System * s, const char * fileName); //Initialise System (array of bodies) from a text file +void System_Compute(System * s); +void System_Forces(System * s1, System * s2); //Compute forces for bodies in s1 due to bodies in s2 (also updates velocities) +void System_Positions(System * s); //Update positions for bodies in s1 +System * Split_System(System * s, unsigned n); // Splits one system into a number of other systems, returns an array of size n + +void Universe_Cleanup(); //Cleanup universe and write bodies to file + +void DisplayStatistics(); // Print information about steps computed, total (real) runtime, and CPU cycles + + +bool ExitCondition(void); //Checks whether the program is supposed to exit automatically yet (ie: other than the user pressing "quit") + + +typedef enum {RUN, QUIT, QUIT_ERROR} RUNSTATE; +extern RUNSTATE runstate; // Set runstate to QUIT or QUIT_ERROR and the simulation will stop on the next step. +// This is fairly redundant in the single threaded version, but useful for making sure *all* threads know to exit in the multi-threaded version + + +extern System universe; // The main array of bodies; global variable. +extern Options options; // Parameters passed to program + + + +#endif //_NBODY_H -- 2.20.1