source: proiecte/SolarSim/C/Parallel/src/SolarSim.c @ 152

Last change on this file since 152 was 152, checked in by (none), 14 years ago
File size: 2.6 KB
Line 
1/* SolarSim.c */
2#include "SolarSim.h"
3
4typedef struct task {
5        int i, j;
6} Task, *TaskPtr;
7
8int sock = 0;
9long SIM_STEPS, DT, SAVE;
10
11void
12do_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
85int
86main( 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
Note: See TracBrowser for help on using the repository browser.