Discussion :

[Matthieu Brucher]

Bon j'ai testé de mon coté, et j'ai les même comportement que toi avec le compilo icc. Inquiétant que pour la portabilité ça ne soit pas le cas, car avec g++ ça ne fonctionne pas! à moins que le __restrict__ ne soit pas reconnu dans le g++ car c'est de la norme C99 et que ce soit vaguement rempalcé par autre chose?

Tu as compilé avec le flag -std=c99 ?

J'ai du mal à voir en quoi __restrict__ permet de gagner un facteur 2 sur ce genre de boucle. As tu regardé le code assembleur généré dans les deux cas ?

Avec Visual2005, __restrict ne change rien, par contre compiler le code en SSE(2) fait perdre un facteur 2 !

Il y a moins de load d'effectués et il peut vectoriser la boucle interne en conséquence, puisqu'il sait que la modification de P ne touchera pas Q.

[loufoque]

Déjà ça irait pas plus vite si tu fusionnais les boucles et que tu évitais d'appeler ta macro XYZ tout le temps ?
Voir comment fait Adobe.GIL pour itérer les pixels par exemple.

[white_tentacle]

Déjà ça irait pas plus vite si tu fusionnais les boucles et que tu évitais d'appeler ta macro XYZ tout le temps ?

En l'occurrence, il m'a semblé que remplacer le n*n par 10000 avait un gain mesurable, de l'ordre de 0,2 s.

J'en suis très étonné, parce que j'ai toujours pensé que ce type d'optimisations (propagation de constantes) était géré correctement par les compilos (d'autant plus quand la constante est connu à la compilation).

[loufoque]

Si compute est pas inlinée, tu peux pas trop propager les constantes...

[white_tentacle]

En l'occurrence, N est constant dans compute, donc N*N est une constante qui pourrait tout à fait être propagée...

Changer la signature de compute (pour passer en const int) ne change d'ailleurs rien à ce niveau là.

[reptils]

de toute façon que ce soit en C/C++ ou Fortran ça change rien par rapport au NxN.

En revanche, je me demande si au niveau des accès mémoire il est plus pertinent de le faire comme tu le fais via ta macro, ou en passant par des pointeurs de pointeurs avec une bloc alloué de manière contiguë (N*N*N) et le dé-référencement qui va bien en [i][j][k]?

[Matthieu Brucher]

Déjà ça irait pas plus vite si tu fusionnais les boucles et que tu évitais d'appeler ta macro XYZ tout le temps ?
Voir comment fait Adobe.GIL pour itérer les pixels par exemple.

Sachant que ce sont les mêmes performances en Fortran, chez qui la boucle est on ne peut plus optimisée, le compilateur fait grosso modo ce qu'il faut en Fortran.
Maintenant, fusionner la boucle, sachant que le stencil est 3D, qu'en réalité, la boucle est bien plus complexe, comme je l'ai dit, c'est inutile et une perte de temps (au moins là, on voit ce que le stencil fait). Et une complexification qui va engendrer une perte de temps (déjà testé sur un autre projet avec des boucles plus petites, le recours à des itérateurs est une perte de temps épouvantable, puisque le problème se situe au niveau du chargement des données qui doit être optimisé).

Encore une fois, le problème n'est pas dans le calcul, mais dans le chargement des données qui doit être optimisé pour se rapprocher du Fortran. Le changement de boucle va te faire gagner, allez, 10% dans le meilleur des cas, sachant que l'utilisation correcte de la mémoire te procure 50%. Comme toujours, mieux vaut se focaliser sur ce qui peut faire gagner le plus, d'autant qu'on sait qu'on peut arriver à ces 50% (démontré avec la version Fortran qui fait strictement la même chose).

de toute façon que ce soit en C/C++ ou Fortran ça change rien par rapport au NxN.

En revanche, je me demande si au niveau des accès mémoire il est plus pertinent de le faire comme tu le fais via ta macro, ou en passant par des pointeurs de pointeurs avec une bloc alloué de manière contiguë (N*N*N) et le dé-référencement qui va bien en [i][j][k]?

A mon avis, c'est plus rapide ainsi, puisqu'il voit bien les strides dans le stencil. Dans le cas d'un déréférencement, il doit chercher pour chaque stencil les nouveaux éléments dans les 3 directions, au cas où, ce qui n'est pas le cas ici puisqu'il connaît le stride.

[reptils]

A mon avis, c'est plus rapide ainsi, puisqu'il voit bien les strides dans le stencil. Dans le cas d'un déréférencement, il doit chercher pour chaque stencil les nouveaux éléments dans les 3 directions, au cas où, ce qui n'est pas le cas ici puisqu'il connaît le stride.

J'ai testé et je passe de 1.7 à 1.6 en passant de ton stride au dé-référencement
Ceci dit j'ai un peu lutté pour bien faire mon dé-référencement en mode tenseur.

EDIT: je pense que l'accès en i+n*j+n*n*k est constant quoi qu'il arrive pour N, en revanche en ce qui concerne les pointeurs de pointeurs de pointeurs, les perfs risquent de s'effondrer si N augmente

[Matthieu Brucher]

Bon, a priori, il est nécessaire de passer par un appel à une fonction dans laquelle on met les boucles et dont les arguments sont les différents tableaux utilisés.
Reste la question de la portabilité du mot-clé "restrict". Est-ce qu'il est prévu qu'il passe dans la norme du C++ sous une forme ou une autre ?

[Jean-Marc.Bourguet]

Est-ce qu'il est prévu qu'il passe dans la norme du C++ sous une forme ou une autre ?

Pas a ma connaissance (je ne l'ai pas trouve dans le CD de C++0X).

[Arzar]

J'ai testé et je passe de 1.7 à 1.6 en passant de ton stride au dé-référencement
Ceci dit j'ai un peu lutté pour bien faire mon dé-référencement en mode tenseur.

J'ai moi aussi essayé, avec ce code (adapté des NRC, je garantis pas qu'il soit juste)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 
template<typename T>
T*** alloc_3DArray(int nrow, int ncol, int ndep)
{
	typedef T*** ptr3;
	typedef T** ptr2;
	typedef T*  ptr;
 
	//allocate pointers to pointers to rows 
	ptr3 m = new ptr2[nrow];
 
	// allocate pointers to rows and set pointers to them
	m[0]= new ptr[nrow * ncol];
 
	// allocate rows and set pointers to them 
	m[0][0]= new T[nrow * ncol * ndep];
 
	for(int j = 1 ; j < ncol ; j++)
		m[0][j] = m[0][j-1] + ndep;
 
	for(int i = 1 ; i < nrow ;i++)
	{
		m[i] = m[i-1] + ncol;
		m[i][0] = m[i-1][0] + ncol * ndep;
		for(int j = 1; j < ncol; j++)
			m[i][j] = m[i][j-1] + ndep;
	}
	return m;
}
 
template <typename T>
void release_3DArray(T*** m)
{
	delete m[0][0];
	delete m[0];
	delete m;
	return;
}

Resultat : passage de 800ms à 880ms.

[Matthieu Brucher]

Pas a ma connaissance (je ne l'ai pas trouve dans le CD de C++0X).

Je trouve ça dommage. Au moment où on ajoute dans les compilateurs des outils pour vectoriser le code automatiquement, ne pas permettre d'ajouter cette indication, c'est une erreur. Sans elle, le compilateur doit être bien plus prudent.
Une fois de plus, il va falloir s'appuyer sur une extension proposée par certains compilateurs pour gagner un facteur 2. Et un facteur 2, c'est énorme. Pas la peine de se creuser la tête à paralélliser le code ou le vectoriser quand on peut facilement gagner un tel facteur (et comme les compilateurs ont souevnt aussi un frontend Fortran, ce n'est pas ajouter ce mot-clé qui va les pénaliser dans leurs développements).

[yan]

ce ne serait pas un problème de ce type?

http://www.developpez.net/forums/d53...acces-memoire/

[Matthieu Brucher]

Pas tout à fait, non Dans ton cas, je ne suis pas persuadé des effets de l'aliasing.

[screetch]

si je fais

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
	float* P = new float[N*N*N];
	float* Q = new float[N*N*N+16]+16;

dans ta premiere fonction compute, le code est presque 2 fois plus rapide.

[Arzar]

En effet: 800ms -> 466ms.
Monsieur screetch, vous comprenez bien que vous ne vous en sortirez pas comme ça, il va falloir se mettre à table maintenant : Qu'est ce que c'est que cette diablerie !?

[screetch]

Store to load forward dependency

CF le post de mongaulois
l'endroit ou l'on veut ecrire en mémoire est aligné avec l'endroit ou on veut lire
du coup le CPU ne sait pas, a cause de l'alignement, si la mémoire a cet endroit est lisible ou pas.
C'est du au fait que le calcul dans la boucle est très simple (basiquement, lecture/ecriture) et cela ne veut pas dire que ca arriverait dans un cas réel.

[reptils]

aucun changement de mon coté?

[screetch]

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
 
static void compute(float* P, float* Q, int N)
{
	float* buffer = (float*) alloca(N*sizeof(float));
	for(int k = 4; k < N-4; ++k)
	{
		for(int j = 4; j < N-4; ++j)
		{
			for(int i = 4; i < N-4; ++i)
			{
				buffer[i] = P[XYZ(i, j, k, N)]
				+ (P[XYZ(i+1, j, k, N)] - P[XYZ(i-1, j, k, N)])
					+ (P[XYZ(i+2, j, k, N)] - P[XYZ(i-2, j, k, N)])
					+ (P[XYZ(i+3, j, k, N)] - P[XYZ(i-3, j, k, N)])
					+ (P[XYZ(i+4, j, k, N)] - P[XYZ(i-4, j, k, N)])
					+ (P[XYZ(i, j+1, k, N)] - P[XYZ(i, j-1, k, N)])
					+ (P[XYZ(i, j+2, k, N)] - P[XYZ(i, j-2, k, N)])
					+ (P[XYZ(i, j+3, k, N)] - P[XYZ(i, j-3, k, N)])
					+ (P[XYZ(i, j+4, k, N)] - P[XYZ(i, j-4, k, N)])
					+ (P[XYZ(i, j, k+1, N)] - P[XYZ(i, j, k-1, N)])
					+ (P[XYZ(i, j, k+2, N)] - P[XYZ(i, j, k-2, N)])
					+ (P[XYZ(i, j, k+3, N)] - P[XYZ(i, j, k-3, N)])
					+ (P[XYZ(i, j, k+4, N)] - P[XYZ(i, j, k-4, N)]);
			}
			for(int i = 4; i < N-4; ++i)
			{
				Q[XYZ(i,j,k,N)] = buffer[i];
			}
		}
	}
}

ca donne quoi chez vous ca ?
chez moi ca résout le probleme d'aliasing et de STLFD

[Matthieu Brucher]

Effectivement, c'est intéressant, j'obtiens 1.6s (pas de facteur 2 en revanche par rapport au Fortran). Le problème, c'est qu'il n'est pas envisageable d'allouer à la volée à chaque fois des tableaux de taille très importante (plusieurs dizaines de Mo), surtout qu'il n'y en a plus qu'un en réalité

Il me semble qu'il s'agit d'un problème de dépendance (c'est ce que j'ai appelé banque dans l'autre thread) qu'on a déjà observé chez nous (reptils> demande à SGI, ils te montreront un exemple très frappant ).