Discussion :

[gobelin88]

Bonjour tous,

Voila j'aimerai optimiser le code ci dessous, pour donner une idée une multiplication en moins dans la 3eme boucle for c'est 1s de temps de calcul mon rêve serait de pouvoir supprimer la racine carre ^^ sachant que je connais a^2 je voudrais a^3

Notes:
- les vecteurs sont des classes de Eigen3
- TmpParamForBS c'est des constantes pré-calculées
- Cette fonction réalise une intégration en type double.

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
39
40
41
42
43
44
45
46
47
48
 
void Truc::get_precalc_B_optimized(const Vector3d & m, Vector3d & b,const TmpParamForBS & tmp)
{
    Vector3d OM=Mrot.transpose()*(m-pos);
    b=Vector3d(0,0,0);
 
    double omx=OM[0];
    double omy=OM[1];
    double omz=OM[2];
    double omx_dl2_L2=omx-tmp.dl2_L2;
 
    std::vector<double> rx_list(tmp.rk);
    for(int k=0;k<tmp.rk;k++) rx_list[k]=(-k*tmp.dl+omx_dl2_L2);
    std::vector<double> rx2_list(tmp.rk);
    for(int k=0;k<tmp.rk;k++) rx2_list[k]=rx_list[k]*rx_list[k];
 
    for(int i=0;i<tmp.ri;i++)
    {
        double alpha=i*tmp.dtheta+tmp.dtheta2;
        double co=cos(alpha);
        double si=sin(alpha);
 
        for(int j=0;j<tmp.rj;j++)
        {
            Vector3d Btemp(0,0,0);
            double ray = j*tmp.de+tmp.R_de2;
            double omy_ray_c = omy-ray*co;
            double omz_ray_s = omz-ray*si;
            double si_r2_co_r1 = si*omz_ray_s+co*omy_ray_c;
            double omy_ray_c2_omz_ray_s2 = omy_ray_c*omy_ray_c+omz_ray_s*omz_ray_s;
 
            for(int k=0;k<tmp.rk;k++)
            {
                double rx=rx_list[k];
                double nr2=rx2_list[k]+omy_ray_c2_omz_ray_s2;
                double nr= sqrt(nr2);
                //if(nr>quantom_zero)
                {
                    Btemp+= Vector3d(si_r2_co_r1, -co*rx, -si*rx)/(nr2*nr);
                }
            }
 
            b+= (tmp.tmp1*j + tmp.tmp2)*Btemp;
        }
    }
    b*=Mu0_Pi4*tmp.J ;
    b=Mrot*b;
}

Je suis preneur de toute idée

merci d'avance pour vos suggestions.

[pepito3364]

Salut,

Je te conseillerai déjà de déclarer toutes tes variables double (ray, omy_ray_c, ...) en-dehors des boucles. Comme ça tu ne les déclare qu'une seule fois et pas, dans le pire des cas, tmp^3.

J'essayerai de calculer les indices 2 et 3 de Vector3d plus tot. Je sais pas si ce que je te propose est correct, j'ai fait vite fait.

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
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
 
void Truc::get_precalc_B_optimized(const Vector3d & m, Vector3d & b,const TmpParamForBS & tmp)
{
    Vector3d OM=Mrot.transpose()*(m-pos);
    b=Vector3d(0,0,0);
 
    double omx=OM[0];
    double omy=OM[1];
    double omz=OM[2];
    double omx_dl2_L2=omx-tmp.dl2_L2;
 
    std::vector<double> rx_list(tmp.rk);
    for(int k=0;k<tmp.rk;k++) 
		rx_list[k]=(-k*tmp.dl+omx_dl2_L2);
 
    std::vector<double> rx2_list(tmp.rk);
 
	for(int k=0;k<tmp.rk;k++) 
		rx2_list[k]=rx_list[k]*rx_list[k];
 
	//++++++
	std::vector<double> vecTmp1(tmp);
	std::vector<double> vecTmp2(tmp);
 
	double omy_ray_c = 0.0;
	double omz_ray_s = 0.0;
	double omy_ray_c2_omz_ray_s2 = 0.0;
	double ray = 0.0;
 
	double rx = 0.0;
	double nr2 = 0.0;
	double nr = 0.0;
        double ray = 0.0;
        double si_r2_co_r1 = 0.0;
	//++++++
 
	for(int i=0;i<tmp.ri;i++)
	{
		double alpha=i*tmp.dtheta+tmp.dtheta2;
		double co=cos(alpha);
		double si=sin(alpha);
 
		//++++++
		for(int k=0;k<tmp.rk;k++)
		{
			vecTmp1[k] = -co*rx_list[k];
			vecTmp2[k] = -si*rx_list[k];
		}		
		//++++++
 
		for(int j=0;j<tmp.rj;j++)
		{
			Vector3d Btemp(0,0,0);
 
			ray = j * tmp.de + tmp.R_de2;
 
			omy_ray_c = omy - ray * co;
			omz_ray_s = omz - ray * si;
			omy_ray_c2_omz_ray_s2 = omy_ray_c * omy_ray_c + omz_ray_s * omz_ray_s;
 
			si_r2_co_r1 = si * omz_ray_s + co * omy_ray_c;
 
			for(int k=0;k<tmp.rk;k++)
			{
				nr2 = rx2_list[k] + omy_ray_c2_omz_ray_s2;
				nr = sqrt(nr2);
 
				//if(nr>quantom_zero)
				{
					Btemp += Vector3d(si_r2_co_r1, vecTmp1[k], vecTmp2[k]) / (nr2*nr);
				}
			}
			bFirst = false;
 
            b+= (tmp.tmp1*j + tmp.tmp2)*Btemp;
        }
    }
    b*=Mu0_Pi4*tmp.J ;
    b=Mrot*b;
}

[Bousk]

Je te conseillerai déjà de déclarer toutes tes variables double (ray, omy_ray_c, ...) en-dehors des boucles. Comme ça tu ne les déclare qu'une seule fois et pas, dans le pire des cas, tmp^3.

Je suis pas sûr qu'il y ait un quelconque intérêt pour un type primitif.
C'est l'appel à un constructeur/destructeur qui est pénalisant.

[pepito3364]

ça va dépendre du compilateur et notamment des options de compilation. Mais j'ai déjà observé des gains de perf avec ce genre de petites optim.
Je suis pas expert, mais j'ai l'impression que le compilateur n'a pas besoin de remettre la variable sur un registre et vu le nombre de variable dans la fonction, on doit pouvoir utiliser un registre pour une variable. Certains compilateur doivent gérer ça correctement je suppose.
Il faudrait voir le code asm dans les 2 cas.

[Astraya]

A ce stade là, tu obtiendras plus de performance en travaillant sur l'organisation mémoire de tes données (bloc contigüe en mémoire, faire le moins de saut possible dans le bloc, surtout valable avec les boucles for dans une autre boucle for, fais un maximum de traitement sur des mémoires contigües)
Sinon un petit truc mais bon tu vas pas avoir de grande différence je pense;
Ce code la:

Code :

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

1
2
3
4
5
6
7
8
std::vector<double> rx_list(tmp.rk);
for(int k=0;k<tmp.rk;k++) 
    rx_list[k]=(-k*tmp.dl+omx_dl2_L2);
 
std::vector<double> rx2_list(tmp.rk);
 
for(int k=0;k<tmp.rk;k++) 
    rx2_list[k]=rx_list[k]*rx_list[k];

Deviens:

Code :

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

1
2
3
4
5
6
for(int k=0;k<tmp.rk;++k)
{
    double val = (-k*tmp.dl+omx_dl2_L2);
    rx_list[k]=val;
    rx2_list[k]=val*val;
}

Et enlève tout tes x++ et mets les en ++x. Certain compilateur sont assez intelligent pour optimiser la première version mais ne joue pas avec le feu

[gobelin88]

...

J'essayerai de calculer les indices 2 et 3 de Vector3d plus tot. Je sais pas si ce que je te propose est correct, j'ai fait vite fait.

Effectivement c'est une optimisation, merci bien vu!

J'ai déjà essayer de sortir la déclaration des variables des boucles mais cela ne change rien ...

Et enlève tout tes x++ et mets les en ++x. Certain compilateur sont assez intelligent pour optimiser la première version mais ne joue pas avec le feu

C'est vrai ... pris bonne note , mais cela n'a pas changé grand chose.

[pepito3364]

C'est déjà ça, sinon je voit pas comment gagner sur la racine..
Tu as pensé à utiliser la parallélisation? il faudrait faire attention à l'écriture sur "b" mais avec des mutex ça passe.
Bon courage pour la suite

[oodini]

Fais un reserve() sur Btemp.

[gobelin88]

C'est déjà ça, sinon je voit pas comment gagner sur la racine..
Tu as pensé à utiliser la parallélisation? il faudrait faire attention à l'écriture sur "b" mais avec des mutex ça passe.
Bon courage pour la suite

Oui tu as raison, mais ca complique pas mal le code ... j'ai pas trop envie de m'engagé sur cette voie sur cette fonction ... surtout que le programme est deja multithreadé .

[bretus]

A ce stade là, tu obtiendras plus de performance en travaillant sur l'organisation mémoire de tes données (bloc contigüe en mémoire, faire le moins de saut possible dans le bloc, surtout valable avec les boucles for dans une autre boucle for, fais un maximum de traitement sur des mémoires contigües)

A ce titre, fusionner ceci :

Code :

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

1
2
3
 
std::vector<double> vecTmp1(tmp);
std::vector<double> vecTmp2(tmp);

en

Code :

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

1
2
 
std::vector< std::pair< double, double > > vecTmp12 ;

ça aiderait en plein cœur des boucles ou ça ne sert à rien?

Code :

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

1
2
3
4
5
6
 
// avant vecTmp1[k] et vecTmp2[k] ne sont pas forcément proche
Btemp += Vector3d(si_r2_co_r1, vecTmp1[k], vecTmp2[k]) / (nr2*nr);
 
// après vecTmp1[k].first et vecTmp2[k].second se suivent en mémoire
Btemp += Vector3d(si_r2_co_r1, vecTmp1[k].first, vecTmp2[k].second) / (nr2*nr);

[Emmanuel Deloget]

Je ne sais pas quel est la précision recherchée ni l'intervalle sur lequel ces données s'étendent, mais si tu restes sur un espace assez réduit, alors tu peux certainement gagner en passant toutes tes valeurs flottantes en nombres à virgule fixe (int64, avec par exemple 48 bits pour la partie entière et 16 bits pour la partie décimale). Le gain peut être assez important.

Sinon, si la précision perdue ne te gène pas trop, tu peux peut être sortir une partie du calcul de ray dans la seconde boucle (j est un int, donc à chaque itération ray est augmenté de tmp.de+tmp.R_de2).

Pour le calcul de nr à partir de nr2, tu peux peut-être t'amuser avec la fonction fast invsqrt (http://en.wikipedia.org/wiki/Fast_inverse_square_root) dont il existe une version traitant des double quelque part. Tu obtiendra 1/sqrt(nr2), donc la valeur que tu cherches sera nr2*nr2*invsqrt(nr2). Le temps perdu sur la triple multiplication sera largement inférieur au temps gagné sur fast invsqrt (que tu feras inline, bien évidemment).

Et si tu a besoin de continuer à utiliser des flottant mais que la précision ne te gène pas trop, tu peux peut-être passer en float (tu gagneras pas mal en temps de traitement).

[imperio]

Pour les calculs parallelises, tu peux toujours utiliser la bibliotheque OpenMP. Je l'ai trouve assez facile a mettre en place, tu devrais peut-etre y jeter un coup d'oeil.

[gobelin88]

Merci pour vos propositions, alors je retiens:

Pour les calculs parallelises, tu peux toujours utiliser la bibliotheque OpenMP. Je l'ai trouve assez facile a mettre en place, tu devrais peut-etre y jeter un coup d'oeil.

Oui ca fait plusieurs fois que je regarde, j'ai encore pas fais le pas. Je vais essayer c'est une bonne occasion.

Pour le calcul de nr à partir de nr2, tu peux peut-être t'amuser avec la fonction fast invsqrt (http://en.wikipedia.org/wiki/Fast_inverse_square_root) dont il existe une version traitant des double quelque part. Tu obtiendra 1/sqrt(nr2), donc la valeur que tu cherches sera nr2*nr2*invsqrt(nr2). Le temps perdu sur la triple multiplication sera largement inférieur au temps gagné sur fast invsqrt (que tu feras inline, bien évidemment).

Magnifique lien! je regarde, reste a trouver/implémenter une version en double. De plus il suffit d'ajuster le nombre d’itérations à la précision voulue. je regarde ca ...

... Virgule fixe + Et si tu a besoin de continuer à utiliser des flottant mais que la précision ne te gène pas trop, tu peux peut-être passer en float (tu gagneras pas mal en temps de traitement).

Malheureusement je dois rester en double et la précision à une grande importance, car il s'agit d'une somme de petits termes donc l'erreur dérive vite.

std::vector< std::pair< double, double > > vecTmp12 ;

pas de gains ou très peu, merci

[gobelin88]

J'ai trouve ca :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
inline double invSqrt( const double& x )
{
    double y = x;
    double xhalf = ( double )0.5 * y;
    long long i = *( long long* )( &y );
    i = 0x5fe6ec85e7de30daLL - ( i >> 1 );//LL suffix for (long long) type for GCC
    y = *( double* )( &i );
    y = y * ( ( double )1.5 - xhalf * y * y );
 
    return y;
}

Voici un test performance precision :

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
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
 
inline double invSqrt_It1( const double& x )
{
    double y = x;
    double xhalf = ( double )0.5 * y;
    long long i = *( long long* )( &y );
    i = 0x5fe6ec85e7de30daLL - ( i >> 1 );//LL suffix for (long long) type for GCC
    y = *( double* )( &i );
    y = y * ( ( double )1.5 - xhalf * y * y );
    return y;
}
 
inline double invSqrt_It2( const double& x )
{
    double y = x;
    double xhalf = ( double )0.5 * y;
    long long i = *( long long* )( &y );
    i = 0x5fe6ec85e7de30daLL - ( i >> 1 );//LL suffix for (long long) type for GCC
    y = *( double* )( &i );
    y = y * ( ( double )1.5 - xhalf * y * y );
    y = y * ( ( double )1.5 - xhalf * y * y );
    return y;
}
 
inline double invSqrt_It3( const double& x )
{
    double y = x;
    double xhalf = ( double )0.5 * y;
    long long i = *( long long* )( &y );
    i = 0x5fe6ec85e7de30daLL - ( i >> 1 );//LL suffix for (long long) type for GCC
    y = *( double* )( &i );
    y = y * ( ( double )1.5 - xhalf * y * y );
    y = y * ( ( double )1.5 - xhalf * y * y );
    y = y * ( ( double )1.5 - xhalf * y * y );
    return y;
}
 
int main(int argc,char * argv[])
{
   QElapsedTimer timer;
   double x=0.0;
 
   timer.start();
   x=0.0;
   for(unsigned int i=0;i<10000000;i++){x+=1.0/sqrt(100.0);}
   std::cout<<"math="<<timer.elapsed()<<" sum="<<x<<std::endl;
 
   timer.start();
   x=0.0;
   for(unsigned int i=0;i<10000000;i++){x+=invSqrt_It1(100.0);}
   std::cout<<"math="<<timer.elapsed()<<" sum="<<x<<std::endl;
 
   timer.start();
   x=0.0;
   for(unsigned int i=0;i<10000000;i++){x+=invSqrt_It2(100.0);}
   std::cout<<"math="<<timer.elapsed()<<" sum="<<x<<std::endl;
 
   timer.start();
   x=0.0;
   for(unsigned int i=0;i<10000000;i++){x+=invSqrt_It3(100.0);}
   std::cout<<"math="<<timer.elapsed()<<" sum="<<x<<std::endl;
 
   return 0;
}

sur mon pc (core i5) ca donne :

math=8(ms) sum=1e+006
math=123(ms) sum=998430
math=176(ms) sum=999996
math=229(ms) sum=1e+006

bof

Donc cc semble pas rentable pour les double

[bretus]

Bonjour,

Merci d'avoir testé la pair.

Sinon, j'aurais tendance à faire voler au maximum les abstractions Vector3d au cœur de cette fonction.

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
 
//Vector3d Btemp(0,0,0);
double BtempX = 0.0 ;
double BtempY = 0.0 ;
double BtempZ = 0.0 ;
 
//Btemp += Vector3d(si_r2_co_r1, vecTmp1[k], vecTmp2[k]) / (nr2*nr);
double inv_nr2_nr = 1.0 / nr2 * nr ;
BtempX += si_r2_co_r1 * inv_nr2_nr ;
BtempY += vecTmp1[k] * inv_nr2_nr ;
BtempZ += vecTmp2[k] * inv_nr2_nr ;
 
//b+= (tmp.tmp1*j + tmp.tmp2)*Btemp;
double coef = tmp.tmp1*j + tmp.tmp2 ;
b[0] += coef * BtempX ;
b[1] += coef * BtempY ;
b[2] += coef * BtempZ ;

[gobelin88]

Bon je suis passé sous OpenMP et ca dechire! ...

Probleme resolu!