[152] | 1 | /* SolarSim.c */ |
---|
| 2 | #include "SolarSim.h" |
---|
| 3 | |
---|
| 4 | typedef struct task { |
---|
| 5 | int i, j; |
---|
| 6 | } Task, *TaskPtr; |
---|
| 7 | |
---|
| 8 | int sock = 0; |
---|
| 9 | long SIM_STEPS, DT, SAVE; |
---|
| 10 | |
---|
| 11 | void |
---|
| 12 | do_work(AState state) { |
---|
| 13 | int i, j, t, ntasks, tid, nthreads; |
---|
| 14 | long sim_step; |
---|
| 15 | int n = state->particles; |
---|
| 16 | vector *r = state->positions; |
---|
| 17 | vector *v = state->velocities; |
---|
| 18 | vector *a = state->accelerations; |
---|
| 19 | vector *ac; |
---|
| 20 | real *m = state->masses; |
---|
| 21 | TaskPtr tasks; |
---|
| 22 | |
---|
| 23 | tasks = malloc( n*(n-1)/2 * sizeof(Task) ); |
---|
| 24 | t = 0; |
---|
| 25 | for( i = 0; i < n-1; i++ ) |
---|
| 26 | for( j = i+1; j < n; j++ ) { |
---|
| 27 | tasks[t].i = i; |
---|
| 28 | tasks[t].j = j; |
---|
| 29 | t++; |
---|
| 30 | } |
---|
| 31 | ntasks = t; |
---|
| 32 | |
---|
| 33 | #pragma omp parallel shared(nthreads) |
---|
| 34 | { |
---|
| 35 | #pragma omp master |
---|
| 36 | nthreads = omp_get_num_threads(); |
---|
| 37 | } |
---|
| 38 | ac = malloc(nthreads * n * sizeof(vector)); |
---|
| 39 | |
---|
| 40 | sim_step = 0; |
---|
| 41 | while( sim_step < SIM_STEPS ) { |
---|
| 42 | if( (SAVE != 0) && (sim_step % SAVE == 0) ) { |
---|
| 43 | write_state(sock, state); |
---|
| 44 | } |
---|
| 45 | |
---|
| 46 | #pragma omp parallel for schedule(runtime) shared(ac,n,nthreads) private(i) default(none) |
---|
| 47 | for( i = 0; i < n*nthreads; i++ ) { |
---|
| 48 | ac[i].x = 0; |
---|
| 49 | ac[i].y = 0; |
---|
| 50 | ac[i].z = 0; |
---|
| 51 | } |
---|
| 52 | |
---|
| 53 | // compute accelerations. each thread uses ac[tid] vector |
---|
| 54 | #pragma omp parallel for schedule(runtime) shared(tasks, ntasks, r, v, ac, m, n) private(i,j,t,tid) default(none) |
---|
| 55 | for( t = 0; t < ntasks; t++ ) { |
---|
| 56 | i = tasks[t].i; |
---|
| 57 | j = tasks[t].j; |
---|
| 58 | tid = omp_get_thread_num(); |
---|
| 59 | // printf("%d computed (%d,%d)\n", tid, i, j); |
---|
| 60 | compute_acc(r+i, r+j, v+i, v+j, ac+tid*n+i, ac+tid*n+j, m[i], m[j]); |
---|
| 61 | } |
---|
| 62 | |
---|
| 63 | // add up the accelerations computed by each thread |
---|
| 64 | #pragma omp parallel for schedule(runtime) shared(a,ac,n,nthreads) private(i,j) default(none) |
---|
| 65 | for( i = 0; i < n; i++ ) { |
---|
| 66 | v_init(a+i,0,0,0); |
---|
| 67 | for( j = 0; j < nthreads; j++ ) { |
---|
| 68 | a[i].x += ac[j*n+i].x; |
---|
| 69 | a[i].y += ac[j*n+i].y; |
---|
| 70 | a[i].z += ac[j*n+i].z; |
---|
| 71 | } |
---|
| 72 | } |
---|
| 73 | |
---|
| 74 | // update positions and speeds |
---|
| 75 | #pragma omp parallel for schedule(runtime) shared(r,v,a,n,DT) private(i) default(none) |
---|
| 76 | for( i = 0; i < n; i++ ) { |
---|
| 77 | update_pos_vel(r+i, v+i, a+i, DT); |
---|
| 78 | } |
---|
| 79 | sim_step ++; |
---|
| 80 | } |
---|
| 81 | if( SAVE != 0 ) |
---|
| 82 | write_state(sock, state); |
---|
| 83 | } |
---|
| 84 | |
---|
| 85 | int |
---|
| 86 | main( int argc, char **argv ) { |
---|
| 87 | TState state; |
---|
| 88 | double before, after; |
---|
| 89 | int i; |
---|
| 90 | |
---|
| 91 | SIM_STEPS = atol(argv[3]); |
---|
| 92 | DT = atol(argv[4]); |
---|
| 93 | SAVE = atol(argv[5]); |
---|
| 94 | |
---|
| 95 | if( SAVE != 0 ) { |
---|
| 96 | sock = open( argv[2], O_WRONLY | O_CREAT | O_TRUNC, 644 ); |
---|
| 97 | } |
---|
| 98 | |
---|
| 99 | read_initial( &state, argv[1] ); |
---|
| 100 | |
---|
| 101 | before = omp_get_wtime(); |
---|
| 102 | do_work(&state); |
---|
| 103 | after = omp_get_wtime(); |
---|
| 104 | |
---|
| 105 | if( SAVE != 0 ) { |
---|
| 106 | i = 0; |
---|
| 107 | write( sock, &i, sizeof(int) ); |
---|
| 108 | close(sock); |
---|
| 109 | } |
---|
| 110 | |
---|
| 111 | #pragma omp parallel |
---|
| 112 | #pragma omp master |
---|
| 113 | printf("%d %lf\n", omp_get_num_threads(), after-before); |
---|
| 114 | |
---|
| 115 | return EXIT_SUCCESS; |
---|
| 116 | } |
---|
| 117 | |
---|