Discussion :

[un passant]

j'ai programme un code qui me resoud une equation de type a*x=b avec a etant une matrice carree et x et b des vecteurs.
pour ce faire j'utilise une resolution par la methode du pivot de gauss.
au depart je faisai avec des exemple de petite matrice (200*200)
mais maintenant je dois resoudre des systemes avec des matrices de plus de 1500*1500. le resulta est toujours bon, le seul probleme est que le temps de calcul est eleve et jaimerai savoir si il y a des methodes pour l'optimiser.

[contremaitre]

Tu es dans le forum C, donc il faudrait montrer ton code pour voir ce qu'ily a à améliorer...
Mais tu aurais peut etre plus de chance dans le forum algorithme

[un passant]

j'ai mis l'essentiel du code;
je savais pas qu'il y avait un forun algorithme.

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
main()
{
float **K, *q, *f;
int nnoeuds;//lu ds le fichier d'entree
 
//Allocation de memoire
    FILE * fp;
    fp=fopen("erreur.txt","w");
    K = (float **)calloc ( nnode , sizeof(float*) );
    if( K == NULL )
    {
        fprintf(fp,"Allocation de K impossible");
        exit(EXIT_FAILURE);
    }
    for(int j=0 ;j<nnode; j++ )
    {
        K[j] =(float *)calloc(nnode, sizeof(float));
        if( K[j] == NULL )
        { 
            fprintf(fp,"Allocation de a impossible");
            exit(EXIT_FAILURE);
        }
    }
    q =(float *)calloc(nnode, sizeof(float));
    if( q == NULL )
    {
        fprintf(fp,"Allocation de q impossible");
        exit(EXIT_FAILURE);
    }
    f =(float *)calloc(nnode, sizeof(float));
    if( f == NULL )
    {
        fprintf(fp,"Allocation de f impossible");
        exit(EXIT_FAILURE);
    }
 
//creation de la matrice de rigidite et du vecteur chargement
    matrig(K);
    vectfor(f);
 
//Resolution de K*Q=F
    gauss(K, q, f, nnoeuds);
    return 0;
}
 
void gauss(float **a,float *x,float *b,int n)
{
	float p,s;
    int i,j,k;
	for(k=0;k<n-1;k++)
    {
		if (a[k][k]==0)
        {
		    printf("\n\n * Un pivot nul ! => methode de Gauss non applicable\n\n");
		    break;
		}
		for(i=k+1;i<n;i++)  //triangularisation de la matrice
        {
			p=a[i][k]/a[k][k];
			for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
			b[i]=b[i]-p*b[k];
		}
	}
	for(i=n-1;i>=0;i--) //Résolution du systeme
    {
		s=0;
		for(j=i+1;j<n;j++)s=s+a[i][j]*x[j];
		x[i]=(b[i]-s)/a[i][i];
	}
}

[je®ome]

Un petit truc en regardant vite fait. Il faut savoir que en langage C, les tableaux sont stoqués ligne par ligne
Cela veut dire que les accès au tableau peuvent être plus ou moins rapides dépendant de la manière dont ont les accède.
P.ex.:

Code :

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

1
2
3
for(i=0;i<N;i++)
for(j=0;j<N;j++)
a[j][i]=a[j-1][i]+i;

Sera plus lent que

Code :

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

1
2
3
for(i=0;i<N;i++)
for(j=0;j<N;j++)
a[i][j]=a[i-1][j]+i;

Car la mémoire est lue par bloc (une partie de la ligne).
Dans le premier code, chaque tour de boucle demande une lecture, par contre, dans le deuxième bout de code, la donnée nécessaire se trouvera déjà quelques fois dans le cache.

Si je n'ai pas été clair, n'hésites pas, mais je pense que cela fera déjà speeder un peu plus.

[Vincent Rogier]

Un petit truc en regardant vite fait. Il faut savoir que en langage C, les tableaux sont stoqués ligne par ligne
Cela veut dire que les accès au tableau peuvent être plus ou moins rapides dépendant de la manière dont ont les accède.
P.ex.:

Code :

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

1
2
3
for(i=0;i<N;i++)
for(j=0;j<N;j++)
a[j][i]=a[j-1][i]+i;

Sera plus lent que

Code :

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

1
2
3
for(i=0;i<N;i++)
for(j=0;j<N;j++)
a[i][j]=a[i-1][j]+i;

Car la mémoire est lue par bloc (une partie de la ligne).
Dans le premier code, chaque tour de boucle demande une lecture, par contre, dans le deuxième bout de code, la donnée nécessaire se trouvera déjà quelques fois dans le cache.

Si je n'ai pas été clair, n'hésites pas, mais je pense que cela fera déjà speeder un peu plus.

Sauf que le compilo aura généré un code différent ou lui même fait ses propres optimisations.

Le code généré n'est jamais une traduction du code écrit !

Ne pas se prendre la tête sur ce genre de détails, laisser le compilo faire ses optimisations tout seul comme un grand.

Mais effectivement, voir plutôt du coté du forum algo...

[souviron34]

une première optimisation de code qui t'accélérera pas mal (sans parler de l'algo) :

sort tout ce qui peut être sortie des boucles ..

Exemple :

Code :

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

1
2
3
4
5
6
7
8
9
    b1 = b[k] /* 1 calcul en moins à chaque boucle i */
    c1 = a[k][k] /* 1 calcul en moins à chaque boucle i */
    for(i=k+1;i<n;i++)  //triangularisation de la matrice
        {
            p=a[i][k]/c1;
            for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
            b[i]=b[i]-p*b1 ;
       }

[un passant]

merci jvai essayer de le faire
mais sans parler de l'algo que je pensais optimiser mon prog.
par exemple je declare la matrice K comme un float qui est code sur 32 bits
si ma geometrie a par exemple 41205 noeuds,
la taille de K sera de 41205*41205 soit 1 697 852 025
que je multiplie par 32 bits =54 331 264 800 bits seront alloues pour la matrice K
ca fait beaucoup.
et mes temps de calculs sont tres eleves environ une dizaine de minutes

je voudrais savoir si il n'y aurais pas des types d'elements differents ou autres choses pour gagner du temps.

[souviron34]

non

tu fais des divisions..

Le MINIMUM est float.

Mais en plus tu as des trucs faux :

Code :

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

if (a[][] == 0 )

ne marchera pas.

c'est

Code :

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

if ( fabs(a[][]) < FLT_EPSILON )

(et encore, fabs est pour les double, pour les float je me souviens plus le nom).

Le type naturel en C des réels est double.

Je crois nettement plus que c'est ton algo, ou alors que tu as surestimé ce que tu serais capable de faire sur ton ordi. (c'est pas pour rien qu'il y a des ordis spécialisés calculs).

[nicolas.sitbon]

(et encore, fabs est pour les double, pour les float je me souviens plus le nom).

Code :

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

1
2
3
4
#include <math.h>
double fabs(double x);
float fabsf(float x);
long double fabsl(long double x);

sauf si l'on invoque la macro générique "fabs" de tgmath.h dispo depuis C99.
sur ma machine par exemple :

Code :

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

#define fabs(Val) __TGMATH_UNARY_REAL_IMAG_RET_REAL (Val, fabs, cabs)

[un passant]

Code :

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

if ( fabsf(a[][]) < FLT_EPSILON )

Si je met ca, quel est la valeur de FLT_ESPILON

et la partie
if ( a[k][k]==0);
marche, je l'ai teste
peut etre c'est la chance du debutant

[souviron34]

normalement la valeur se trouve dans float.h (ou limits.h pour les anciens compilos).

En fait, c'est tout simplement dû à la précision. Un flottant ayant n décimales, avoir 0.0000000000000000000000000000000000 est pratiquement impossible. Donc, à une certaine précision près, tout nombre très petit est EGAL à zéro à toutes fin pratiques (y compris dans une division).

C'est la différence entre les maths et la physique..


En général, pour les double c'est de l'ordre de 10^-13 à 10^-17...

pour les float c'est environ la moitié (10^-7 à 10^-9)

[un passant]

ok merci
et est ce que mon calcul sur la taille allouee pour la matrice K est correct.
sachant que 1 octet=8bits je dois diviser mon resultat par 8 pour avoir la taille occupee en octets.