cell.c 8.18 KB
Newer Older
1
/*******************************************************************************
2
 * This file is part of SWIFT.
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 * 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"
36
#include "timers.h"
37
38
#include "part.h"
#include "cell.h"
39
40
#include "error.h"
#include "inline.h"
41
42
43
44
45
46
47
48
49
50
51
52
53
54


/**
 * @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. */
55
    if ( c->hold || lock_trylock( &c->lock ) != 0 ) {
56
        TIMER_TOC(timer_locktree);
57
58
59
60
61
62
63
64
65
66
67
        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. */
68
        TIMER_TOC(timer_locktree);
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
        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 ) {
91
        TIMER_TOC(timer_locktree);
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
        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. */
107
        TIMER_TOC(timer_locktree);
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
        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 );
        
134
    TIMER_TOC(timer_locktree);
135
136
137
138
139
140
141
142
143
144
145
146
        
    }
    
    
/**
 * @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  ) {

147
    int i, j, k;
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
    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;
            }
        }
Pedro Gonnet's avatar
Pedro Gonnet committed
167
    /* for ( k = 0 ; k <= j ; k++ )
168
169
170
171
        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] )
Pedro Gonnet's avatar
Pedro Gonnet committed
172
            error( "cell_split: sorting failed." ); */
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
    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;
                }
            }
Pedro Gonnet's avatar
Pedro Gonnet committed
188
        /* for ( int kk = left[k] ; kk <= j ; kk++ )
189
190
191
192
            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)." );
                }
Pedro Gonnet's avatar
Pedro Gonnet committed
193
        for ( int kk = i ; kk <= right[k] ; kk++ )
194
            if ( parts[kk].x[1] < pivot[1] )
Pedro Gonnet's avatar
Pedro Gonnet committed
195
                error( "sorting failed (right)." ); */
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
        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;
                }
            }
Pedro Gonnet's avatar
Pedro Gonnet committed
212
        /* for ( int kk = left[k] ; kk <= j ; kk++ )
213
214
215
216
            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)." );
                }
Pedro Gonnet's avatar
Pedro Gonnet committed
217
        for ( int kk = i ; kk <= right[k] ; kk++ )
218
219
220
            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)." );
Pedro Gonnet's avatar
Pedro Gonnet committed
221
                } */
222
223
224
225
226
227
228
229
230
231
        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;
        c->progeny[k]->parts = &c->parts[ left[k] ];
        }
        
Pedro Gonnet's avatar
Pedro Gonnet committed
232
233
234
235
236
237
238
239
240
    /* Verify that _all_ the parts have been assigned to a cell. */
    /* for ( k = 1 ; k < 8 ; k++ )
        if ( &c->progeny[k-1]->parts[ c->progeny[k-1]->count ] != c->progeny[k]->parts )
            error( "Particle sorting failed (internal consistency)." );
    if ( c->progeny[0]->parts != c->parts )
        error( "Particle sorting failed (left edge)." );
    if ( &c->progeny[7]->parts[ c->progeny[7]->count ] != &c->parts[ c->count ] )
        error( "Particle sorting failed (right edge)." ); */
        
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
    /* 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)." ); */

    }