Parallel Programming - Tidied things up
authorSam Moore <sam@daedalus.(none)>
Sat, 8 Sep 2012 14:41:53 +0000 (22:41 +0800)
committerSam Moore <sam@daedalus.(none)>
Sat, 8 Sep 2012 14:41:53 +0000 (22:41 +0800)
Since most of the single thread code can be reused, I put in symbolic links for the mthread version.
I added comments. The programs can be run without graphics. They can also output performance info.

Also wrote python script to get performance information and plot it.
I was going to do the OpenMP version, but got carried away with all the "tidying", so I haven't yet.

Curses.

13 files changed:
course/semester2/pprog/assignment1/get_data.py [new file with mode: 0755]
course/semester2/pprog/assignment1/mthread/Makefile [changed from file to symlink]
course/semester2/pprog/assignment1/mthread/graphics.c [changed from file to symlink]
course/semester2/pprog/assignment1/mthread/graphics.h [changed from file to symlink]
course/semester2/pprog/assignment1/mthread/main.c [changed from file to symlink]
course/semester2/pprog/assignment1/mthread/nbody.c
course/semester2/pprog/assignment1/mthread/nbody.h
course/semester2/pprog/assignment1/single-thread/Makefile
course/semester2/pprog/assignment1/single-thread/graphics.c
course/semester2/pprog/assignment1/single-thread/graphics.h
course/semester2/pprog/assignment1/single-thread/main.c
course/semester2/pprog/assignment1/single-thread/nbody.c
course/semester2/pprog/assignment1/single-thread/nbody.h

diff --git a/course/semester2/pprog/assignment1/get_data.py b/course/semester2/pprog/assignment1/get_data.py
new file mode 100755 (executable)
index 0000000..8ef8241
--- /dev/null
@@ -0,0 +1,76 @@
+#!/usr/bin/python -u
+
+# @file get_data.py
+# @purpose Extract statistics from repeated runs of N-Body simulators
+# @author Sam Moore (20503628) - 2012
+
+import sys
+import os
+import subprocess
+
+import Gnuplot, Gnuplot.funcutils
+plot = Gnuplot.Gnuplot(persist=0)
+
+programs = {
+       "single-thread" : "./single-thread/nbody -g -v 5"
+}
+for n in range(1, 5):
+       programs.update({"mthread"+str(n) : "./mthread/nbody -g -v 5 -n "+str(n)})
+
+
+
+def RunProgram(string):
+       p = subprocess.Popen(string.split(" "), stdout=subprocess.PIPE)
+       results = []
+       for line in p.stdout:
+               results.append(map(float, line.strip(" \r\n\t").split("\t")))
+       results.sort(key = lambda e : e[0])
+       return results
+
+def RunForSteps(program, field, steps):
+       results = []
+       #print(str(program) + " -s "+str(steps) + " " + str(field))
+       results.append(RunProgram(str(program) + " -s "+str(steps) + " " + str(field)))
+       return results
+
+def VaryField(program, steps, fields):
+       results = []
+       for f in fields:
+               results.append(RunProgram(program + " -s " +str(steps) + " " + str(field)))
+       return results
+
+def VaryTime(program, field, time_min, time_max, time_inc=10):
+       results = []
+       for t in range(time_min, time_max, time_inc):
+               results.append(RunProgram(program + " -t "+str(t) + " " + str(field)))
+       return results
+
+def main(argv):
+       plot("set xlabel \"steps\"")
+       plot("set ylabel \"real time (s)\"")
+       plot("set title \"Time .vs. Steps Computed\"")
+       plot("set key outside right")
+       data = {}
+       data.update(programs)
+
+       fields = map(int, os.listdir("fields"))
+       fields.sort(key = lambda e : e)
+       
+       for n in fields:
+               for p in data.keys():
+                       data[p] = [RunForSteps(programs[p], "fields/"+str(n), 200), n]
+                       #print(str(data[p]))
+                       if (len(data[p]) > 0):
+                               plot.replot(Gnuplot.Data(data[p][0], title=str(p)+":"+str(n), with_="lp linecolor "+str(fields.index(n))))
+
+       #print(str(data.items()))
+       # Score the programs
+       
+       #for p in sorted(data.items(), key = lambda e : e[1][0][len(e[1][0])-1][1]):
+       #       print(str(p[0])+":"+str(n)+"\t"+str(p[1][0][len(p[1][0])-1][1]))
+               
+       sys.stdin.readline()
+       return 0
+
+if __name__ == "__main__":
+       sys.exit(main(sys.argv))
deleted file mode 100644 (file)
index 0f46414b020bc86fca08721f6c4f4eae18f646fc..0000000000000000000000000000000000000000
+++ /dev/null
@@ -1,26 +0,0 @@
-#Makefile for nbody program - pthread version
-
-CXX = gcc
-LIBRARIES = -lm -lGL -lglut -lGLU -lpthread
-FLAGS = --std=c99 -Wall -pedantic -g
-PREPROCESSOR_FLAGS = 
-SINGLE_THREAD_OBJ = main.o nbody.o graphics.o
-LINK_OBJ = main.o nbody.o graphics.o
-
-BIN = nbody
-
-$(BIN) : $(LINK_OBJ) 
-       $(CXX) -o $(BIN) $(LINK_OBJ) $(LIBRARIES)
-
-nbody : 
-
-%.o : %.c
-       $(CXX) $(FLAGS) $(PREPROCESSOR_FLAGS) -c $<
-
-clean :
-       $(RM) $(BIN) $(OBJ) $(LINK_OBJ)
-
-clean_full: #cleans up all backup files
-       $(RM) $(BIN) $(OBJ) $(LINK_OBJ)
-       $(RM) *.*~
-       $(RM) *~
new file mode 120000 (symlink)
index 0000000000000000000000000000000000000000..b0d614daa37de3663da7d410a0e3c913699d764b
--- /dev/null
@@ -0,0 +1 @@
+../single-thread/Makefile
\ No newline at end of file
deleted file mode 100644 (file)
index 6f043ad8ec19338d523d32519bdd675292a0ae99..0000000000000000000000000000000000000000
+++ /dev/null
@@ -1,184 +0,0 @@
-/**
- * @file graphics.c
- * @author Sam Moore (20503628) 2012 - adapted from template program provided by UWA
- * @purpose Definition of graphics related functions
- * NOTE: The graphics functions are run in the main thread.
- */
-
-#include "graphics.h" //Function declarations
-#include "nbody.h" 
-#include <stdlib.h>
-#include <stdio.h>
-#include <time.h> // Needed for clock
-#include <math.h> // Maths functions (sin, cos)
-/**
- * Variables
- */
-double previousTime, eyeTheta, eyePhi, eyeRho;
-float look[3];
-int windowWidth, windowHeight, upY;
-
-double scale = 1.0;
-
-
-/**
- * Initialise the graphics
- */
-void Graphics_Run(int argc, char ** argv) 
-{
-
-       glutInit(&argc, argv);  
-       glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
-       glutInitWindowSize(WIDTH, HEIGHT);
-       glutInitWindowPosition(POSITION_X, POSITION_Y);
-       glutCreateWindow("N-Body : Multi-Threaded (pthread)");
-       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);  
-
-
-       printf("Use:\n X - exit\n I, J, K, M - rotate\n W, Z, A, S - move to view"
-               " point\n ./, - zoom in/out\n +/- - scaled zoom in/out\n");
-
-       glutMainLoop();   
-}
-
-/**
- * This function redraws the screen after the positions of particles 
- * have been updated
- */
-void Graphics_Display() 
-{
-       //Check whether the runstate has been set to quit the program
-       switch (runstate)
-       {
-               case RUN:
-                       break;
-               case QUIT:
-                       exit(EXIT_SUCCESS);
-                       break;
-               case QUIT_ERROR:
-                       exit(EXIT_FAILURE);
-                       break;
-       }
-
-       //Check if window exists, quit if it doesn't
-       if (glutGetWindow() == 0)
-       {
-               exit(EXIT_SUCCESS);
-       }
-
-
-       glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
-       glMatrixMode(GL_MODELVIEW);
-       glLoadIdentity();
-       gluLookAt(eyeRho * sin(eyePhi) * sin(eyeTheta), eyeRho * cos(eyePhi),
-               eyeRho * sin(eyePhi) * cos(eyeTheta),look[0], look[1], look[2], 0, upY, 0);
-
-       for (int i = 0; i < universe.N; ++i) 
-       {
-               Body * b = universe.body+i;
-               glColor3f(0.0f, b->mass/1e11*100, 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();
-}
-
-/**
- * This function is to manipulate with the image
- */
-void Graphics_Keyboard(unsigned char theKey, int mouseX, int mouseY) 
-{
-       if (theKey == 'x' || theKey == 'X') 
-       {
-               exit(EXIT_SUCCESS);
-       }
-
-               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;
-}
-
-/**
- * Resize the view 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);
-}
-
-
-
new file mode 120000 (symlink)
index 0000000000000000000000000000000000000000..6fc2284c096c9d499ccfd419b3b909ba89ff91ae
--- /dev/null
@@ -0,0 +1 @@
+../single-thread/graphics.c
\ No newline at end of file
deleted file mode 100644 (file)
index 2fe01084ed24090cc35483688ceb0c763bf456d4..0000000000000000000000000000000000000000
+++ /dev/null
@@ -1,41 +0,0 @@
-#ifndef _GRAPHICS_H
-#define _GRAPHICS_H
-
-/**
- * @file graphics.h
- * @author Sam Moore (20503628) 2012 - adapted from template program provided by UWA
- * @purpose N-Body simulator - declarations of all graphics related functions
- * NOTE: I prefer to keep graphics entirely seperate from the simulation as much as possible, hence seperate files
- */
-
-#include <GL/gl.h>
-#include <GL/glut.h>
-
-#include "nbody.h"
-
-#define WIDTH 800
-#define HEIGHT 800
-#define POINT_SIZE 1
-#define POSITION_X 112
-#define POSITION_Y 20
-#define WORLD_LEFT -10000
-#define WORLD_RIGHT 10000
-#define WORLD_BOTTOM -10000
-#define WORLD_TOP 10000
-#define VIEW_ANGLE 45
-#define RHO 100
-#define WORLD_NEAR 0.1
-#define WORLD_FAR 1000000
-#define BALL_SIZE 0.5
-#define REFRESH_RATE 0.001
-#define LINE_SIZE 1000
-
-void Graphics_Run(int argc, char ** argv);
-void Graphics_Display(void);
-void Graphics_Keyboard(unsigned char key, int mouse_x, int mouse_y);
-void Graphics_Reshape(int width, int height);
-
-
-#endif //_GRAPHICS_H
-
-//EOF
new file mode 120000 (symlink)
index 0000000000000000000000000000000000000000..aebee138b9a8f60cf2a97e9f1fbe7daaf1e73c25
--- /dev/null
@@ -0,0 +1 @@
+../single-thread/graphics.h
\ No newline at end of file
deleted file mode 100644 (file)
index 531624616e1c2a271d827fc3e47fe60fedc8e742..0000000000000000000000000000000000000000
+++ /dev/null
@@ -1,59 +0,0 @@
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <pthread.h>
-
-#include "nbody.h"
-#include "graphics.h"
-
-unsigned numberOfProcessors;
-
-System universe;
-pthread_t compute_thread;
-pthread_mutex_t mutex_runstate;
-RUNSTATE runstate = RUN;
-
-void Thread_Cleanup(void)
-{
-       if (runstate == RUN)
-               QuitProgram(false);
-       pthread_join(compute_thread, NULL);
-}
-
-// This is main function. Do not change it.
-int main(int argc, char** argv)
-{
-       if (argc < 2) 
-       {
-               printf("Please provide the filename, i.e. \'%s bodiesfield.dat\'\n", argv[0]);
-               exit(EXIT_FAILURE);
-       }
-       if (argc == 2)
-       {
-               System_Init(&universe,argv[1]);
-               atexit(Universe_Cleanup);
-       }
-       if (argc == 3) 
-       {
-               numberOfProcessors = atoi(argv[2]);
-       }
-       
-       atexit(Thread_Cleanup);
-       
-       // Create a thread to handle computation loop
-       
-       if (pthread_create(&compute_thread, NULL, Compute_Thread, (void*)&universe) != 0)
-       {
-               perror("Error creating compute thread");
-               exit(EXIT_FAILURE);
-       }
-       
-       Graphics_Run(argc, argv); // Run graphics loop in the main thread
-       
-       printf("Close!\n");
-       //We get to this point when the "close" button is clicked in the graphics window
-       QuitProgram(false);
-       pthread_join(compute_thread, NULL);
-       
-}
-
new file mode 120000 (symlink)
index 0000000000000000000000000000000000000000..37e5985321bbd9f9ade691e08d60d5122cd91e5f
--- /dev/null
@@ -0,0 +1 @@
+../single-thread/main.c
\ No newline at end of file
index cd0656a..afd1d60 100644 (file)
 /**
  * @file nbody.c
- * @author Sam Moore (20503628) 2012
- * @purpose N-Body simulator - Definition of simulation functions; single threaded version
+ * @purpose Implentation of multi-threaded N-Body simulator using pthreads
+ * @author Sam Moore (20503628) - 2012
  */
 
-#include <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-#include <time.h>
-#include <string.h>
-#include <stdbool.h>
-#include <stdlib.h>
-#include <unistd.h>
-
-#include "nbody.h"
-
-//Stuff to do with worker threads
-// The worker threads (if used) are spawned by the computation thread, not the main thread
-// The main thread (which handles initialisation, graphics, user input, cleanup and exit) does not need to know about these.
-#if NUM_THREADS > 1
-       
-       pthread_t worker_thread[NUM_THREADS]; //Worker threads responsible for Force and Position updates
-       System sub_system[NUM_THREADS]; //Array of System's used to divide up the main "universe" System for worker threads
-       pthread_mutex_t mutex_workers;
-       pthread_cond_t workers_done_cv;
-       unsigned worker_count;
-#endif 
+#include "nbody.h" // Declarations
+#include "../single-thread/nbody.c" // Include all functions from the single threaded version
 
-/**
- * Prints the body on screen
- */
-void Body_Print(Body * a)
-{
-       printf("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]);
-}
+// --- Variable declarations --- //
 
-/**
- * Computing forces
- */
-void Body_Force(Body * a, System * s) 
-{
-       double distance;
-       double con;
-       double gd;
+pthread_t compute_thread; // The thread responsible for computations; it spawns worker threads
+       
+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
+pthread_mutex_t mutex_workers; // Mutex used for the barrier between Force and Position updates
+pthread_cond_t workers_done_cv; // Conditional used for the barrier between Force and Position updates
+unsigned workers_busy; // Number of workers currently doing something
 
-       for (unsigned i = 0; i < DIMENSIONS; ++i)
-               a->F[i] = 0;
+pthread_mutex_t mutex_runstate; // Mutex around the runstate
 
-       for (unsigned index = 0; index < s->N; ++index)
-       {
-               Body * b = s->body+index;
-               if (b == a)
-                       continue;
-               
-               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)
-                       a->F[i] += gd * (b->x[i] - a->x[i]);
-       }
-}
 
-/**
- * Compute velocities
- */
-void Body_Velocity(Body * a) 
-{
-       for (unsigned i = 0; i < DIMENSIONS; ++i)
-               a->v[i] += a->F[i] / a->mass * DELTA_T;
-}
 
 /**
- * Compute positions
+ * @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 (ie: &universe)
  */
-void Body_Position(Body * a) 
-{
-       for (unsigned i = 0; i < DIMENSIONS; ++i)
-               a->x[i] += a->v[i] * DELTA_T;
-}
-
-/**
- * Compute forces on each object in System s1 due to all bodies in System s2
- */
-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_Velocity(s1->body+i);
-       }
-       //printf("Compute positions for %p - Done!\n", (void*)s1);
-}
-
-/**
- * Compute positions of all objects in System
- */
-void System_Positions(System * s)
+void * Compute_Thread(void * arg)
 {
-       //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);
-}
-
-
 
+       System * s = (System*)(arg); //cast argument to a System*
 
-/*
- * This function reads an input file. You can change it if you choose a 
- * different file format
- */
-#define LINE_SIZE BUFSIZ
-void System_Init(System * s, 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));
 
-       //printf("----------------------Initial field-------------------------------\n");
+       // If no number of threads provided, use the default value, unless someone changed that to a stupid value
+       if (options.num_threads <= 0)
+               options.num_threads = (DEFAULT_WORKING_THREADS > 1) ? DEFAULT_WORKING_THREADS : 1;
 
-       for (unsigned i = 0; i < s->N; ++i)
+       // Do a sanity check; there is no point spawning more threads than bodies.
+       if (options.num_threads > s->N)
        {
-               if (fgets(line, LINE_SIZE, file) != NULL)
-               {
-                       Body * a = s->body+i;
-                       token = strtok(line, ",; ");
-                       a->mass = atof(token);
-                       
-                       for (unsigned j = 0; j < DIMENSIONS; ++j)
-                       {
-                               token = strtok(NULL, ",; ");
-                               a->x[j] = atof(token);
-                       }
-                       for (unsigned j = 0; j < DIMENSIONS; ++j)
-                       {
-                               token = strtok(NULL, ",; ");
-                               a->v[j] = atof(token);
-                       }
-                       //Body_Print(a);
-               }
+               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;
        }
-
-       //printf("--------------------------------------------------------------\n");
        
-       fclose(file);
-}
-
-/**
- * Cleans up the universe by freeing all memory
- */
-void Universe_Cleanup()
-{
-       free(universe.body);
-       
-}
-
-/**
- * Thread - Continuously computes steps for a system of bodies
- */
-void * Compute_Thread(void * arg)
-{
+       if (options.num_threads > 1) // Allocate worker threads and sub systems, as long as there would be more than 1
+       {
+               worker_thread = (pthread_t*)(calloc(options.num_threads, sizeof(pthread_t)));
+               if (worker_thread == NULL)
+               {
+                       perror("Couldn't allocate array of worker threads");
+                       QuitProgram(true);
+                       pthread_exit(NULL);
+               }
+               sub_system = (System*)(calloc(options.num_threads, sizeof(System)));
+               if (sub_system == NULL)
+               {
+                       perror("Couldn't allocate array of systems for worker threads to use");
+                       QuitProgram(true);
+                       pthread_exit(NULL);
+               }
 
-       System * s = (System*)(arg); //cast argument to a System*
-       // First set up the "sub_system" array
-       #if NUM_THREADS > 1
-               unsigned bodies_per_system = (s->N) / NUM_THREADS;
-               unsigned remainder = (s->N) % NUM_THREADS;
-               for (unsigned i = 0; i < NUM_THREADS; ++i)
+               // Divide up the Body array owned by s into options.num_threads arrays, one for each worker thread
+               unsigned bodies_per_system = (s->N) / options.num_threads;
+               unsigned remainder = (s->N) % options.num_threads;
+               for (unsigned i = 0; i < options.num_threads; ++i)
                {
                        sub_system[i].body = (s->body)+(i*bodies_per_system);
                        sub_system[i].N = bodies_per_system;
-                       if (i == NUM_THREADS - 1)
-                               sub_system[i].N += remainder;
+                       sub_system[i].steps = 0;
+                       if (i == options.num_threads - 1)
+                               sub_system[i].N += remainder; // The last thread gets the remainder
+
                }
+       }
+
 
-               pthread_attr_t attr; //thread attribute for the workers
-               pthread_attr_init(&attr);
-               pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_DETACHED);
+       pthread_attr_t attr; //thread attribute for the workers. 
+       pthread_attr_init(&attr);
+       pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_DETACHED); //Needs to be detached, so that memory can be reused.
 
-       #endif
+       // The main computation loop
        while (true)
        {
                
                if (runstate != RUN) pthread_exit(NULL); //Check whether the main thread wants to exit
-       
+
 
        
-               #if NUM_THREADS <= 1
-                       //No worker threads; do everything in this thread
+               //Check whether the program should quit due to steps being computed, or a timeout
+               if (options.timeout > 0.0)
+               {
+                       if ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)
+                               QuitProgram(false);
+               }
+
+               if (options.num_steps > 0 && universe.steps > options.num_steps)
+                       QuitProgram(false);
+
+               if (options.draw_graphics == false && options.verbosity != 0 
+                       && universe.steps % options.verbosity == 1)
+               {
+                       DisplayStatistics();
+               }
+
+
+               if (options.num_threads <= 1)
+               {
+                       // Just do everything in this thread 
                        System_Forces(s, s);
                        System_Positions(s);
-               #else
+                       continue;
+               }
 
                
 
                
-               worker_count = NUM_THREADS; //All threads working
+               workers_busy = options.num_threads; //All threads working
 
                //Compute forces
-               for (unsigned i = 0; i < NUM_THREADS; ++i)
+               for (unsigned i = 0; i < options.num_threads; ++i)
                {
                        if (pthread_create(worker_thread+i, &attr, Force_Thread, (void*)(sub_system+i)) != 0)
                        {
@@ -219,16 +133,18 @@ void * Compute_Thread(void * arg)
 
 
 
-               //Wait for forces to be computed
+               //Barrier - Wait for forces to be computed
                pthread_mutex_lock(&mutex_workers);
-               while (worker_count > 0)
+               while (workers_busy > 0)
                        pthread_cond_wait(&workers_done_cv, &mutex_workers);
                pthread_mutex_unlock(&mutex_workers);
+
+               //All the forces are now computed
                
-               worker_count = NUM_THREADS; //All threads working
+               workers_busy = options.num_threads; //All threads working
 
                //Compute positions
-               for (unsigned i = 0; i < NUM_THREADS; ++i)
+               for (unsigned i = 0; i < options.num_threads; ++i)
                {
                        if (pthread_create(worker_thread+i, &attr, Position_Thread, (void*)(sub_system+i)) != 0)
                        {
@@ -240,58 +156,120 @@ void * Compute_Thread(void * arg)
 
                //Wait for positions to be computed
                pthread_mutex_lock(&mutex_workers);
-               while (worker_count > 0)
+               while (workers_busy > 0)
                        pthread_cond_wait(&workers_done_cv, &mutex_workers);
-               pthread_mutex_unlock(&mutex_workers);           
+               pthread_mutex_unlock(&mutex_workers);
+
+               //Update number of steps computed
+               universe.steps += 1;
+
 
-               
-               #endif
 
        }
 }
 
+/**
+ * @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_Forces((System*)s, &universe);
+       System_Forces((System*)s, &universe); //Simple wrapper
+
        pthread_mutex_lock(&mutex_workers);
-       worker_count -= 1;
-       if (worker_count == 0)
+       workers_busy -= 1;      // Count this thread as done
+       if (workers_busy == 0)
        {
-               pthread_cond_signal(&workers_done_cv);
+               pthread_cond_signal(&workers_done_cv); // All threads done; wake up the compute_thread
        }
        pthread_mutex_unlock(&mutex_workers);
        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);
+       System_Positions((System*)s); // Simple wrapper
+
        pthread_mutex_lock(&mutex_workers);
-       worker_count -= 1;
-       if (worker_count == 0)
+       workers_busy -= 1; // Count this thread as done
+       if (workers_busy == 0)
        {
-               pthread_cond_signal(&workers_done_cv);
+               pthread_cond_signal(&workers_done_cv); //All threads done; wake up the compute_thread
        }
        pthread_mutex_unlock(&mutex_workers);
        return NULL;
 }      
 
 /**
- * Child threads can call this to signal the main thread to quit the program
- * The main thread also uses this to tell child threads that the program is quitting
- *     Note that the function doesn't call exit(), that is done by the main thread
+ * @function QuitProgram
+ * @purpose This function can either be called by the main thread in order to signal other threads
+ *             that it wants to exit. The main thread then calls pthread_join before exiting.
+ *     It can also be called by a child thread to request the main thread to exit.
+ *     It is only used this way if there is an unrecovarable error (ie: Can't allocate memory in a child thread)
  */
 void QuitProgram(bool error)
 {
        
-       pthread_mutex_lock(&mutex_runstate);
-       if (error)
+       pthread_mutex_lock(&mutex_runstate); // aquire mutex
+       if (error) // set the runstate
                runstate = QUIT_ERROR;
        else
                runstate = QUIT;
-       pthread_mutex_unlock(&mutex_runstate);
+       pthread_mutex_unlock(&mutex_runstate); //release mutex
+}
+
+/**
+ * @function Thread_Cleanup
+ * @purpose Will be called in the main thread when exit() is called
+ *     Automatically tells all other threads to quit (if they haven't already been told) 
+ *     Then waits for them to finish.
+ *     Also frees memory associated with the worker threads.   
+ */
+void Thread_Cleanup(void)
+{
+       if (runstate == RUN) // If this is true, as far as child threads are concerned, the simulation is still running
+               QuitProgram(false); // So call QuitProgram which will set runstate, and cause child threads to exit
+       pthread_join(compute_thread, NULL);
+       free(worker_thread);
+       free(sub_system);
 }
 
 
+/**
+ * @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 compute thread
+ */
+void Simulation_Run()
+{
+       atexit(Thread_Cleanup);
+
+       if (options.draw_graphics == false)
+       {
+               Compute_Thread((void*)(&universe));
+
+               // If there are no graphics, run the computation loop in the main thread
+               // The call to pthread_exit(NULL) in the computation loop ends the main thread
+
+               // This means the main function never calls the Graphics_Run function
+               // Without this condition here, Graphics_Display would essentially be running a busy loop in the main thread
+               
+               fprintf(stderr,"Should not see this\n");
+               exit(EXIT_FAILURE);
+       }
+       
+       // Create a thread to handle computation loop
+       if (pthread_create(&compute_thread, NULL, Compute_Thread, (void*)&universe) != 0)
+       {
+               perror("Error creating compute thread");
+               exit(EXIT_FAILURE);
+       }
+}
index b9cac6f..9b76bc7 100644 (file)
@@ -1,87 +1,38 @@
-#ifndef _NBODY_H
-#define _NBODY_H
+#ifndef _NBODY_MTHREAD_H
+#define _NBODY_MTHREAD_H
 
-/**
- * @file nbody.h
- * @author Sam Moore (205030628)
- * @purpose N-Body simulator: declarations of simulation related parameters
- */
-
-
-#define NUM_THREADS 5 //The number of *worker* threads (not including the computation thread) to use. 
-                       // If this is 1, the computation thread will not spawn any workers (faster to just do the work itself), 
-                       //      otherwise it will spawn NUM_THREADS workers.
 
+#include "../single-thread/nbody.h" //Use original simulation code
 #include <pthread.h>
-#include <stdbool.h>
-
-#define M_PI        3.14159265358979323846264338327950288   /* pi */
-#define G 6.67428E-11
-#define DELTA_T 0.05
-#define DIMENSIONS 3
-#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 in a single System. 
- * @param N - Size of the array
- * @param body - The array of bodies
- */
-typedef struct
-{
-       unsigned N;
-       Body * body;
 
-} System;
+#undef SINGLE_THREADED
+#define PTHREADED
 
-void Body_Print(Body * a); //Print body a
-void Body_Force(Body * a, System * s); //Compute force on body a due to system of bodies s
-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, char * fileName); //Initialise System (array of bodies) from a text file
-
-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
-void Universe_Cleanup();
-
-extern System universe; // The main array of bodies; global variable.
+#define DEFAULT_WORKING_THREADS 2 
 
+//Undefine default macros, replace with functions
+#undef Simulation_Run
+void Simulation_Run(void);
+#undef QuitProgram
+void QuitProgram(bool error);
 
-/**
- * Multithreading stuff below here
- */
 
-extern pthread_t compute_thread; // ID of the thread that runs Compute_Thread. Set in main()
 void * Compute_Thread(void * system); //Thread - Continuously perform computations for a System of bodies. May spawn additional worker 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
 
-typedef enum {RUN, QUIT, QUIT_ERROR} RUNSTATE;
-extern RUNSTATE runstate;
-extern pthread_mutex_t mutex_runstate; //Mutex around the "runstate" variable
 
-void QuitProgram(bool error);
 
+void Thread_Cleanup(void); //Called at program exit to safely join computation thread
+
+extern pthread_t compute_thread; // ID of the thread that runs Compute_Thread. 
+extern pthread_mutex_t mutex_runstate; //Mutex around the "runstate" variable
 
+extern pthread_t * worker_thread; //Array of worker threads responsible for Force and Position updates
+extern System * sub_system; //Array of Systems used to divide up the main "universe" System for worker threads
+extern pthread_mutex_t mutex_workers;
+extern pthread_cond_t workers_done_cv;
+extern unsigned workers_busy;
 
-#endif //_NBODY_H
+#endif //_NBODY_MTHREAD_H
index fbece84..f34eab2 100644 (file)
@@ -1,9 +1,10 @@
-#Makefile for nbody program - single threaded version
-
+# Makefile for nbody program
+# Compiles on Ubuntu 12.04 and Debian 6.0.4
+# NOTE: This file is identical for both the single-threaded and multi-threaded versions of the program
 CXX = gcc
-LIBRARIES = -lm -lGL -lglut -lGLU
+LIBRARIES = -lm -lGL -lglut -lGLU -lpthread
 FLAGS = --std=c99 -Wall -pedantic -g
-PREPROCESSOR_FLAGS = 
+PREPROCESSOR_FLAGS = -fopenmp
 SINGLE_THREAD_OBJ = main.o nbody.o graphics.o
 LINK_OBJ = main.o nbody.o graphics.o
 
index d3a8b3c..e3722e9 100644 (file)
@@ -1,28 +1,64 @@
 /**
  * @file graphics.c
- * @author Sam Moore (20503628) 2012 - adapted from template program provided by UWA
+ * @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 "graphics.h" //Function declarations
+
+#include <stdlib.h>
+#include <stdio.h>
 #include <time.h> // Needed for clock
 #include <math.h> // Maths functions (sin, cos)
-/**
- * Variables
- */
+
+#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;
 
 
 /**
- * Initialise the graphics
+ * @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_Init(void
+void Graphics_Run(int argc, char ** argv
 {
-    
+       if (options.draw_graphics == false)
+       {
+               while (true)
+                       Graphics_Display();
+                       
+               return;
+       }       
+
+       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);
@@ -58,14 +94,56 @@ void Graphics_Init(void)
        look[1] = 0;
        look[2] = 0;
        gluPerspective(VIEW_ANGLE, (double)WIDTH/(double)HEIGHT, WORLD_NEAR, WORLD_FAR);  
+
+       glutMainLoop();   
 }
 
 /**
- * This function redraws the screen after the positions of particles 
- * have been updated
+ * @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.verbosity != 0 && universe.steps % options.verbosity == 1)
+               DisplayStatistics();
+
+       #ifdef SINGLE_THREADED
+               System_Compute(&universe);
+
+               //Check whether the program should quit due to steps being computed, or a timeout
+               if (options.timeout > 0.0)
+               {
+                       if ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)
+                               exit(EXIT_SUCCESS);
+               }
+               if (options.num_steps > 0 && universe.steps > options.num_steps)
+                       exit(EXIT_SUCCESS);
+
+
+       #endif
+
+
+       //Check whether the runstate has been set to quit the program
+       switch (runstate)
+       {
+               case RUN:
+                       break;
+               case QUIT:
+                       exit(EXIT_SUCCESS);
+                       break;
+               case QUIT_ERROR:
+                       exit(EXIT_FAILURE);
+                       break;
+       }
+
+       if (options.draw_graphics == false)
+               return;
+       
+
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
        glMatrixMode(GL_MODELVIEW);
        glLoadIdentity();
@@ -85,7 +163,10 @@ void Graphics_Display()
 }
 
 /**
- * This function is to manipulate with the image
+ * @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) 
 {
@@ -126,7 +207,9 @@ void Graphics_Keyboard(unsigned char theKey, int mouseX, int mouseY)
 }
 
 /**
- * Resize the view window
+ * @function Graphics_Reshape
+ * @purpose Resize the view window
+ * @param width, height - New size of the window
  */
 void Graphics_Reshape(int width, int height) 
 {
@@ -140,12 +223,4 @@ void Graphics_Reshape(int width, int height)
 }
 
 
-/**
- * This function is called repeatedly by graphics library. You can consider 
- * it as main loop in the program.
- */
-void Graphics_Animate(void) 
-{
-       System_Compute(&universe); //Compute and update new positions for the time step
-       Graphics_Display(); // Display needs to be called to redraw the screen
-}
+
index b3f9dc5..2fbc0a6 100644 (file)
@@ -5,7 +5,7 @@
  * @file graphics.h
  * @author Sam Moore (20503628) 2012 - adapted from template program provided by UWA
  * @purpose N-Body simulator - declarations of all graphics related functions
- * NOTE: I prefer to keep graphics entirely seperate from the simulation as much as possible, hence seperate files
+ * NOTE: I prefer to keep graphics seperate from the simulation as much as possible, hence seperate files
  */
 
 #include <GL/gl.h>
 #define REFRESH_RATE 0.001
 #define LINE_SIZE 1000
 
-void Graphics_Init(void);
+void Graphics_Run(int argc, char ** argv);
 void Graphics_Display(void);
 void Graphics_Keyboard(unsigned char key, int mouse_x, int mouse_y);
 void Graphics_Reshape(int width, int height);
-void Graphics_Animate(void);
+
 
 #endif //_GRAPHICS_H
 
index 63f0d2e..f9f3b52 100644 (file)
+/**
+ * @file main.c
+ * @author Sam Moore (20503628) - 2012
+ * @purpose Contains main function, argument handling
+ * NOTE: This file is identical for both the single-threaded and multi-threaded versions of the program
+ */
 
 #include <stdlib.h>
 #include <stdio.h>
+#include <time.h>
+#include <signal.h>
 
 #include "nbody.h"
 #include "graphics.h"
 
-unsigned numberOfProcessors;
+// --- Variable definitions --- //
+System universe; // global variable declared in "nbody.h" - The universe of bodies
+Options options; // global variable declared in "nbody.h" - Options passed to program through arguments
 
-System universe;
 
-// This is main function. Do not change it.
+// --- Function forward declarations --- //
+void HandleArguments(int argc, char ** argv); //Interprets program arguments and sets up the "options" variable
+void UnsignedArgument(unsigned * i, int argc, char ** argv, unsigned * store); //Helper function to get switch value for unsigned
+void FloatArgument(unsigned * i, int argc, char ** argv, float * store); //Helper function to get switch value for float
+void DisplayStatistics(); //Called on exit of program, displays information about computation time, steps computed, etc
+void Interrupt(int dummy); // Interrupt handler function, called when SIGINT (Ctrl-C) is sent to program
+
+// --- Function implementations --- //
+
+/**
+ * @function main
+ * @purpose The main function. Calls HandleArguments to set up options, 
+ *             sets up functions to call at exit, starts simulation and graphics loops
+ * @param argc - Number of arguments
+ * @param argv - The argument strings
+ */
 int main(int argc, char** argv)
 {
-       glutInit(&argc, argv);
+       //Set default options values here
+       options.program = argv[0];
+       options.input = NULL;
+       options.output = NULL;
+       options.num_threads = 0;
+       options.num_steps = 0;
+       options.timeout = 0;
+       options.draw_graphics = true;
+       options.print_positions = false;
+       options.verbosity = 0;
 
-       if (argc < 2) 
+       //If there are no arguments, print information about usage
+       if (argc == 1)
        {
-               printf("Please provide the filename, i.e. \'%s bodiesfield.dat\'\n", argv[0]);
+               printf("Usage: %s [OPTIONS] field\n", argv[0]);
+               printf("Controls:\n X - exit\n I, J, K, M - rotate\n W, Z, A, S - move to view"
+                       " point\n ./, - zoom in/out\n +/- - scaled zoom in/out\n");
                exit(EXIT_SUCCESS);
        }
-       if (argc == 2)
+
+       //Parse the arguments
+       HandleArguments(argc, argv);
+       
+       //Check there is an initial field.
+       if (options.input == NULL)
        {
-               System_Init(&universe,argv[1]);
-               atexit(Universe_Cleanup);
+               fprintf(stderr, "Usage: %s [OPTIONS] field\n", argv[0]);
+               fprintf(stderr, " (You did not provide a file for the initial field of bodies)\n");
+               exit(EXIT_FAILURE);
        }
-       if (argc == 3) 
+
+       signal(SIGINT, Interrupt); //Handle SIGINT signals
+       atexit(Universe_Cleanup); //On exit, cleanup universe (and write positions of bodies to file if supplied).
+       atexit(DisplayStatistics); //On exit, print information about the computations done
+       System_Init(&universe,options.input); //Initialise the universe from the initial field
+       
+       //Setup the time of day at which the simulation starts
+       if (gettimeofday(&(options.start_time), NULL) != 0)
        {
-               numberOfProcessors = atoi(argv[2]);
+               perror("Couldn't get time of day");
+               exit(EXIT_FAILURE);
        }
+       
+       options.start_clock = clock(); //Get CPU cycles executed before simulation starts
+       Simulation_Run(); // Start the simulation
+       Graphics_Run(argc, argv); //Start the graphics loop
 
+       exit(EXIT_FAILURE); //Should never get to this line
+}
+
+
+/**
+ * @function HandleArguments
+ * @purpose Fill the "options" variable by parsing command line arguments
+ * @param argc - Number of arguments
+ * @param argv - Argument strings
+ */
+void HandleArguments(int argc, char ** argv)
+{
+       for (unsigned i = 1; i < argc; ++i) //For all arguments
+       {
+               if (argv[i][0] != '-') //If the argument is not a switch value...
+               {
+                       //Interpret first non switch as input file, and second as output file.
+                       //If there are too many non switch arguments, give an error.
+                       if (options.input == NULL) 
+                               options.input = argv[i];
+                       else if (options.output == NULL)
+                               options.output = argv[i];
+                       else
+                       {
+                               fprintf(stderr,"Usage: %s [OPTIONS] field [output]\n", argv[0]);
+                               fprintf(stderr," (You provided too many file names)\n");
+                               exit(EXIT_FAILURE);
+                       }
+                       continue;
+               }
+               switch (argv[i][1]) //The argument is a switch if we get here
+               {
+                       case 'n': //Number of threads switch
+                               UnsignedArgument(&i, argc, argv, &(options.num_threads));
+                               #ifdef SINGLE_THREADED
+                               fprintf(stderr, "Warning: -%c switch has no effect in the single-threaded program.\n", argv[i][1]);
+                               options.num_threads = 1;
+                               #endif //SINGLE_THREADED
+                               break;
+                       case 's': //Number of steps switch
+                               UnsignedArgument(&i, argc, argv, &(options.num_steps));
+                               break;
+                       case 't': //Timeout switch (in seconds)
+                               UnsignedArgument(&i, argc, argv, &(options.timeout));
+                               break;
+                       case 'g': //Graphics switch
+                               options.draw_graphics = !options.draw_graphics;
+                               break;
+                       case 'v': //Verbosity switch
+                               UnsignedArgument(&i, argc, argv, &(options.verbosity));
+                               break;
+
+                       default:
+                               fprintf(stderr, "Unrecognised switch -%c\n", argv[i][1]);
+                               exit(EXIT_FAILURE);
+                               break;
+               }
+               
+       }
+}
 
+/**
+ * @function UnsignedArgument
+ * @purpose Helper function to get an unsigned integer following a argument switch
+ * @param i - position in the argument array, will be updated after this function
+ * @param argc - number of arguments
+ * @param argv - argument strings
+ * @param store - pointer to unsigned to store result in
+ */
+void UnsignedArgument(unsigned * i, int argc, char ** argv, unsigned * store)
+{
+       if (*i >= argc-1)
+       {
+               fprintf(stderr,"Supply a positive integer for the -%c switch.\n", argv[*i][1]);
+               exit(EXIT_FAILURE);
+       }
+       int val = atoi(argv[*i+1]);
+       if (val <= 0)
+       {
+               fprintf(stderr,"Supply a positive integer for the -%c switch. %s is invalid.\n", argv[*i][1], argv[*i+1]);
+               exit(EXIT_FAILURE);
+       }
+       *store = val;
+       *i += 1;
+}
+/**
+ * @function FloatArgument
+ * @purpose Helper function to get a float following a argument switch
+ * @param i - position in the argument array, will be updated after this function
+ * @param argc - number of arguments
+ * @param argv - argument strings
+ * @param store - pointer to float to store result in
+ */
+void FloatArgument(unsigned * i, int argc, char ** argv, float * store)
+{
+       if (*i >= argc-1)
+       {
+               fprintf(stderr,"Supply a float for the -%c switch.\n", argv[*i][1]);
+               exit(EXIT_FAILURE);
+       }
+       *store = atof(argv[*i+1]);
+       if (*store == 0.0 && argv[*i+1][0] != '0')
+       {
+               fprintf(stderr,"Supply a float for the -%c switch.\n", argv[*i][1]);
+               exit(EXIT_FAILURE);
+       }
+       *i += 1;
+}
 
-       glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
-       glutInitWindowSize(WIDTH, HEIGHT);
-       glutInitWindowPosition(POSITION_X, POSITION_Y);
-       glutCreateWindow("N-Body : Single Threaded");
-       glutDisplayFunc(Graphics_Display);
-       glutIdleFunc(Graphics_Animate);
-       glutKeyboardFunc(Graphics_Keyboard);
-       glutReshapeFunc(Graphics_Reshape);
-       Graphics_Init();
 
 
-       printf("Use:\n X - exit\n I, J, K, M - rotate\n W, Z, A, S - move to view"
-               " point\n ./, - zoom in/out\n +/- - scaled zoom in/out\n");
+/**
+ * @function Interrupt
+ * @purpose Handle SIGINT signal; quit the program gracefully, so that statistics are still displayed
+ * @param dummy - a dummy
+ */
 
-       glutMainLoop();
-       
+void Interrupt(int dummy)
+{
+       QuitProgram(false);
 }
+
index 1dedde7..47820d1 100644 (file)
@@ -9,20 +9,31 @@
 #include <math.h>
 #include <time.h>
 #include <string.h>
+#include <stdbool.h>
+#include <stdlib.h>
+#include <unistd.h>
 
 #include "nbody.h"
 
+RUNSTATE runstate = RUN;
+
 /**
- * Prints the body on screen
+ * @function Body_Print
+ * @purpose Print the body to a file (or stdout)
+ * @param a - Body to print'
+ * @param out - FILE* to print to
  */
-void Body_Print(Body * a)
+void Body_Print(Body * a, FILE * out)
 {
-       printf("Body %p M=%f X=%f Y=%f Z=%f Fx=%f Fy=%f Fz=%f Vx=%f Vy=%f Vz=%f\n", 
+       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]);
 }
 
 /**
- * Computing forces
+ * @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
  */
 void Body_Force(Body * a, System * s) 
 {
@@ -30,28 +41,31 @@ void Body_Force(Body * a, System * s)
        double con;
        double gd;
 
-       for (unsigned i = 0; i < DIMENSIONS; ++i)
+       for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero
                a->F[i] = 0;
 
-       for (unsigned index = 0; index < s->N; ++index)
+       for (unsigned index = 0; index < s->N; ++index) // Step through all bodies in the System
        {
-               Body * b = s->body+index;
+               Body * b = s->body+index; //Current body
                if (b == a)
-                       continue;
+                       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)
+               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]);
        }
 }
 
 /**
- * Compute velocities
+ * @function Body_Velocity
+ * @purpose Compute velocity of a body
+ * @param a - The Body
  */
 void Body_Velocity(Body * a) 
 {
@@ -60,7 +74,9 @@ void Body_Velocity(Body * a)
 }
 
 /**
- * Compute positions
+ * @function Body_Position
+ * @purpose Compute position of a body
+ * @param a - The Body
  */
 void Body_Position(Body * a) 
 {
@@ -69,32 +85,59 @@ void Body_Position(Body * a)
 }
 
 /**
- * Main compute function
+ * @function System_Compute
+ * @purpose Compute forces and then positions for bodies in System
+ * NOTE: Only used in the single threaded version of the program
  */
-void System_Compute(System * s) 
+void System_Compute(System * s)
 {
-//     clock_t start, finish;
-
-//     start = clock();
+       System_Forces(s, s);
+       System_Positions(s);
+}
 
-       for (unsigned i = 0; i < s->N; ++i)
+/**
+ * @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
+ */
+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(s->body+i, s);
-               Body_Velocity(s->body+i);
-               Body_Position(s->body+i);
+               Body_Force(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
+ */
+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);
+       s->steps += 1;
 }
 
 
 
 
-/*
- * This function reads an input file. You can change it if you choose a 
- * different file format
+/**
+ * @function System_Init
+ * @purpose Initialise a System from an input file
+ * @param s - The System
+ * @param fileName - The input file
  */
 #define LINE_SIZE BUFSIZ
-void System_Init(System * s, char *fileName) 
+void System_Init(System * s, const char *fileName) 
 {
        char line[LINE_SIZE];
        char * token;
@@ -103,6 +146,7 @@ void System_Init(System * s, char *fileName)
        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");
 
@@ -134,10 +178,51 @@ void System_Init(System * s, char *fileName)
 }
 
 /**
- * Cleans up the universe by freeing all memory
+ * @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.
  */
 void Universe_Cleanup()
 {
-       free(universe.body);
+       if (options.output != NULL) // An output file was specified...
+       {
+               FILE * save = fopen(options.output, "w");
+               if (save == NULL)
+                       perror("Couldn't save universe to file.");
+               else
+               {
+                       fprintf(save, "# Final field:\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
+ */
+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);
 }
 
index e288c6d..b8c9813 100644 (file)
@@ -7,8 +7,20 @@
  * @purpose N-Body simulator: declarations of simulation related parameters
  */
 
+#include <stdlib.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <time.h>
+#include <sys/time.h> //POSIX time
+
+          
+
 #define SINGLE_THREADED
 
+//These macros will be undefined by the multithreaded versions
+#define Simulation_Run() //Sets up the simulation; in multithreaded versions, will spawn threads
+#define QuitProgram(error) (runstate = (error == true) ? QUIT : QUIT_ERROR) //Prepares to exit program, is thread safe in multithreaded versions
+
 #define M_PI        3.14159265358979323846264338327950288   /* pi */
 #define G 6.67428E-11
 #define DELTA_T 0.05
@@ -35,29 +47,63 @@ typedef struct
 
 /**
  * Structure to store an array of bodies, along with the size of the array.
- * The universe is represented in a single System. 
+ * 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;
-       Body * body;
+       unsigned N; // Number of bodies in the System
+       Body * body; // Array of bodies
+
+       unsigned steps; //Number of steps simulated
 
 } System;
 
-void Body_Print(Body * a); //Print body a
+/**
+ * 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
+       unsigned num_threads; // number of worker threads to spawn (must be greater than 1 for any to be spawned)
+       unsigned num_steps; // number of steps to run before stopping (run indefinately if equal to zero)
+       unsigned timeout; // number of seconds to run before stopping (run indefinately if equal to zero)
+       bool draw_graphics; // whether or not to actually draw graphics
+       bool print_positions; // print positions of bodies to stdout on every step
+       unsigned 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_Velocity(Body * a); //Compute velocity of body a
 void Body_Position(Body * a); //Compute position of body a
 
-void System_Init(System * s, char * fileName); //Initialise System (array of bodies) from a text file
+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
+
+
+void Universe_Cleanup(); //Cleanup universe and write bodies to file
+
+void DisplayStatistics(); // Print information about steps computed, total (real) runtime, and CPU cycles
 
-void System_Compute(System * system); //Perform a single computation step for a System of bodies
 
-void Universe_Cleanup();
+
+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

UCC git Repository :: git.ucc.asn.au