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

/* SolarSim.c */
#include "SolarSim.h"

typedef struct task {
        int i, j;
} Task, *TaskPtr;

int sock = 0;
long SIM_STEPS, DT, SAVE;

void
do_work(AState state) {
        int tid, nthreads;
        long i, j, t, ntasks, sim_step;
        int n = state->particles;
        vector *r = state->positions;
        vector *v = state->velocities;
        vector *a = state->accelerations;
        vector *ac;
        real *m = state->masses;
        TaskPtr tasks;

        tasks = malloc( n*(n-1)/2 * sizeof(Task) );
        t = 0;
        for( i = 0; i < n-1; i++ )
                for( j = i+1; j < n; j++ ) {
                        tasks[t].i = i;
                        tasks[t].j = j;
                        t++;
                }
        ntasks = t;

        #pragma omp parallel shared(nthreads)
        #pragma omp master
        nthreads = omp_get_num_threads();

        ac = malloc(nthreads * n * sizeof(vector));
        
        #pragma omp parallel shared(r,v,a,ac,m,n,tasks,ntasks,nthreads,DT,sim_step,SIM_STEPS,SAVE) private(i,j,t,tid) default(none)
        sim_step = 0;
        while( sim_step < SIM_STEPS ) {
                if( (SAVE != 0) && (sim_step % SAVE == 0) ) {
                        write_state(sock, state);
                }
                
                #pragma omp for schedule(static,n) //shared(ac,n,nthreads) private(i) default(none) nowait
                for( i = 0; i < n*nthreads; i++ ) {
                        ac[i].x = 0;
                        ac[i].y = 0;
                        ac[i].z = 0;
                }

                // compute accelerations. each thread uses ac[tid] vector
                #pragma omp for schedule(static) //shared(tasks, ntasks, r, v, ac, m, n) private(i,j,t,tid) default(none) nowait
                for( t = 0; t < ntasks; t++ ) {
                        i = tasks[t].i;
                        j = tasks[t].j;
                        tid = omp_get_thread_num();
                        compute_acc(r+i, r+j, v+i, v+j, ac+tid*n+i, ac+tid*n+j, m[i], m[j]);
                }

                // add up the accelerations computed by each thread
                #pragma omp parallel for schedule(static,1) //shared(a,ac,n,nthreads) private(i,j) default(none) nowait
                for( i = 0; i < n; i++ ) {
                        v_init(a+i,0,0,0);
                        for( j = 0; j < nthreads; j++ ) {
                                a[i].x += ac[j*n+i].x;
                                a[i].y += ac[j*n+i].y;
                                a[i].z += ac[j*n+i].z;
                        }
                }

                // update positions and speeds
                #pragma omp parallel for schedule(static,1) //shared(r,v,a,n,DT) private(i) default(none) nowait
                for( i = 0; i < n; i++ ) {
                        update_pos_vel(r+i, v+i, a+i, DT);
                }

                #pragma omp master
                sim_step ++;

                #pragma omp barrier
        }
        if( SAVE != 0 )
                write_state(sock, state);
}

int
main( int argc, char **argv ) {
        TState state;
        double before, after;
        int i, n=8;

        SIM_STEPS = atol(argv[3]);
        DT = atol(argv[4]);
        SAVE = atol(argv[5]);

        if( SAVE != 0 ) {
                sock = open( argv[2], O_WRONLY | O_CREAT | O_TRUNC, 644 );
        }

//      read_initial( &state, argv[1] );
        rand_data( &state, 100 );

        before = omp_get_wtime();
        do_work(&state);
        after = omp_get_wtime();

        if( SAVE != 0 )
                close(sock);

        #pragma omp parallel
        #pragma omp master
        printf("%d %lf\n", omp_get_num_threads(), after-before);

        return EXIT_SUCCESS;
}