Discussion :

[Avatar36]

Bonjour,

Je requiers votre aide aujourd'hui car j'ai réalisé un algorithme pour résoudre un système d'équation basé sur le pivot de Gauss. Mon problème est qu'il ralentit arrivé à 40% (voir l'émission du signal "progressChanged()" pour l'avancement). Ma matrice "matCopy" fait environ 1600 lignes et 1600 colonnes. Donc certes c'est long mais pourquoi ça ralentit à 40%??? Entre 0 et 40%, on avance d' 1% par secondes mais à 40% on est à 1% toutes les 5-10 sec. Sur quoi est-ce que je peux jouer pour l'accélérer un peu?

Voici le code:

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
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
MyThread::results MyThread::solve(QVector<QVariant>& boundariesCdt, QVector<double>& forces)
{
    results result; //Variable qui va contenir les résultats
 
    //On copie la matrice de rigidité dans une autre matrice que l'on pourrat modifier
    matCopy.setSize(mat.rowCount(),mat.columnCount());
    for(int i=0;i<mat.rowCount();i++)
    {
        for(int j=0;j<mat.columnCount();j++)
        {
            matCopy(i,j) = mat(i,j);
        }
    }
 
    //On copie les forces dans resultat.displacements (sera transformé par le pivot de gauss pour la résolution)
    for(int j=0;j<forces.size();j++)
    {
        result.displacements.push_back(forces.at(j));
    }
 
    //Comme certaines réactions aux connecteurs sont inconnues mais on connait les déplacement, on ne peut pas appliquer le pivot de
    //Gauss à certaines lignes. La seule solution est de la résoudre à part et d'appliquer le pivot aux lignes où le membre de droite
    //est connu. Une fois qu'on aura calculé les déplacements, il suffira de résoudre les équations mises à part
    //pour les efforts aux connecteurs.
 
    for(int j=0;j<boundariesCdt.size();j++)
    {
        if(boundariesCdt.at(j)!=QVariant()) //Si le déplacement (ou la rotation) vaut une valeur, on modifie la matrice.
            //Sinon on ne fait rien.
        {
            if(j<6) //Si on est dans les 6 premières lignes de la matrice (le noeud 1)
            {
                for(int i=0; i<matCopy.columnCount();i++)
                {
                    if(i==j) //si i=j alors on incrit un 1, sinon un 0. Permet de faire une équation xi*1=0
                    {
                        matCopy(j,i) = 1;
                    }
                    else
                    {
                        matCopy(j,i) = 0;
                    }
                }
                result.displacements[j] = boundariesCdt[j].toDouble();
            }
            else //Si on est dans les 6 dernières lignes de la matrice (le dernier noeud)
            {
                for(int i=0; i<matCopy.columnCount();i++)
                {
                    if(i==matCopy.rowCount()-12+j)//si i=n-12+j alors on incrit un 1, sinon un 0. Permet de faire une équation xi*1=0
                    {
                        matCopy(matCopy.rowCount()-12+j,i) = 1;
                    }
                    else
                    {
                        matCopy(matCopy.rowCount()-12+j,i) = 0;
                    }
                }
                result.displacements[result.displacements.size()-12+j] = boundariesCdt[j].toDouble();
            }
        }
    }
 
    //###############################################################################################################################
    //Pivot de Gauss appliqué à "matCopy" et "result.displacement"
    int r = 0;
    int k = 0;
    double maxPivot = 0;
 
    for(double j=0; j<matCopy.columnCount();j++)
    {
        emit progressChanged(floor(100*j/(matCopy.columnCount()-1)));
        maxPivot = 0; //On met "maxPivot" à une valeur minimale
 
        //On cherche le coef max de la colonne j et on le stocke dans "maxPivot". L'indice de la ligne du pivot est k.
        for(int i=r;i<matCopy.rowCount();i++)
        {
            if(matCopy(i,j)*matCopy(i,j)>maxPivot*maxPivot)
            {
                maxPivot = matCopy(i,j);
                k = i;
            }
        }
 
        if(matCopy(k,j)!=0)//Si le pivot est non nul, on procède aux opérations. Sinon, cela veut dire que toute la colonne est nulle
        {
            //############ On divise la ligne par le pivot (1 en début de ligne) ##############################
            for(int i = 0; i<matCopy.columnCount();i++)
            {
                matCopy(k,i) /= maxPivot;
            }
            result.displacements[k] /= maxPivot;
 
            double  coef = 0;
 
            //############ Echange des lignes r et k ##############################
            for(int i = 0; i<matCopy.columnCount();i++)
            {
                coef = matCopy(k,i);
                matCopy(k,i) = matCopy(r,i); //On modifie la ligne k
                matCopy(r,i) = coef; //On modifie la ligne r
            }
            coef = forces[r];
            result.displacements[r] = result.displacements[k];
            result.displacements[k] = coef;
 
            for(int i=0; i<matCopy.rowCount();i++)
            {
                if(i!=r)
                {
                    double coefLigne = matCopy(i,j);
 
                    for(int l=0;l<matCopy.columnCount();l++)
                    {
                        matCopy(i,l) -= matCopy(r,l)*coefLigne; // [Li] = [Li] - [Lr]*coefLigne
                    }
 
                    result.displacements[i] -= result.displacements[r]*coefLigne;
                }
            }
            r++;
        }
    }
    //###############################################################################################################################
    //Maintenant on va calculer les efforts aux connecteurs
    for(int i=0;i<mat.rowCount();i++)
    {
        result.reactionLoads.push_back(0);
 
        for(int j=0;j<mat.columnCount();j++)
        {
            result.reactionLoads[i] += mat(i,j)*result.displacements[j];
        }
    }
 
    //"forces"              est un tableau contenant toutes les solutions des déplacements
    //"connectorReaction"   est un tableau contenant toutes les solutions des efforts
 
    return result;
}

Cordialement

[Iradrille]

Hello,

Dès qu'on parles optimisation il n'y à qu'une réponse valide : profile et vois où ça bloque.

Quelques pistes à creuser

Code :

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

1
2
3
4
5
6
7
8
9
//On copie la matrice de rigidité dans une autre matrice que l'on pourrat modifier
matCopy.setSize(mat.rowCount(),mat.columnCount());
for(int i=0;i<mat.rowCount();i++)
{
	for(int j=0;j<mat.columnCount();j++)
	{
		matCopy(i,j) = mat(i,j);
	}
}

Utilise plutôt un ctor par copie qui pourra copier la matrice efficacement
Via std::vector::assign ou std::copy (ou équivalent Qt)

Code :

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

1
2
3
4
5
//On copie les forces dans resultat.displacements (sera transformé par le pivot de gauss pour la résolution)
for(int j=0;j<forces.size();j++)
{
	result.displacements.push_back(forces.at(j));
}

Pourquoi pas un simple result.displacements.assign(forces.begin(), forces.end()); ? (Au minimum reserve pour éviter les multiples réallocations)

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
//Comme certaines réactions aux connecteurs sont inconnues mais on connait les déplacement, on ne peut pas appliquer le pivot de
//Gauss à certaines lignes. La seule solution est de la résoudre à part et d'appliquer le pivot aux lignes où le membre de droite
//est connu. Une fois qu'on aura calculé les déplacements, il suffira de résoudre les équations mises à part
//pour les efforts aux connecteurs.
 
for(int j=0;j<boundariesCdt.size();j++)
{
	if(boundariesCdt.at(j)!=QVariant()) //Si le déplacement (ou la rotation) vaut une valeur, on modifie la matrice.
		//Sinon on ne fait rien.
	{
		if(j<6) //Si on est dans les 6 premières lignes de la matrice (le noeud 1)
		{
			for(int i=0; i<matCopy.columnCount();i++)
			{
				if(i==j) //si i=j alors on incrit un 1, sinon un 0. Permet de faire une équation xi*1=0
				{
					matCopy(j,i) = 1;
				}
				else
				{
					matCopy(j,i) = 0;
				}
			}
			result.displacements[j] = boundariesCdt[j].toDouble();
		}
		else //Si on est dans les 6 dernières lignes de la matrice (le dernier noeud)
		{
			for(int i=0; i<matCopy.columnCount();i++)
			{
				if(i==matCopy.rowCount()-12+j)//si i=n-12+j alors on incrit un 1, sinon un 0. Permet de faire une équation xi*1=0
				{
					matCopy(matCopy.rowCount()-12+j,i) = 1;
				}
				else
				{
					matCopy(matCopy.rowCount()-12+j,i) = 0;
				}
			}
			result.displacements[result.displacements.size()-12+j] = boundariesCdt[j].toDouble();
		}
	}
}

généralement il vaut mieux avoir

Code :

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

1
2
3
4
5
6
7
8
9
10
if(...) {
   for(...);
else {
   for(...);
}
// Que
for(...) {
   if(...);
   else;
}

Le if interne peut être remplacé par un opérateur ternaire (pour aider le compilo à générer un cmov), voir même remplacé par

Code :

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

1
2
3
4
5
6
7
8
9
10
11
if(i==matCopy.rowCount()-12+j)//si i=n-12+j alors on incrit un 1, sinon un 0. Permet de faire une équation xi*1=0
{
	matCopy(matCopy.rowCount()-12+j,i) = 1;
}
else
{
	matCopy(matCopy.rowCount()-12+j,i) = 0;
}
// devient
auto const jj = matCopy.rowCount()-12+j;
matCopy(jj,i) = 1 - ((jj - i) & 1);

Quand tu travailles sur un membre d'une classe, copie le localement et travailles sur la copie si possible, le compilo sera heureux.
Exemple

Code :

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

1
2
3
4
5
6
7
8
9
//On cherche le coef max de la colonne j et on le stocke dans "maxPivot". L'indice de la ligne du pivot est k.
for(int i=r;i<matCopy.rowCount();i++)
{
	if(matCopy(i,j)*matCopy(i,j)>maxPivot*maxPivot)
	{
		maxPivot = matCopy(i,j);
		k = i;
	}
}

Beaucoup d'appels de fonctions dans cette boucle (le compilo va très certainement les inliner, mais ça représente quand même du calcul inutile.
matCopy(i,j) -> Est probablement du genre matCopy.m_data[j * m_width + i]. Le calcul de l'index tu n'as pas besoin de le faire plusieurs fois.

Le fait de travailler sur des variables locales autorise aussi le compilo à faire des optimisations plus aggressives.

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
auto const rCount = matCopy.rowCount();
auto maxPivotSq = maxPivot;
 
for(int i=r;i<rCount;i++)
{
	auto val = matCopy(i, j);
	if(val * val > maxPivotSq)
	{
		maxPivot = val;
		maxPivotSq = val * val;
		k = i;
	}
}

Il est aussi peut être intéressant de fournir un pointeur sur le contenu de tes matrices, pour itérer dessus simplement sans avoir à recalculer l'index à chaque fois.

Vérifie aussi que tes matrices soient "dans le bon sens".

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
struct Mat4a {
   double data[16];
   double& operator()(unsigned x, unsigned y) { return data[y * 4 + x]; }
};
 
struct Mat4b {
   double data[16];
   double& operator()(unsigned x, unsigned y) { return data[x * 4 + y]; }
};
 
Mat4a ma;
Mat4b mb;
for(unsigned x=0; x<4; ++x) {
   for(unsigned y=0; y<4; ++y) {
      ma(x, y) = 42.0; // lent
      mb(x, y) = 42.0; // rapide
   }
}
 
for(unsigned y=0; y<4; ++y) {
   for(unsigned x=0; x<4; ++x) {
      ma(x, y) = 42.0; // rapide
      mb(x, y) = 42.0; // lent
   }
}

Bref, profile et vois où ça bloque pour améliorer ça.

[Avatar36]

Les optimisation que tu m'as proposées, je les ai mises en place mais ça n'a pas amélioré beaucoup le process étant donné que ce ne sont pas les endroits qui peinent. En revanche ça rend le code plus clair et plus simple. Merci pour cette réponse.
Après avoir fait 2-3 manip, j'en suis arrivé à la conclusion que la partie "lente" de l'algorithme provient de cette boucle:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
            for(int i=0; i<rCount;i++)
            {
                if(i!=r)
                {
                    double coefLigne = matCopy(i,j);
 
                    for(int l=0;l<cCount;l++)
                    {
                        matCopy(i,l) -= matCopy(r,l)*coefLigne; // [Li] = [Li] - [Lr]*coefLigne
                    }
 
                    result.displacements[i] -= result.displacements[r]*coefLigne;
                }
            }

C'est en fait la dernière étape de l'algorithme de Gauss. Le problème c'est que je vois pas quoi faire de plus ici.

Je ne pense pas pouvoir faire mieux étant donné que j'opère sur toutes les cases de la matrice sauf sur la ligne r. Donc matrice 1600*1600-1600=2558400 opérations. Et je pense pas pouvoir raccourcir ça. Sauf si quelqu'un a une idée .

[Avatar36]

J'ai trouvé un petit trick pour raccourcir les calculs nécessaires dans cette boucle!

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
            for(int i=0; i<rCount;i++)
            {
                if(i!=r)
                {
                    double coefLigne = matCopy(i,j);
 
                    for(int l=r;l<cCount;l++)
                    {
                        matCopy(i,l) -= matCopy(r,l)*coefLigne; // [Li] = [Li] - [Lr]*coefLigne
                    }
 
                    result.displacements[i] -= result.displacements[r]*coefLigne;
                }
            }

Regardez bien la boucle avec l'incrément l. L'indice démarre à r. En effet, de la colonne l à r-1, les coefficients de la matrice valent 0 sur la ligne du pivot. Par cette méhode je retire un paquets d'opérations de soustraction de 0 à des coefficients. Maintenant ça s'execute nettement plus vite qu'avant mais toujours pas assez.
Avez-vous d'autres idées?

[darkman19320]

Salut,

Quel est le type de conteneur que tu utilises pour matCopy?

[jo_link_noir]

Je me demande s'il ne serait pas plus rapide de faire une boucle pour [i,r[ et une autre pour ]r, rCount] ?
Comme la dit Iradrille, essaye de modifier le calcul de l'indice pour ne plus avoir matCopy(i,l) mais matCopyLine[i], je ne suis pas sûr que dans l'état actuel le compilateur fasse bien son boulot d'optimisation.
Peut-être aussi mettre le calcul sur result.displacements dans une boucle à part.

[Iradrille]

Donc certes c'est long mais pourquoi ça ralentit à 40%??? Entre 0 et 40%, on avance d' 1% par secondes mais à 40% on est à 1% toutes les 5-10 sec.

Bon j'étais curieux, j'ai implémenté les morceaux manquants pour tester, j'ai pas observé de ralentissements et j'ai pas réussi à reproduire ça. Peut être une histoire de nombres dénormalisés.

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
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
#include <vector>
#include <iostream>
#include <memory>
#include <random>
#include <algorithm>
#include <chrono>
#include <string>
 
template <class Function>
void timeIt(Function function, std::string const& desc) {
	auto startTime = std::chrono::high_resolution_clock::now();
	function();	
	auto endTime = std::chrono::high_resolution_clock::now();
	auto t = std::chrono::duration_cast<std::chrono::duration<double>>(endTime - startTime);
	std::cout << desc << " " << t.count() << " seconds." << std::endl;
}
 
struct Result {
	std::vector<double> displacements;
	std::vector<double> reactionLoads;
};
 
struct Mat {
 
	Mat(int x, int y): m_width(x), m_height(y), m_data(x*y) { }
	Mat(Mat const& other): m_width(other.m_width), m_height(other.m_height), m_data(other.m_data) { }
 
	std::vector<double> m_data;
 
	int m_width;
	int m_height;
 
	int rowCount() const { return m_height; }
	int columnCount() const { return m_width; }
 
	double& operator()(int x, int y) { return m_data[y * m_width + x]; }
};
 
Result solve(Mat &mat, std::vector<double>& boundariesCdt, std::vector<double>& forces) {
	Result result; //Variable qui va contenir les résultats
 
	Mat matCopy = mat;
 
	using std::begin;
	using std::end;
 
	result.displacements.resize(forces.size());
	std::copy(begin(forces), end(forces), std::begin(result.displacements));
 
	//Comme certaines réactions aux connecteurs sont inconnues mais on connait les déplacement, on ne peut pas appliquer le pivot de
	//Gauss à certaines lignes. La seule solution est de la résoudre à part et d'appliquer le pivot aux lignes où le membre de droite
	//est connu. Une fois qu'on aura calculé les déplacements, il suffira de résoudre les équations mises à part
	//pour les efforts aux connecteurs.
 
	for(int j=0;j<boundariesCdt.size();j++) {
		//if(boundariesCdt.at(j)!=QVariant()) { //Si le déplacement (ou la rotation) vaut une valeur, on modifie la matrice.
			//Sinon on ne fait rien.
			if(j<6) { //Si on est dans les 6 premières lignes de la matrice (le noeud 1)
				for(int i=0; i<matCopy.columnCount();i++)
				{
					if(i==j) //si i=j alors on incrit un 1, sinon un 0. Permet de faire une équation xi*1=0
					{
						matCopy(j,i) = 1;
					}
					else
					{
						matCopy(j,i) = 0;
					}
				}
				//result.displacements[j] = boundariesCdt[j].toDouble();
				result.displacements[j] = boundariesCdt[j];
			}
			else { //Si on est dans les 6 dernières lignes de la matrice (le dernier noeud)
				for(int i=0; i<matCopy.columnCount();i++) {
					if(i==matCopy.rowCount()-12+j) { //si i=n-12+j alors on incrit un 1, sinon un 0. Permet de faire une équation xi*1=0
						matCopy(matCopy.rowCount()-12+j,i) = 1;
					}
					else {
						matCopy(matCopy.rowCount()-12+j,i) = 0;
					}
				}
				result.displacements[result.displacements.size()-12+j] = boundariesCdt[j];
			}
		//}
	}
 
	//###############################################################################################################################
	//Pivot de Gauss appliqué à "matCopy" et "result.displacement"
	int r = 0;
	int k = 0;
	double maxPivot = 0;
 
	for(double j=0; j<matCopy.columnCount(); j++) {
		//std::cout << floor(100*j/(matCopy.columnCount()-1)) << std::endl;
 
		maxPivot = 0; //On met "maxPivot" à une valeur minimale
 
		//On cherche le coef max de la colonne j et on le stocke dans "maxPivot". L'indice de la ligne du pivot est k.
		for(int i=r;i<matCopy.rowCount();i++)
		{
			if(matCopy(int(i), int(j))*matCopy(int(i), int(j)) > maxPivot*maxPivot) {
				maxPivot = matCopy(int(i), int(j));
				k = i;
			}
		}
 
		if(matCopy(int(k), int(j))!=0)//Si le pivot est non nul, on procède aux opérations. Sinon, cela veut dire que toute la colonne est nulle
		{
			//############ On divise la ligne par le pivot (1 en début de ligne) ##############################
			for(int i = 0; i<matCopy.columnCount();i++)
			{
				matCopy(k,i) /= maxPivot;
			}
			result.displacements[k] /= maxPivot;
 
			double  coef = 0;
 
			//############ Echange des lignes r et k ##############################
			for(int i = 0; i<matCopy.columnCount();i++)
			{
				coef = matCopy(k,i);
				matCopy(k,i) = matCopy(r,i); //On modifie la ligne k
				matCopy(r,i) = coef; //On modifie la ligne r
			}
			coef = forces[r];
			result.displacements[r] = result.displacements[k];
			result.displacements[k] = coef;
 
			for(int i=0; i<matCopy.rowCount();i++)
			{
				if(i!=r)
				{
					double coefLigne = matCopy(int(i), int(j));
 
					for(int l=0;l<matCopy.columnCount();l++)
					{
						matCopy(i,l) -= matCopy(r,l)*coefLigne; // [Li] = [Li] - [Lr]*coefLigne
					}
 
					result.displacements[i] -= result.displacements[r]*coefLigne;
				}
			}
			r++;
		}
	}
	//###############################################################################################################################
	//Maintenant on va calculer les efforts aux connecteurs
	for(int i=0;i<mat.rowCount();i++)
	{
		result.reactionLoads.push_back(0);
 
		for(int j=0;j<mat.columnCount();j++)
		{
			result.reactionLoads[i] += mat(i,j)*result.displacements[j];
		}
	}
 
	//"forces"              est un tableau contenant toutes les solutions des déplacements
	//"connectorReaction"   est un tableau contenant toutes les solutions des efforts
 
	return std::move(result);
}
 
int main(int, char**) {
 
	using std::begin;
	using std::end;
 
	Mat m(1600, 1600);
	std::vector<double> bounds(12); // size ?
	std::vector<double> forces(m.rowCount()); // size ?
 
	std::default_random_engine generator;
	std::uniform_real_distribution<double> distribution(-42.0, 42.0);
	auto gen = [&]() { return distribution(generator); };
 
	std::generate(begin(m.m_data), end(m.m_data), gen);
	std::generate(begin(bounds), end(bounds), gen);
	std::generate(begin(forces), end(forces), gen);
 
	timeIt([&]() {
		auto ret = solve(m, bounds, forces);
		std::cout << ret.displacements[42] << std::endl;
	}, "solve");
 
	return 0;
}

Bref, après un rapide profiling, le point qui pose problème est évident (temps d'exec : 51.4089 sec, donc j'imagine que j'ai implémenté le reste de la même façon que toi).



97.7% du temps total passé sur cette ligne, rien que ça.
Cette ligne représente ~4 milliards (1600^3) multiplications / soustractions de doubles
-> ça devrait représenter quelques secondes maximum et pas quasiment une minute
-> temps d'exec plombé par un nombre inimaginable de cache misses.

En changeant simplement la façon d'accéder aux données dans ta matrice tu peux réduire drastiquement le nombre de cache misses.

Code :

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

1
2
3
4
struct Mat {
	// ...
	double& operator()(int x, int y) { return m_data[x * m_height + y]; }
};

Résultat

Code :

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

solve 3.71921 seconds.

[Avatar36]

Pour ce qui est de l'accès à la matrice, j'ai ça:

Code :

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

1
2
    double operator()(int rowId, int columnId) const;
    double& operator()(int rowId, int columnId);

Code :

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

1
2
3
4
5
6
7
8
9
double Matrix::operator()(int rowId, int columnId) const
{
    return data[column*rowId+columnId];
}
 
double& Matrix::operator()(int rowId, int columnId)
{
    return data[column*rowId+columnId];
}

Ce qui est la même chose que toi. Le conteneur est un QVector. Qu'est-ce que je dois changer du coup?

@jo_link_noir: j'ai déjà essayé en retirant le if associé du coup => aucune amélioration notable.

[Iradrille]

Ce qui est la même chose que toi. Le conteneur est un QVector. Qu'est-ce que je dois changer du coup?

Hum, étonnant de voir si grande différence de perf du coup.

Question conne mais tu compile bien en release (ou -O3 avec GCC / MinGW / clang) ? L'utilisation de la STL (et c'est probablement pareil pour Qt) est extrêmement lente en debug (plein de tests "inutiles" qui rendent le debug plus simple mais qui plombent le temps d'exec).

Est-ce que tu as essayé de remplacer l'operator[] par QVector::at() ? D'après la doc operator[] peut provoquer une deep copy (le COW me dépasse "un peu", donc j'ai aucune idée de si tu te trouves dans un cas qui va provoquer ça).

Note that using non-const operators can cause QVector to do a deep copy.

De même tes getters sont biens inline ? Essaie de forcer l'inlining pour voir.

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
#if defined(__GNUC__) || defined(__MINGW32__) // GCC et MinGW
#    define FORCE_INLINE __attribute__((always_inline))
#elif defined(_MSC_VER) // MSVC
#    define FORCE_INLINE __forceinline
#else // autres
#    define FORCE_INLINE inline
#endif // FORCE_INLINE
 
FORCE_INLINE double operator()(int rowId, int columnId) const {
    return data.at(column*rowId+columnId);
}
 
FORCE_INLINE double& operator()(int rowId, int columnId) {
    return data.at(column*rowId+columnId);
}

(M'enfin, ya pas de raison que ces fonctions ne soient pas inline, peu de chances que ce soit ça...).

Tu as profiler le programme pour voir ce qui bloque ?
Je t'avoue être de plus en plus curieux.

[Avatar36]

Je compile bien en release et j'ai assayé le IN_LINE: sans améliorations.
Comment est-ce que tu fais le profiling du code?
J'ai utilisé QElapsedTimer mais il n'est pas assez précis.
J'arrive néanmoins à déterminer que c'est bien ce morceau qui ralentit le processus:

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
            for(int i=0; i<rCount;i++)
            {
 
                if(i!=r)
                {
                    double coefLigne = matCopy(i,j);
 
 
 
                    for(int l=r;l<cCount;l++)
                    {
                        matCopy(i,l) -= matCopy(r,l)*coefLigne; // [Li] = [Li] - [Lr]*coefLigne
                    }
 
                    result.displacements[i] -= result.displacements[r]*coefLigne;
                }
 
            }

Par contre je ne sais pas comment faire pour déterminer quelle ligne est la plus lente.
Ce qui est bizarre, c'est que d'un seul coup, je passe de 14 millisecondes à 174 ms sur ce bout de code.
Est-ce que ça peut être lié au contenu de la matrice? Par exemple au début on divise des nombres entiers puis au fur et à mesure les nombres décimaux sont de plus en plus longs (au niveau de la précision).
Je vais tester en changeant le contenu de la matrice (on sait jamais).

[ternel]

Un profiler est un programme qui observe l'exécution de ton code.
Un assez populaire est gprofile, mais il est très lié à linux.

Au besoin, si tu cherche gprofile sur le forum, tu trouveras le nom de celui habituellement conseillé pour windows.

Après exécution de ton programme, il te fera un rapport te disant énormément de chose, dont les fonctions les plus utilisées, et les plus gourmandes.

[Avatar36]

Je viens d'essayer en changeant le contenu et effectivement c'est allé jusqu'au bout en 5 secondes. J'ai utilisé ça pour remplir la matrice:

Code :

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

1
2
3
4
5
6
7
    for(int i=0;i<mat.rowCount();i++)
    {
        for(int j=0;j<mat.columnCount();j++)
        {
            mat(i,j) = i+j;
        }
    }

Avant j'avais ça:

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
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
    double Iy = M_PI*(pow(p_pipe->OD,4)-pow(p_pipe->OD-2*p_pipe->WT,4))/64; //Moment quadratique
    double Iz = M_PI*(pow(p_pipe->OD,4)-pow(p_pipe->OD-2*p_pipe->WT,4))/64; //Moment quadratique
    double J  = M_PI*(pow(p_pipe->OD,4)-pow(p_pipe->OD-2*p_pipe->WT,4))/32; //Moment quadratique polaire
    double S  = M_PI*(pow(p_pipe->OD,2)-pow(p_pipe->OD-2*p_pipe->WT,2))/4; //Section du pipe
    double G  = p_pipe->E/(2*(1+p_pipe->poissonRatio)); //Module de cisaillement
    double E  = p_pipe->E; //Module d'Young
 
    meshPipe(); //Remplissage de "nodes" (création des noeuds)
 
    //###################################### CONSTRUCTION DE LA MATRICE DE RAIDEUR ###################################################
 
    mat.setSize(6*nodes.size(),6*nodes.size());
 
    Matrix *matElement = new Matrix; //On définit une matrice d'un élément (12x12) (adaptée pour chaque élément)
    matElement->setSize(12,12);
 
    double L = 0; //Longueur de l'élément (adaptée pour chaque élément)
 
    for(int k=0;k<nodes.size()-1;k++)//On itère sur les éléments (nombre de noeuds-1)
    {
        L = nodes.at(k).data.distanceToPoint(nodes.at(k+1).data); //Longueur de l'élément
 
matElement->resetToZero(); //Remplit matElement de 0
 
        //Remplissage du triangle supérieur de la matrice
        (*matElement)(0,0) = E*S/L;
        (*matElement)(0,6) = -E*S/L;
        (*matElement)(1,1) = 12*E*Iy/pow(L,3);
        (*matElement)(1,5) = 6*E*Iz/pow(L,2);
        (*matElement)(1,7) = -12*E*Iz/pow(L,3);
        (*matElement)(1,11) = 6*E*Iz/pow(L,2);
        (*matElement)(2,2) = 12*E*Iy/pow(L,3);
        (*matElement)(2,4) = -6*E*Iy/pow(L,2);
        (*matElement)(2,8) = -12*E*Iy/pow(L,3);
        (*matElement)(2,10) = -6*E*Iy/pow(L,2);
        (*matElement)(3,3) = G*J/L;
        (*matElement)(3,9) = -G*J/L;
        (*matElement)(4,4) = 4*E*Iy/L;
        (*matElement)(4,8) = 6*E*Iy/pow(L,2);
        (*matElement)(4,10) = 2*E*Iy/L;
        (*matElement)(5,5) = 4*E*Iz/L;
        (*matElement)(5,7) = -6*E*Iz/pow(L,2);
        (*matElement)(5,11) = 2*E*Iz/L;
        (*matElement)(6,6) = E*S/L;
        (*matElement)(7,7) = 12*E*Iz/pow(L,3);
        (*matElement)(7,11) = -6*E*Iz/pow(L,2);
        (*matElement)(8,8) = 12*E*Iy/pow(L,3);
        (*matElement)(8,10) = 6*E*Iy/pow(L,2);
        (*matElement)(9,9) = G*J/L;
        (*matElement)(10,10) = 4*E*Iy/L;
        (*matElement)(11,11) = 4*E*Iz/L;
 
        //Remplissage de la moitié inférieure par symétrie
        for(int i=0;i<matElement->rowCount();i++)
        {
            for(int j=0;j<i;j++)
            {
                (*matElement)(i,j) = (*matElement)(j,i);
            }
        }
        //On crée une matrice rotation pour mettre la matrice de rigidité de l'élément dans le repère global.
        //On utilisera [K]=[R]t[K'][R] pour changer de repère (avec [R]t = [R]-1, transposée egal à l'inverse car transformation
        //orthogonale).
 
        Matrix rotation;
        rotation.setSize(12,12);
        rotation(0,0) = (nodes.at(k+1).data.x-nodes.at(k).data.x)/nodes.at(k).data.distanceToPoint(nodes.at(k+1).data);
        rotation(0,1) = (nodes.at(k+1).data.y-nodes.at(k).data.y)/nodes.at(k).data.distanceToPoint(nodes.at(k+1).data);
        rotation(0,2) = (nodes.at(k+1).data.z-nodes.at(k).data.z)/nodes.at(k).data.distanceToPoint(nodes.at(k+1).data);
        rotation(1,0) = -rotation(0,1)/(pow(pow(rotation(0,1),2)+pow(rotation(0,0),2),0.5));
        rotation(1,1) = rotation(0,0)/(pow(pow(rotation(0,1),2)+pow(rotation(0,0),2),0.5));
        rotation(1,2) = 0;
        rotation(2,0) = -rotation(0,2)*rotation(1,1);
        rotation(2,1) = -rotation(0,2)*rotation(1,0);
        rotation(2,2) = rotation(0,0)*rotation(1,1)-rotation(1,0)*rotation(0,1);
 
        //On remplit le reste de la matrice pour obtenir une matrice de rotation 12x12 (taille de la matrice rigidité)
        for(int l=0;l<4;l++)
        {
            for(int j=0;j<3;j++)
            {
                for(int i=0 ; i<3 ; i++)
                {
                    rotation(l*3+i,l*3+j) = rotation(i,j);
                }
            }
        }
 
        Matrix rotationTranspose;
        rotationTranspose.setSize(12,12);
 
        for(int i=0;i<12;i++)
        {
            for(int j=0;j<12;j++)
            {
                rotationTranspose(i,j) = rotation(i,j);
            }
        }
        rotationTranspose.transpose();
 
        //Changement de repère de la matrice grâce à la matrice de rotation
        (*matElement) = rotationTranspose*(*matElement)*rotation;
 
        //Assemblage avec la matrice globale de rigidité
        for(int i=0;i<matElement->rowCount();i++)
        {
            for(int j=0;j<matElement->columnCount();j++)
            {
                mat(k*6+i,k*6+j) += (*matElement)(i,j);
            }
        }
    }

Est-ce qu'il y a quelque chose que je ne dois pas faire?

[Avatar36]

Ah oui effectivement, à partir d'un certain nombre de matrices, les 0 ne sont plus des 0 mais des nan.
Je vais essayer de voir d'où ça vient.

EDIT: En fait je modifiais les cases de la matrice sans jamais les remettre à 0. Du coup quand je réutilisais matElement pour l'élément suivant, j'utilisais en partie les valeurs de la matrice précédente. J'ai ajouté une fonction resetToZero à ma classe Matrix pour remettre la valeur des cases à 0. Bon par contre ça ralentit encore à 42%...

EDIT²: Il y a encore des nan dans la matrice. Je vais tâcher de savoir d'où ils viennent.

[ternel]

Essaie une taille plus petite (genre 50 par 50), et vérifie que les calculs sont bons.
Si tes formules ne sont pas exactes, ca ne sert à rien d'optimiser, il faudra changer le code (et donc probablement détruire les optimisations).

Vite mais bien.

[Avatar36]

C'est bon, j'ai trouvé le problème!
C'était dans la déclaration de la matrice "rotation", il pouvait survenir des cas où je divisais par 0. J'ai géré ce cas avec un boucle if. Maintenant il n'y a plus de ralentissements à 42%.
Bon même si c'était encore un problème pourri (comme 90% de mes problèmes), ça m'a au moins permis d'améliorer plusieurs trucs dans cet algorithme. Maintenant il s’exécute en une petite dizaine de secondes (parfait pour moi).
Merci de votre aide.

Cordialement

[Bousk]

Ah oui effectivement, à partir d'un certain nombre de matrices, les 0 ne sont plus des 0 mais des nan.

Tu dois avoir une division par un float ou double qui vaut 0 quelque part pour avoir un truc pareil. A partir de là, l'escalade de la violence peut débuter.