Optimisation d'addition de deux tableaux*constante

Effectivement le heu d’instruction SSE peut aider. Mais je ne suis pas convaincu.
La tendance actuel pou accélérer les calculs est de les paralléliser ,pour tirer partie des architectures multi-core. Une solution serait d'affecter un thread a une particule ou un groupes de particules, et de faire tourner les calculs en parallèles.
Je me permets de poster des fichiers (.c et .h) concernant une simulation de trajectoires orbitale (utilisant les lois de Newton) d'un ou plusieurs projectiles (ce qui doit etre proche du calcul particulaire).

Code:

positionX[i] = positionX[i] + vitesseX[i]*t

correspond à la partie /*INTREGATION*/ de la méthode projectile_compute_next_state.
Dans mon programme un thread est associe à un projectile (pthread_create() est appelé avec projectile_thread_computation en paramètre). Les calculs sont lancés pour un nombre fini d’itérations. Les données calculées sont sauvegardées, puis les threads sont relancés.
voici projectile.c (je sais je code en Franglais...)

Code:

1
2
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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
/*----------------------------------------------------------------------------*/
/*============================================================================*/
/*            INCLUDE                                                         */
/*============================================================================*/
/*----------------------------------------------------------------------------*/
#include <math.h>
#include <stdlib.h>
#include "projectile.h"
#include "space.h"
#include "simu.h"
#include "cte.h"
/*----------------------------------------------------------------------------*/
/*============================================================================*/
/*            PUBLIC METHODES IMPLEMENATION                                   */
/*============================================================================*/
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
/*  i_init_projectile                                                         */
/*----------------------------------------------------------------------------*/
int i_init_projectile(projectile_t * pst_p, void * pst_sp, uint64_t u64_id)
{
  if(pst_p == NULL || pst_sp == NULL){return -1;}
  pst_p->pst_sp = (space_t*)pst_sp;
  pst_p->u64_id = u64_id;
  pst_p->computed_vect = pst_p->computed_data;
  pst_p->current_vect = pst_p->current_data;
  pst_p->prev_vect = pst_p->prev_data;
  pst_p->pprev_vect = pst_p->pprev_data;
  return EXIT_SUCCESS;
}/* fin i_init_projectile*/
/*----------------------------------------------------------------------------*/
/* projectile_compute_next_state                                              */
/*----------------------------------------------------------------------------*/
void projectile_compute_next_state(projectile_t * pst_p, double dt_sur_24)
{
 
  register double current_x;
  register double current_y;
 
 
  double dist_s_p2;
  double dist_s_p;
  double cos_teta_S;
  double sin_teta_S;
 
  double px_relatif_earth;
  double py_relatif_earth;
  double dist_e_p2;
  double dist_e_p;
  double cos_teta_E;
  double sin_teta_E;
 
  double px_relatif_moon;
  double py_relatif_moon;
  double dist_m_p2;
  double dist_m_p;
  double cos_teta_M;
  double sin_teta_M;
 
  double px_relatif_mars;
  double py_relatif_mars;
  double dist_mars_p2 ;
  double dist_mars_p;
  double cos_teta_Mars;
  double sin_teta_Mars;
 
  double px_proj;
  double py_proj;
  double dist_proj_2 ;
  double dist_proj_p;
  double cos_teta_proj;
  double sin_teta_proj;
 
  double computed_ax;
  double computed_ay;
 
  double pp_vy;
  double p_vy;
  double cur_vy;
 
  double pp_vx;
  double p_vx;
  double cur_vx;
 
  double pp_ay;
  double p_ay ;
  double cur_ay;
 
  double pp_ax;
  double p_ax;
  double cur_ax;
 
  double computed_vx;
  double computed_vy;
 
  space_t * sp = (space_t *)(pst_p->pst_sp);
 
  current_x = projectile_getX(pst_p);
  current_y = projectile_getY(pst_p);
 
 
  /*SOLEIL*/
  /*distance projectile soleil au carré en m²*/
  dist_s_p2 = current_x *current_x + current_y*current_y;
  dist_s_p = sqrt(dist_s_p2);
  cos_teta_S = current_x /dist_s_p;
  sin_teta_S = current_y /dist_s_p;
 
  /*TERRE*/
  /*distance projectile terre au carré en m²*/
  px_relatif_earth = current_x - sp->st_earth.current_vect[S_X];
  py_relatif_earth = current_y - sp->st_earth.current_vect[S_Y];
  dist_e_p2 =   px_relatif_earth * px_relatif_earth +
                py_relatif_earth * py_relatif_earth;
  dist_e_p = sqrt(dist_e_p2);
  cos_teta_E = px_relatif_earth /dist_e_p;
  sin_teta_E = py_relatif_earth /dist_e_p;
 
  /*MOON*/
 
  /*distance projectile terre au carré en m²*/
  px_relatif_moon = current_x - sp->st_moon.current_vect[S_X];
  py_relatif_moon = current_y - sp->st_moon.current_vect[S_Y];
  dist_m_p2 = px_relatif_moon*px_relatif_moon +
              py_relatif_moon*py_relatif_moon ;
  dist_m_p = sqrt(dist_m_p2);
  cos_teta_M = px_relatif_moon / dist_m_p;
  sin_teta_M = py_relatif_moon / dist_m_p;
 
  /*MARS*/
  /*distance projectile terre au carré en m²*/
  px_relatif_mars = current_x - sp->st_mars.current_vect[S_X];
  py_relatif_mars = current_y - sp->st_mars.current_vect[S_Y];
  dist_mars_p2 = px_relatif_mars * px_relatif_mars +
                 py_relatif_mars * py_relatif_mars;
  dist_mars_p = sqrt(dist_mars_p2);
  cos_teta_Mars = px_relatif_mars /dist_mars_p;
  sin_teta_Mars = py_relatif_mars /dist_mars_p;
 
 
  /*acceleration du projectil en fonction des satellite x en m/s/s*/
  computed_ax = -G*MS/dist_s_p2*cos_teta_S /*attration du soleil x*/
                -G*ME/dist_e_p2*cos_teta_E /*attration de la terre x*/
                -G*MM/dist_m_p2*cos_teta_M /*attraction de la lune x*/
                -G*MMARS/dist_mars_p2*cos_teta_Mars; /*attration de Mars x*/
 
  computed_ay = -G*MS/dist_s_p2*sin_teta_S /*attration du soleil y*/
                -G*ME/dist_e_p2*sin_teta_E /*attration de la terre y*/
                -G*MM/dist_m_p2*sin_teta_M/*attraction de la lune y*/
                -G*MMARS/dist_mars_p2*sin_teta_Mars; /*attration de Mars y*/
 
  /* interractions entre les projectiles*/
  uint64_t u64_i;
  uint64_t u64_nb_prjectile = sp->u64_nb_projectile;
  for(u64_i = 0 ; u64_i < u64_nb_prjectile; u64_i++)
  {/* parcours des projectiles*/
    if(u64_i != projectile_getID(pst_p))
    {/* projectile !=*/
      const projectile_t * pst_p_next = pst_space_get_projectile(sp,u64_i);
      px_proj = current_x - projectile_getX(pst_p_next);
      py_proj = current_y - projectile_getY(pst_p_next);
      dist_proj_2 = px_proj*px_proj + py_proj*py_proj;
      dist_proj_p = sqrt(dist_proj_2);
      cos_teta_proj = px_proj/dist_proj_p;
      sin_teta_proj = py_proj/dist_proj_p;
      computed_ax -= G*projectile_getMASSE(pst_p_next)/dist_proj_2*cos_teta_proj;
      computed_ay -= G*projectile_getMASSE(pst_p_next)/dist_proj_2*sin_teta_proj;
    }
  }/* fin parcours des projectiles*/
 
  /*INTEGRATION*/
  pp_vy  = pst_p->pprev_vect[P_VY];
  p_vy   = pst_p->prev_vect[P_VY];
  cur_vy = pst_p->current_vect[P_VY];
 
  pp_vx  = pst_p->pprev_vect[P_VX];
  p_vx   = pst_p->prev_vect[P_VX];
  cur_vx = pst_p->current_vect[P_VX];
 
  pp_ay  = pst_p->pprev_vect[P_AY];
  p_ay   = pst_p->prev_vect[P_AY];
  cur_ay = pst_p->current_vect[P_AY];
 
  pp_ax  = pst_p->pprev_vect[P_AX];
  p_ax   = pst_p->prev_vect[P_AX];
  cur_ax = pst_p->current_vect[P_AX];
 
  /*Vitesse du projectile suivant x en m/s*/
  computed_vx = pst_p->current_vect[P_VX] +
    (9.0*computed_ax + 19.0*cur_ax - 5.0*p_ax + pp_ax)*dt_sur_24;
  /*Vitesse du projectile suivant y en m/s*/
  computed_vy = pst_p->current_vect[P_VY] +
    (9.0*computed_ay + 19.0*cur_ay - 5.0*p_ay + pp_ay)*dt_sur_24;
  /*Position du projectile suivant x en m*/
  pst_p->computed_vect[P_X] = current_x +
    (9.0*computed_vx + 19.0*cur_vx - 5.0*p_vx + pp_vx)*dt_sur_24;
  /*Position du projectile suivant y en m*/
  pst_p->computed_vect[P_Y] = current_y +
    (9.0*computed_vy + 19.0*cur_vy - 5.0*p_vy + pp_vy)*dt_sur_24;
 
 
  pst_p->computed_vect[P_VX] = computed_vx;
  pst_p->computed_vect[P_VY] = computed_vy;
  pst_p->computed_vect[P_AX] = computed_ax;
  pst_p->computed_vect[P_AY] = computed_ay;
 
}
/*----------------------------------------------------------------------------*/
/*  projectile_setX                                                           */
/*----------------------------------------------------------------------------*/
void projectile_setX(projectile_t * p, double x)
{
    p->current_vect[P_X] = x;
    p->prev_vect[P_X] = x;
    p->pprev_vect[P_X] = x;
    projectile_setVX(p,0.0);
}/* fin projectile_setX*/
/*----------------------------------------------------------------------------*/
/*  projectile_setY                                                           */
/*----------------------------------------------------------------------------*/
void projectile_setY(projectile_t * p, double y)
{
    p->current_vect[P_Y] = y;
    p->prev_vect[P_Y] = y;
    p->pprev_vect[P_Y] = y;
    projectile_setVY(p,0.0);
}/* fin projectile_setY */
/*----------------------------------------------------------------------------*/
/*  projectile_setVX                                                          */
/*----------------------------------------------------------------------------*/
void projectile_setVX(projectile_t * p, double vx)
{
    p->current_vect[P_VX] = vx;
    p->prev_vect[P_VX] = vx;
    p->pprev_vect[P_VX] = vx;
}/* fin projectile_setVX*/
/*----------------------------------------------------------------------------*/
/*  i_set_projectile_state                                                    */
/*----------------------------------------------------------------------------*/
void projectile_setVY(projectile_t * p, double vy)
{
    p->current_vect[P_VY] = vy;
    p->prev_vect[P_VY] = vy;
    p->pprev_vect[P_VY] = vy;
}/* projectile_setVY*/
/*----------------------------------------------------------------------------*/
/* void projectile_setMASSE(projectile_t * p, double masse);                  */
/*----------------------------------------------------------------------------*/
void projectile_setMASSE(projectile_t * p, double masse)
{
  p->current_vect[P_MASSE] = masse;
  p->prev_vect[P_MASSE] = masse;
  p->pprev_vect[P_MASSE] = masse;
}/* projectile_setMASSE*/
/*----------------------------------------------------------------------------*/
/*  set_projectile_state                                                      */
/*----------------------------------------------------------------------------*/
void set_projectile_state(const projectile_t * pst_p,
                          st_projectile_state_t * pst_st)
{
  pst_st->u64_id  = projectile_getID(pst_p);
  pst_st->d_X     = projectile_getX(pst_p);
  pst_st->d_Y     = projectile_getY(pst_p);
  pst_st->d_vX    = projectile_getvX(pst_p);
  pst_st->d_vY    = projectile_getvY(pst_p);
}/* fin i_set_projectile_state*/
 
/*----------------------------------------------------------------------------*/
/* projectile_set                                                             */
/*----------------------------------------------------------------------------*/
void projectile_set(projectile_t * pst_p,
                    double x,
                    double vx,
                    double y,
                    double vy)
{
  pthread_log_t * pst_log = pst_simu_get_log(pst_space_get_sim(
    pst_projectile_get_space(pst_p)));
  projectile_setX(pst_p,x);
  projectile_setVX(pst_p,vx);
  projectile_setY(pst_p,y);
  projectile_setVY(pst_p,vy);
  i_log(pst_log,"Projectile %lu set x:%E m, vx:%E m/s, y:%E m, vy:%E m/s",
          projectile_getID(pst_p),
          projectile_getX(pst_p),projectile_getvX(pst_p),
          projectile_getY(pst_p),projectile_getvY(pst_p));
}/*fin projectile_set*/
/*----------------------------------------------------------------------------*/
/* projectile_thread_computation                                              */
/*----------------------------------------------------------------------------*/
void * projectile_thread_computation(void *pst_args)
{
  projectile_t * pst_p = (projectile_t *)pst_args;
  simu_t * pst_sim = pst_space_get_sim(pst_projectile_get_space(pst_p));
  double dt_sur_24 = pst_sim->dt_sur_24;
  uint64_t u64_nb_iter = pst_sim->u64_nb_iter;
  pthread_t th_this = pst_sim->pst_th_computation[projectile_getID(pst_p)+1];
  uint64_t u64_i;
  //bool test;
  for(u64_i = 0; u64_i < u64_nb_iter ; u64_i++)
  {
    //met la priorité au max
    pthread_setschedprio(th_this,sched_get_priority_max(SCHED_FIFO));
    projectile_compute_next_state(pst_p,dt_sur_24);
 
    //met la priorité au min
    pthread_setschedprio(th_this,sched_get_priority_min(SCHED_FIFO));
    simu_threads_wait_synchro(pst_sim);
 
    //met la priorité au max
    pthread_setschedprio(th_this,sched_get_priority_max(SCHED_FIFO));
    projectile_update_vect(pst_p);
 
    //met la priorité au min
    pthread_setschedprio(th_this,sched_get_priority_min(SCHED_FIFO));
    simu_threads_wait_synchro(pst_sim);
 
  }/* fin des iterations*/
   pthread_exit(NULL);
}/* fin projectile_thread_computation*/

et projectile .h

Code:

1
2
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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
/* 
 * File:   projectile.h
 * Author: mathieu
 *
 * Created on 29 mai 2010, 16:33
 */
 
#ifndef _PROJECTILE_H
#define	_PROJECTILE_H
 
#ifdef	__cplusplus
extern "C" {
#endif
/*----------------------------------------------------------------------------*/
/*============================================================================*/
/*              INCLUDE                                                       */
/*============================================================================*/
/*----------------------------------------------------------------------------*/
#include <stdint.h>
#include <pthread.h>
 
/*----------------------------------------------------------------------------*/
/*============================================================================*/
/*              CONSTANTES                                                    */
/*============================================================================*/
/*----------------------------------------------------------------------------*/
#define P_X 0
#define P_VX 1
#define P_Y 2
#define P_VY 3
#define P_AX 4
#define P_AY 5
#define P_MASSE 6
#define P_VECT_LENGTH 7
/*----------------------------------------------------------------------------*/
/*============================================================================*/
/*              TYPES AND DEFINTIONS                                          */
/*============================================================================*/
/*----------------------------------------------------------------------------*/
typedef struct
{
    double    pprev_data[P_VECT_LENGTH];/*le vecteur d'etat t-2*/
    double    prev_data[P_VECT_LENGTH];/*le vecteur d'etat t-1*/
    double    current_data[P_VECT_LENGTH];/*le vecteur d'etat t courant*/
    double    computed_data[P_VECT_LENGTH];/*le vecteur d'etat calculé*/
    uint64_t  u64_id;
    double *  pprev_vect;/* reference sur le vecteur d'etat t-2*/
    double *  prev_vect;/* reference sur le vecteur d'etat t-1*/
    double *  current_vect;/* reference sur le vecteur d'etat t courant*/
    double *  computed_vect;/* reference sur le vecteur d'etat calculé*/
    void   *  pst_sp;/* reference sur l'espace de simulation space*/
} projectile_t;
 
typedef struct
{
    uint64_t    u64_id;
    double      d_X;
    double      d_vX;
    double      d_Y;
    double      d_vY;
}st_projectile_state_t;
/*----------------------------------------------------------------------------*/
/*============================================================================*/
/*              PUBLIC static inline METHODS IMPLEMENATION                    */
/*============================================================================*/
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
/* projectile_getID                                                           */
/*----------------------------------------------------------------------------*/
static inline uint64_t projectile_getID(const projectile_t * pst_p)
{
    return pst_p->u64_id;
}/* fin projectile_getID*/
/*----------------------------------------------------------------------------*/
/* projectile_getX                                                            */
/*----------------------------------------------------------------------------*/
static inline double projectile_getX(const projectile_t * pst_p)
{
    return pst_p->current_vect[P_X];
}/* fin projectile_getX*/
/*----------------------------------------------------------------------------*/
/* projectile_getY                                                            */
/*----------------------------------------------------------------------------*/
static inline double projectile_getY(const projectile_t * pst_p)
{
    return pst_p->current_vect[P_Y];
}/* fin projectile_getY*/
/*----------------------------------------------------------------------------*/
/* projectile_getvX                                                           */
/*----------------------------------------------------------------------------*/
static inline double projectile_getvX(const projectile_t * pst_p)
{
    return pst_p->current_vect[P_VX];
}/* fin projectile_gevtX*/
/*----------------------------------------------------------------------------*/
/* projectile_getvY                                                           */
/*----------------------------------------------------------------------------*/
static inline double projectile_getvY(const projectile_t * pst_p)
{
    return pst_p->current_vect[P_VY];
}/* fin projectile_getvY*/
/*----------------------------------------------------------------------------*/
/* projectile_getMASSE                                                           */
/*----------------------------------------------------------------------------*/
static inline double projectile_getMASSE(const projectile_t * pst_p)
{
    return pst_p->current_vect[P_MASSE];
}/* fin projectile_getMASSE*/
/*----------------------------------------------------------------------------*/
/* projectile_update_vect                                                     */
/*----------------------------------------------------------------------------*/
static inline void projectile_update_vect(projectile_t * pst_p)
{
    pst_p->pprev_vect   = pst_p->prev_vect;
    pst_p->prev_vect    = pst_p->current_vect;
    pst_p->current_vect = pst_p->computed_vect;
}/* fin projectile_update_vect*/
/*----------------------------------------------------------------------------*/
/*============================================================================*/
/*              PUBLIC METHODS DECLARATION                                    */
/*============================================================================*/
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
/*============================================================================*/
int i_init_projectile(projectile_t * pst_p, void * pst_sp, uint64_t u64_id);
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void set_projectile_state(const projectile_t * pst_p,
                          st_projectile_state_t * pst_st);
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void projectile_compute_next_state(projectile_t * pst_p, double dt_sur_24);
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void projectile_setX(projectile_t * p, double x);
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void projectile_setY(projectile_t * p, double y);
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void projectile_setVX(projectile_t * p, double vx);
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void projectile_setVY(projectile_t * p, double vy);
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void projectile_setMASSE(projectile_t * p, double masse);
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void projectile_set(projectile_t * pst_p,
                    double x,
                    double vx,
                    double y,
                    double vy);
/*----------------------------------------------------------------------------*/
/*============================================================================*/
void * projectile_thread_computation(void *pst_args);
 
#ifdef	__cplusplus
}
#endif
 
#endif	/* _PROJECTILE_H */

Je me permet egalements de faire un hors sujet en ajoutant que les processeurs CISC ca sert vraiment a rien puisqu'un compilateur n'est pas capable de generer une instrucition SSE, est que en plus d'avoir un tas d'instruction inutiles les processeurs CISC n'ont pas d'instruction de type mulacc (multiply and accumulate) qui permet de faire a += b*c, qui est exactement le type d'instruction d'ont ou a besoin pour faire des intregales....

bon courage.