Voilà ma source:
/*******************************************************/
/* Ce programme met en oeuvre la méthode des moindres carrés */
/*******************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
float Npoints;
#define Npoints 351054
/*******************************************************************/
/* La fonction suivante appelle la fonction qui va inverser la matrice passée en */
/* argument et affiche la solution lorsque celle ci existe. */
/*******************************************************************/
void inv_mat( float mat[3][3], float Sol[3])
{
int i, err;
int gauss_pivmax(float mat[3][3], float Sol[3]);
err=gauss_pivmax(mat,Sol);
if(err==1)
{
printf("Solution:\n");
for(i=0; i<=2; i++)
printf("a%d=%f \n",i,Sol[i]);
}
else printf("Erreur, matrice singuliere\n");
}
/*******************************************************************/
/* La fonction suivante renvoie la valeur 1 si la matrice est inversible, 0 sinon. */
/* Quand la valeur 1 est renvoyée, la solution du système se trouve */
/* dans le tableau Sol. */
/*******************************************************************/
int gauss_pivmax(float mat[3][3], float Sol[3])
{
int i, j, k, l, err;
float pivot, coeff, s;
float t[3][4];
double max;
for(i=0; i<=2; i++)
{
for(j=0; j<=2; j++)
{
t[i][j]=mat[i][j];
t[i][3]=Sol[i];
}
}
err=1;
k=1;
while(err==1 && k<2)
{
max=fabs(t[k][k]);
l=k;
for(i=k+1; i<=2; i++)
{
if(max<fabs(t[i][k]))
{
max=fabs(t[i][k]);
l=i;
}
}
if(max!=0)
{
if(l!=k)
{
for(j=k; j<=3; j++)
{
pivot=t[k][j];
t[k][j]=t[l][j];
t[l][j]=pivot;
}
}
pivot=t[k][k];
for(i=k+1; i<=2; i++)
{
coeff=t[i][k]/pivot;
for(j=k+1; j<=3; j++)
t[i][j]-=coeff*t[k][j];
}
}
else err=0;
k++;
}
if(t[2][2]==0) err=0;
if (err==1)
{
Sol[2]=t[2][3]/t[2][2];
for(i=1; i>=0; i--)
{
s=t[i][3];
for(j=i+1; j<=2; j++)
s-=t[i][j]*Sol[j];
Sol[i]=s/t[i][i];
}
}
return(err);
}
/************************Programme principale************************/
/* La première colonne contient les N valeurs de X, la deuxième colonne les N */
/* valeurs de Y et la troisième colonne les N valeurs de Z. */
/*****************************************************************/
void main (void)
{
int i;
float taille;
/**************************************** ***************/
/* La valeur de Npoints doit être réajustée pour chaque fichier */
/************ *******************************************/
float solution[3], tableau[3][3];
float xcarre[Npoints], xfoisy[Npoints], ycarre[Npoints], xfoisz[Npoints], yfoisz[Npoints];
float sumxcarre[Npoints], sumycarre[Npoints], sumxfoisy[Npoints], sumx[Npoints], sumy[Npoints], sumz[Npoints], sumxfoisz[Npoints], sumyfoisz[Npoints];
float coord[Npoints][3];
void inv_mat( float mat[3][3], float Sol[3]);
FILE * fich;
char nomfich[21];
float a,b,c;
/* long taille_fich, lSize;*/
printf("*************Methode des moindres carres*****************\n");
/*************************************/
/* Ouverture du fichier de données: */
/*************************************/
printf("nom du fichier a ouvrir:\n");
scanf("%20s",&nomfich);
/*printf("%s\n",nomfich);*/
fich=fopen(nomfich,"r");
if (fich == 0)
printf("impossible d'ouvrir le fichier\n");
else printf("\nfichier ouvert\n\n");
/*******************************************************************/
/* Les instructions suivantes permettent de lister le contenu */
/* du fichier ouvert precedemment. */
/*******************************************************************/
i=0;
do
{
fscanf(fich,"%f %f %f",&a,&b,&c);
/*printf("%f %f %f\n",a,b,c);*/
coord[i][0]=a;
coord[i][1]=b;
coord[i][2]=c;
i=i+1;
}
while(!feof(fich));
/* printf("\n\n");
fseek(fich,0,SEEK_END);
lSize=ftell(fich);
printf("float :%d\n",sizeof(float));
printf("int :%d\n",sizeof(int));
printf("long :%d\n",sizeof(long));
printf("double :%d\n",sizeof(double));
taille_fich=lSize/sizeof(float);
rewind (fich);
printf("taille du fichier:%hd elements.\n",taille_fich);*/
fclose(fich);
/* printf("\n\n Matrice des coordonnees:\n\n");
for(i=0; i<=Npoints-1;i++)
printf("%f %f %f\n",coord[i][0],coord[i][1],coord[i][2]);*/
/********************************************************************************/
/* Ici on calcule les différentes valeurs de X², Y², XY, XZ, YZ, nécessaire à */
/* la mise en oeuvre de la méthode des moindres carrés. */
/********************************************************************************/
for (i=0; i<=Npoints-1; i++)
{
xcarre[i]=coord[i][0]*coord[i][0];
ycarre[i]=coord[i][1]*coord[i][1];
xfoisy[i]=coord[i][0]*coord[i][1];
xfoisz[i]=coord[i][0]*coord[i][2];
yfoisz[i]=coord[i][1]*coord[i][2];
}
/*************************************************************************************/
/* Maintenant on calcule les différents termes qui interviennent dans la forme */
/* matricielle du problème. On commence par initialiser les valeurs avant de faire */
/* les sommes. */
/*************************************************************************************/
sumxcarre[0]=xcarre[0];
sumycarre[0]=ycarre[0];
sumxfoisy[0]=xfoisy[0];
sumxfoisz[0]=xfoisz[0];
sumyfoisz[0]=yfoisz[0];
sumx[0]=coord[0][0];
sumy[0]=coord[0][1];
sumz[0]=coord[0][2];
for (i=0; i<=Npoints-2; i++)
{
sumxcarre[i+1]=sumxcarre[i]+xcarre[i+1];
sumycarre[i+1]=sumycarre[i]+ycarre[i+1];
sumxfoisy[i+1]=sumxfoisy[i]+xfoisy[i+1];
sumxfoisz[i+1]=sumxfoisz[i]+xfoisz[i+1];
sumyfoisz[i+1]=sumyfoisz[i]+yfoisz[i+1];
sumx[i+1]=sumx[i]+coord[i+1][0];
sumy[i+1]=sumy[i]+coord[i+1][1];
sumz[i+1]=sumz[i]+coord[i+1][2];
}
/* printf("Differents parametres calcules\n");
/*******************************************************************************************/
/* Les sommes contenant tous les termes nécessaires sont donc dans les grandeurs indicées */
/* [Npoints-1]. */
/*******************************************************************************************/
/**************************************************************************************/
/* Maintenant on va construire la matrice du problème ainsi que le vecteur solution: */
/**************************************************************************************/
/****************** ____Construction du vecteur solution: ****************************/
solution[0]=sumz[Npoints-1];
solution[1]=sumxfoisz[Npoints-1];
solution[2]=sumyfoisz[Npoints-1];
/* printf("%f %f %f\n",solution[0],solution[1],solution[2]);
/********************____Construction de la matrice (symétrique!) du problème:*********/
tableau[0][1]=tableau[1][0]=sumx[Npoints-1];
tableau[0][2]=tableau[2][0]=sumy[Npoints-1];
tableau[1][2]=tableau[2][1]=sumxfoisy[Npoints-1];
tableau[0][0]=Npoints;
tableau[1][1]=sumxcarre[Npoints-1];
tableau[2][2]=sumycarre[Npoints-1];
/* printf("Matrice et second membre du probleme calcules\n");
for(i=0;i<=2;i++)
printf("%f %f %f\n",tableau[i][0],tableau[i][1],tableau[i][2]);
/******************************************************************************************/
/* On passe la matrice et le vecteurs solutions à la subroutine chargée d'inverser A et */
/* de donner la solution. */
/******************************************************************************************/
inv_mat(tableau,solution);
printf("Fin de la methode des moindres carres\n");
}
}
(on ne fait pas [3][3] mais seulement [3]) bref, c comme ça 
)
