major refactoring of the code, split each distinct type into its own source file.

Former-commit-id: d291ae4cadb5810b8982c69933dfb08bf0588686

src/Makefile.am

@@ -26,8 +26,8 @@ AM_LDFLAGS = $(LAPACK_LIBS) $(BLAS_LIBS) -version-info 0:0:0
 
 # Build the libgadgetsmp library
 lib_LTLIBRARIES = libgadgetsmp.la
-libgadgetsmp_la_SOURCES = space.c runner.c
+libgadgetsmp_la_SOURCES = space.c runner.c queue.c task.c cell.c
 
 # List required headers
-include_HEADERS = space.h runner.h lock.h
+include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h gadgetsmp.h
 

src/cell.c

+/*******************************************************************************
+ * This file is part of GadgetSMP.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * 
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <pthread.h>
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+
+/* Local headers. */
+#include "cycle.h"
+#include "lock.h"
+#include "task.h"
+#include "part.h"
+#include "cell.h"
+
+
+/* Error macro. */
+#define error(s) { printf( "%s:%s:%i: %s\n" , __FILE__ , __FUNCTION__ , __LINE__ , s ); abort(); }
+
+
+/* The timers. */
+ticks cell_timer[ cell_timer_count ];
+
+
+/* Define the timer macros. */
+#ifdef TIMER_VERBOSE
+    #define TIMER
+#endif
+#ifdef TIMER
+    #define TIMER_TIC ticks tic = getticks();
+    #define TIMER_TOC(t) timer_toc( t , tic )
+    #define TIMER_TIC2 ticks tic2 = getticks();
+    #define TIMER_TOC2(t) timer_toc( t , tic2 )
+    #ifndef INLINE
+    # if __GNUC__ && !__GNUC_STDC_INLINE__
+    #  define INLINE extern inline
+    # else
+    #  define INLINE inline
+    # endif
+    #endif
+    INLINE ticks timer_toc ( int t , ticks tic ) {
+        ticks d = (getticks() - tic);
+        __sync_add_and_fetch( &cell_timer[t] , d );
+        return d;
+        }
+#else
+    #define TIMER_TIC
+    #define TIMER_TOC(t)
+#endif
+
+/**
+ * @brief Lock a cell and hold its parents.
+ *
+ * @param c The #cell.
+ */
+ 
+int cell_locktree( struct cell *c ) {
+
+    struct cell *finger, *finger2;
+    TIMER_TIC
+
+    /* First of all, try to lock this cell. */
+    if ( lock_trylock( &c->lock ) != 0 ) {
+        TIMER_TOC(cell_timer_tree);
+        return 1;
+        }
+        
+    /* Did somebody hold this cell in the meantime? */
+    if ( c->hold ) {
+        
+        /* Unlock this cell. */
+        if ( lock_unlock( &c->lock ) != 0 )
+            error( "Failed to unlock cell." );
+            
+        /* Admit defeat. */
+        TIMER_TOC(cell_timer_tree);
+        return 1;
+    
+        }
+        
+    /* Climb up the tree and lock/hold/unlock. */
+    for ( finger = c->parent ; finger != NULL ; finger = finger->parent ) {
+    
+        /* Lock this cell. */
+        if ( lock_trylock( &finger->lock ) != 0 )
+            break;
+            
+        /* Increment the hold. */
+        __sync_fetch_and_add( &finger->hold , 1 );
+        
+        /* Unlock the cell. */
+        if ( lock_unlock( &finger->lock ) != 0 )
+            error( "Failed to unlock cell." );
+    
+        }
+        
+    /* If we reached the top of the tree, we're done. */
+    if ( finger == NULL ) {
+        TIMER_TOC(cell_timer_tree);
+        return 0;
+        }
+        
+    /* Otherwise, we hit a snag. */
+    else {
+    
+        /* Undo the holds up to finger. */
+        for ( finger2 = c->parent ; finger2 != finger ; finger2 = finger2->parent )
+            __sync_fetch_and_sub( &finger2->hold , 1 );
+            
+        /* Unlock this cell. */
+        if ( lock_unlock( &c->lock ) != 0 )
+            error( "Failed to unlock cell." );
+            
+        /* Admit defeat. */
+        TIMER_TOC(cell_timer_tree);
+        return 1;
+    
+        }
+
+    }
+    
+    
+/**
+ * @brief Unock a cell's parents.
+ *
+ * @param c The #cell.
+ */
+ 
+void cell_unlocktree( struct cell *c ) {
+
+    struct cell *finger;
+    TIMER_TIC
+
+    /* First of all, try to unlock this cell. */
+    if ( lock_unlock( &c->lock ) != 0 )
+        error( "Failed to unlock cell." );
+        
+    /* Climb up the tree and unhold the parents. */
+    for ( finger = c->parent ; finger != NULL ; finger = finger->parent )
+        __sync_fetch_and_sub( &finger->hold , 1 );
+        
+    TIMER_TOC(cell_timer_tree);
+        
+    }
+    
+    
+/**
+ * @brief Sort the parts into eight bins along the given pivots.
+ *
+ * @param c The #cell array to be sorted.
+ */
+ 
+void cell_split ( struct cell *c  ) {
+
+    int i, j, k, kk;
+    struct part temp, *parts = c->parts;
+    int left[8], right[8];
+    double pivot[3];
+    
+    /* Init the pivot. */
+    for ( k = 0 ; k < 3 ; k++ )
+        pivot[k] = c->loc[k] + c->h[k]/2;
+    
+    /* Split along the x-axis. */
+    i = 0; j = c->count - 1;
+    while ( i <= j ) {
+        while ( i <= c->count-1 && parts[i].x[0] <= pivot[0] )
+            i += 1;
+        while ( j >= 0 && parts[j].x[0] > pivot[0] )
+            j -= 1;
+        if ( i < j ) {
+            temp = parts[i]; parts[i] = parts[j]; parts[j] = temp;
+            }
+        }
+    for ( k = 0 ; k <= j ; k++ )
+        if ( parts[k].x[0] > pivot[0] )
+            error( "cell_split: sorting failed." );
+    for ( k = i ; k < c->count ; k++ )
+        if ( parts[k].x[0] < pivot[0] )
+            error( "cell_split: sorting failed." );
+    left[1] = i; right[1] = c->count - 1;
+    left[0] = 0; right[0] = j;
+    
+    /* Split along the y axis, twice. */
+    for ( k = 1 ; k >= 0 ; k-- ) {
+        i = left[k]; j = right[k];
+        while ( i <= j ) {
+            while ( i <= right[k] && parts[i].x[1] <= pivot[1] )
+                i += 1;
+            while ( j >= left[k] && parts[j].x[1] > pivot[1] )
+                j -= 1;
+            if ( i < j ) {
+                temp = parts[i]; parts[i] = parts[j]; parts[j] = temp;
+                }
+            }
+        for ( kk = left[k] ; kk <= j ; kk++ )
+            if ( parts[kk].x[1] > pivot[1] ) {
+                printf( "cell_split: ival=[%i,%i], i=%i, j=%i.\n" , left[k] , right[k] , i , j );
+                error( "sorting failed (left)." );
+                }
+        for ( kk = i ; kk <= right[k] ; kk++ )
+            if ( parts[kk].x[1] < pivot[1] )
+                error( "sorting failed (right)." );
+        left[2*k+1] = i; right[2*k+1] = right[k];
+        left[2*k] = left[k]; right[2*k] = j;
+        }
+
+    /* Split along the z axis, four times. */
+    for ( k = 3 ; k >= 0 ; k-- ) {
+        i = left[k]; j = right[k];
+        while ( i <= j ) {
+            while ( i <= right[k] && parts[i].x[2] <= pivot[2] )
+                i += 1;
+            while ( j >= left[k] && parts[j].x[2] > pivot[2] )
+                j -= 1;
+            if ( i < j ) {
+                temp = parts[i]; parts[i] = parts[j]; parts[j] = temp;
+                }
+            }
+        for ( kk = left[k] ; kk <= j ; kk++ )
+            if ( parts[kk].x[2] > pivot[2] ) {
+                printf( "cell_split: ival=[%i,%i], i=%i, j=%i.\n" , left[k] , right[k] , i , j );
+                error( "sorting failed (left)." );
+                }
+        for ( kk = i ; kk <= right[k] ; kk++ )
+            if ( parts[kk].x[2] < pivot[2] ) {
+                printf( "cell_split: ival=[%i,%i], i=%i, j=%i.\n" , left[k] , right[k] , i , j );
+                error( "sorting failed (right)." );
+                }
+        left[2*k+1] = i; right[2*k+1] = right[k];
+        left[2*k] = left[k]; right[2*k] = j;
+        }
+        
+    /* Store the counts and offsets. */
+    for ( k = 0 ; k < 8 ; k++ ) {
+        c->progeny[k]->count = right[k] - left[k] + 1;
+        if ( c->progeny[k]->count < 0 )
+            abort();
+        c->progeny[k]->parts = &c->parts[ left[k] ];
+        }
+        
+    /* Verify a few sub-cells. */
+    /* for ( k = 0 ; k < c->progeny[0]->count ; k++ )
+        if ( c->progeny[0]->parts[k].x[0] > pivot[0] ||
+             c->progeny[0]->parts[k].x[1] > pivot[1] ||
+             c->progeny[0]->parts[k].x[2] > pivot[2] )
+            error( "Sorting failed (progeny=0)." );
+    for ( k = 0 ; k < c->progeny[1]->count ; k++ )
+        if ( c->progeny[1]->parts[k].x[0] > pivot[0] ||
+             c->progeny[1]->parts[k].x[1] > pivot[1] ||
+             c->progeny[1]->parts[k].x[2] <= pivot[2] )
+            error( "Sorting failed (progeny=1)." );
+    for ( k = 0 ; k < c->progeny[2]->count ; k++ )
+        if ( c->progeny[2]->parts[k].x[0] > pivot[0] ||
+             c->progeny[2]->parts[k].x[1] <= pivot[1] ||
+             c->progeny[2]->parts[k].x[2] > pivot[2] )
+            error( "Sorting failed (progeny=2)." ); */
+
+    }
+
+

src/cell.h

+/*******************************************************************************
+ * This file is part of GadgetSMP.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * 
+ ******************************************************************************/
+
+
+
+/* The queue timers. */
+enum {
+    cell_timer_none = 0,
+    cell_timer_tree,
+    cell_timer_count,
+    };
+extern ticks cell_timer[ cell_timer_count ];
+
+
+/* Structure to store the data of a single cell. */
+struct cell {
+
+    /* The cell location on the grid. */
+    double loc[3];
+    
+    /* The cell dimensions. */
+    double h[3];
+    
+    /* Max radii in this cell. */
+    double r_max;
+    
+    /* The depth of this cell in the tree. */
+    int depth, split;
+    
+    /* Nr of parts. */
+    int count;
+    
+    /* Pointers to the particle data. */
+    struct part *parts;
+    
+    /* Pointers for the sorted indices. */
+    struct entry *sort;
+    
+    /* Number of pairs associated with this cell. */
+    int nr_pairs;
+    
+    /* Pointers to the next level of cells. */
+    struct cell *progeny[8];
+    
+    /* Parent cell. */
+    struct cell *parent;
+    
+    /* The tasks computing this cell's sorts. */
+    struct task *sorts[14];
+    
+    /* Number of tasks this cell is waiting for and whether it is in use. */
+    int wait;
+    
+    /* Is the data of this cell being used in a sub-cell? */
+    int hold;
+    
+    /* Spin lock for various uses. */
+    lock_type lock;
+    
+    /* Linking pointer for "memory management". */
+    struct cell *next;
+
+    } __attribute__((aligned (64)));
+
+
+/* Function prototypes. */
+void cell_split ( struct cell *c  );
+int cell_locktree( struct cell *c );
+void cell_unlocktree( struct cell *c );

src/part.h

+/*******************************************************************************
+ * This file is part of GadgetSMP.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * 
+ ******************************************************************************/
+
+
+/* Some constants. */
+#define part_maxwait                    3
+#define part_maxunlock                  39
+
+
+/* Data of a single particle. */
+struct part {
+
+    /* Particle position. */
+    double x[3];
+    
+    /* Particle cutoff radius. */
+    float r;
+    
+    /* Particle time-step. */
+    float dt;
+    
+    /* Particle ID. */
+    int id;
+    
+    /* Number of pairwise interactions. */
+    double count, count_dh;
+    int icount;
+    
+    } __attribute__((aligned (32)));
+    
+

src/queue.c

+/*******************************************************************************
+ * This file is part of GadgetSMP.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * 
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <pthread.h>
+#include <math.h>
+#include <float.h>
+#include <limits.h>
+#include <omp.h>
+#include <sched.h>
+
+/* Local headers. */
+#include "cycle.h"
+#include "lock.h"
+#include "task.h"
+#include "cell.h"
+#include "queue.h"
+
+/* Error macro. */
+#define error(s) { printf( "%s:%s:%i: %s\n" , __FILE__ , __FUNCTION__ , __LINE__ , s ); abort(); }
+
+/* Define the timer macros. */
+#ifdef TIMER_VERBOSE
+    #define TIMER
+#endif
+#ifdef TIMER
+    #define TIMER_TIC ticks tic = getticks();
+    #define TIMER_TOC(t) timer_toc( t , tic )
+    #define TIMER_TIC2 ticks tic2 = getticks();
+    #define TIMER_TOC2(t) timer_toc( t , tic2 )
+    #ifndef INLINE
+    # if __GNUC__ && !__GNUC_STDC_INLINE__
+    #  define INLINE extern inline
+    # else
+    #  define INLINE inline
+    # endif
+    #endif
+    INLINE ticks timer_toc ( int t , ticks tic ) {
+        ticks d = (getticks() - tic);
+        __sync_add_and_fetch( &queue_timer[t] , d );
+        return d;
+        }
+#else
+    #define TIMER_TIC
+    #define TIMER_TOC(t)
+#endif
+
+
+/* Counter macros. */
+#ifdef COUNTER
+    #define COUNT(c) ( __sync_add_and_fetch( &queue_counter[ c ] , 1 ) )
+#else
+    #define COUNT(c)
+#endif
+
+
+/* The timers. */
+ticks queue_timer[ queue_timer_count ];
+
+/* The counters. */
+int queue_counter[ queue_counter_count ];
+
+        
+
+/** 
+ * @brief Initialize the given queue.
+ *
+ * @param q The #queue.
+ * @param size The maximum size of the queue.
+ * @param tasks List of tasks to which the queue indices refer to.
+ */
+ 
+void queue_init ( struct queue *q , int size , struct task *tasks ) {
+    
+    /* Allocate the task list. */
+    q->size = size;
+    if ( ( q->tid = (int *)malloc( sizeof(int) * size ) ) == NULL )
+        error( "Failed to allocate queue tids." );
+    q->tasks = tasks;
+        
+    /* Init counters. */
+    q->count = 0;
+    q->next = 0;
+    
+    /* Init the queue lock. */
+    if ( lock_init( &q->lock ) != 0 )
+        error( "Failed to init queue lock." );
+
+    }
+
+
+/**
+ * @brief Get a task free of dependencies and conflicts.
+ *
+ * @param q The task #queue.
+ * @param blocking Block until access to the queue is granted.
+ * @param keep Remove the returned task from this queue.
+ */
+ 
+struct task *queue_gettask ( struct queue *q , int blocking , int keep ) {
+
+    int k, tid = -1, qcount, *qtid = q->tid;
+    lock_type *qlock = &q->lock;
+    struct task *qtasks = q->tasks, *res = NULL;
+    TIMER_TIC
+    
+    /* If there are no tasks, leave immediately. */
+    if ( q->next >= q->count ) {
+        TIMER_TOC(queue_timer_gettask);
+        return NULL;
+        }
+
+    /* Main loop, while there are tasks... */
+    while ( q->next < q->count ) {
+    
+        /* Grab the task lock. */
+        // if ( blocking ) {
+            if ( lock_lock( qlock ) != 0 )
+                error( "Locking the task_lock failed.\n" );
+        //     }
+        // else {
+        //     if ( lock_trylock( qlock ) != 0 )
+        //         break;
+        //     }
+            
+        /* Loop over the remaining task IDs. */
+        qcount = q->count;
+        for ( k = q->next ; k < qcount ; k++ ) {
+        
+            /* Put a finger on the task. */
+            res = &qtasks[ qtid[k] ];
+            
+            /* Is this task blocked? */
+            if ( res->wait )
+                continue;
+            
+            /* Different criteria for different types. */
+            if ( res->type == tid_self || (res->type == tid_sub && res->cj == NULL) ) {
+                if ( res->ci->hold || cell_locktree( res->ci ) != 0 )
+                    continue;
+                }
+            else if ( res->type == tid_pair || (res->type == tid_sub && res->cj != NULL) ) {
+                if ( res->ci->hold || res->cj->hold || res->ci->wait || res->cj->wait )
+                    continue;
+                if ( cell_locktree( res->ci ) != 0 )
+                    continue;
+                if ( cell_locktree( res->cj ) != 0 ) {
+                    cell_unlocktree( res->ci );
+                    continue;
+                    }
+                }
+            
+            /* If we made it this far, we're safe. */
+            break;
+        
+            } /* loop over the task IDs. */
+            
+        /* Did we get a task? */
+        if ( k < qcount ) {
+        
+            /* Do we need to swap? */
+            if ( k != q->next )
+                COUNT(queue_counter_swap);
+        
+            /* get the task ID. */
+            tid = qtid[k];
+        
+            /* Remove the task? */
+            if ( keep ) {
+            
+                /* Bubble-up. */
+                q->count = qcount - 1;
+                for ( ; k < qcount - 1 ; k++ )
+                    qtid[k] = qtid[k+1];
+            
+                }
+                
+            /* No, leave it in the queue. */
+            else {
+            
+                TIMER_TIC2
+
+                /* Bubble-down the task. */
+                while ( k > q->next ) {