OPÉRATION SUR MATRICE

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


typedef struct{
	int n, m; /*n->lignes, m->colonnes*/
  	double **P;/*Matrice*/
} *Matrice, structMatrice;


Matrice allouerMatrice(int l,int c){/*allouer une matrice de l lignes et c colonnes*/
	Matrice A;
	int i;
	A=(Matrice)malloc(sizeof(structMatrice));
	A->n=l;A->m=c;
	A->P=(double**)malloc(l*sizeof(double*));
	for(i=0;i<l;i++) {
		A->P[i]=(double *)malloc(c*sizeof(double ));
	}
	return A;
}

void libererMatrice(Matrice A){/*libÃ&#168;re l'espace mémoire occupé par la matrice*/
	int i;
	for(i=0;i<A->n;i++)free(A->P[i]);
}
int deMemeDim(Matrice A, Matrice B){/*teste si les matrices A et B sont de même dimension*/
	if(A->n == B->n && A->m == B->m)return 1;
	return 0;
}

Matrice Mat_add(Matrice A, Matrice B){/*renvoie la matrice issue de l'addition de A et de B*/
	int i,j;
	Matrice C;
	if(!deMemeDim(A,B))return NULL;
	C=allouerMatrice(A->n,A->m);
	for(i=0;i<A->n;i++){
		for(j=0;j<A->m;j++){
			C->P[i][j]=A->P[i][j]+B->P[i][j];
		}
	}
	return C;
}

Matrice Mat_Id(int n){/*renvoie la matrice identité de dimension n x n */
	int i,j;
	Matrice C;
	C=allouerMatrice(n,n);
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			if(i!=j)C->P[i][j]=0;
			else C->P[i][j]=1;
		}
	}
	return C;
}

Matrice Mat_nulle(int n){/*renvoie une matrice nulle*/
	int i,j;
	Matrice C;
	C=allouerMatrice(n,n);
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			C->P[i][j]=0;
		}
	}
	return C;
}
Matrice Mat_mul(Matrice A, double c){/*renvoie le produit de A par le scalaire c*/
	int i,j;
	Matrice C;
	C=allouerMatrice(A->n,A->m);
	for(i=0;i<A->n;i++){
		for(j=0;j<A->m;j++){
			C->P[i][j]=c*A->P[i][j];
		}
	}
	return C;
}



Matrice Mat_dot(Matrice A, Matrice B){/*renvoie le produit de A et B,
	renvoie NULL si l'opération n'est pas possible*/
	int i,j,k;
	Matrice C;
	if(A->m!=B->n)return NULL;
	C=allouerMatrice(A->n,B->m);
	for(i=0;i<C->n;i++){
		for(j=0;j<C->m;j++){
			for(k=0;k<A->m;k++){
			C->P[i][j]+= A->P[i][k]*B->P[k][j];
			}
		}
	}
	
	return C;
}

Matrice Mat_transpose(Matrice A){/*renvoie la transposée de la matrice A*/
	int i,j;
	Matrice C;
	C=allouerMatrice(A->m,A->n);
	for(i=0;i<C->n;i++){
		for(j=0;j<C->m;j++){
			C->P[i][j]=A->P[j][i];
		}
	}
	return C;
}



Matrice Mat_sousMat(Matrice A,int ib,int jb){/*renvoie la matrice A à laquelle à été supprimé
	la ib°ligne et la jb° colonne*/
	int i,j;
	int is=0,js=0;
	Matrice B;
	B=allouerMatrice(A->n-1,A->m-1);
	for(i=0;i<A->n;i++){
		for(j=0;j<A->m;j++){
			if(i!=ib && j!=jb){
			B->P[is][js]=A->P[i][j];
			if(js<B->n-1)js++;
			else{js=0;is++;}	
			}
		}
	}
	return B;
}
double Mat_determinant(Matrice A){/*renvoie le déterminant de la matrice A*/
	int i,j;
	double det=0;
	Matrice B;
	//condition d'arret
	if(A->n==1)return A->P[0][0];
	if(A->n==2)return A->P[0][0]*A->P[1][1]-A->P[0][1]*A->P[1][0];
	//condition de descente dans la récurrence
	else{
		for(j=0;j<A->m;j++){
		B=Mat_sousMat(A,0,j);
		if(j%2==0){det+=A->P[0][j]*Mat_determinant(B);}
		else{det+=(-1)*A->P[0][j]*Mat_determinant(B);}
		libererMatrice(B);
		}
	}
	return det;
}
void Mat_afficherMat(Matrice A){/*affiche la matrice A*/
	int i,j;
	printf("\n Matrice %dx%d :\n",A->n,A->m);
	for(i=0;i<A->n;i++){
		for(j=0;j<A->m;j++){
			printf("%f\t",A->P[i][j]);
		}
		printf("\n");
	}
}
void Mat_remplirMatrice(Matrice A){/*remplie la matrice A par saisie au clavier*/
	int i,j;
	printf("\nRemplissage de matrice %dx%d\n",A->n,A->m);
	for(i=0;i<A->n;i++){
		for(j=0;j<A->m;j++){
			printf("element (%d,%d)?",i,j);
			scanf("%lf",&(A->P[i][j]));
		}
	}
}

Matrice Mat_coMat(Matrice A){/*renvoie la comatrice de A*/
	int i,j;
	Matrice B;
	Matrice S;
	B=allouerMatrice(A->n,A->m);
	//S=allouerMatrice(A->n-1,A->m-1);
	for(i=0;i<B->n;i++){
		for(j=0;j<B->m;j++){
			S=Mat_sousMat(A,i,j);
			if((i+j)%2==0){B->P[i][j]=Mat_determinant(S);}
			else{B->P[i][j]=-1*Mat_determinant(S);}
			libererMatrice(S);
		}
	}
	return B;
}
Matrice Mat_invMat(Matrice A){/*renvoie la matrice inverse de A*/
	double det;
	Matrice coA,t_coA,inv_A;
	det=Mat_determinant(A);
	if(det==0)return NULL;
	coA=Mat_coMat(A);
	t_coA=Mat_transpose(coA);
	inv_A=Mat_mul(t_coA,1/det);
	libererMatrice(coA);
	libererMatrice(t_coA);
	libererMatrice(inv_A);
	return inv_A;
}

int Mat_moindreCarre(Matrice A,Matrice *U,Matrice B){/*résout A.U=B par la méthode de moindre carré
	renvoie 1 si le calcul est possible, 0 sinon*/
	Matrice At,AtA,AtA_inv,AtA_inv_At;
	At=Mat_transpose(A);
	AtA=Mat_dot(At,A);
	AtA_inv=Mat_invMat(AtA);
	if(AtA_inv==NULL)return 0;
	AtA_inv_At=Mat_dot(AtA_inv,At);
	*U=Mat_dot(AtA_inv_At,B);
	libererMatrice(At);
	libererMatrice(AtA);
	libererMatrice(AtA_inv);
	libererMatrice(AtA_inv_At);
	return 1;
}
int main(){
	Matrice A,B,C,D,U;
	double det;
	int i,j;
	/*A=allouerMatrice(3,3);
	Mat_remplirMatrice(A);
	Mat_afficherMat(A);
	getchar();
	B=Mat_mul(A,3.2);
	Mat_afficherMat(B);
	getchar();
	libererMatrice(B);
	B=Mat_Id(3);
	C=Mat_dot(A,B);
	Mat_afficherMat(C);
	getchar();
	det=Mat_determinant(A);
	printf("det(A)=%lf\n",det);
	getchar();
	A=Mat_mul(A,2);
	Mat_afficherMat(A);
	getchar();	
	D=allouerMatrice(2,2);
	Mat_remplirMatrice(D);
	Mat_afficherMat(D);
	det=Mat_determinant(D);
	printf("det=%lf\n",det);
	D=Mat_invMat(D);
	D=Mat_mul(D,det);
	Mat_afficherMat(D);*/
	//Test Méthode de moindre carré
	A=allouerMatrice(10,4);
	B=allouerMatrice(10,1);
	for(i=0;i<A->n;i++){
		B->P[i][0]=(double)rand()*10/RAND_MAX;
		for(j=0;j<A->m;j++){
			A->P[i][j]=(double)rand()*10/RAND_MAX;
		}
	}
	
	Mat_afficherMat(A);
	Mat_afficherMat(B);
	Mat_moindreCarre(A,&U,B);
	Mat_afficherMat(U);
	libererMatrice(A);
	libererMatrice(B);
	libererMatrice(U);
}