#include "malloc.h"
#include "stdio.h"
#include "stdlib.h"

void identite(double *mat,int dim);
void copie_matriceR(double *mat1,double *mat2,int dim);
void gaussR(double *mat1,double *mat2,int dim);

void mat_random(double *mat,int dim) {
 int i,j;
 for(i=0;i<dim;i++) {
  for(j=0;j<dim;j++) {
   mat[i+j*dim]= ((double)rand() / ((double)(RAND_MAX)+(double)(1)) )*10.0;
  }
 }
}

void print_mat(double *mat, int dim) {
 int i,j;
 printf("Matrice %d*%d\n",dim,dim);
 for(i=0;i<dim;i++) {
  for(j=0;j<dim;j++) {
   printf("%g ",mat[i+j*dim]);
  }
  printf("\n");
 }
}

/* fonction qui rempli une matrice */
/* par la matrice identité */
void identite(double *mat,int dim) {
 int i,j;
 for(i=0;i<dim;i++) {
  for(j=0;j<dim;j++) {
   if(i==j) {
    mat[i+j*dim]=1.;
   } else {
    mat[i+j*dim]=0.;
   }
  }
 }
}

/* fonction qui copie une matrice dans une autre */
void copie_matriceR(double *mat1,double *mat2,int dim) {
 int i,j;

 for(i=0;i<dim;i++) {
  for(j=0;j<dim;j++) {
   mat2[i+j*dim]=mat1[i+j*dim];
  }
 }
}

/*la matrice mat1 est celle qui est à inverser
//et l'inverse est contenu dans mat2 */
void gaussR(double *mat1,double *mat2,int dim) {
 double *temp;
 double a,b;
 int c,i,j,k,q,l;
 temp=(double *)malloc(sizeof(double)*dim*dim);
 copie_matriceR(mat1,temp,dim);
 a=0.;
 b=0.;
 c=0;
 identite(mat2,dim);
 for(k=0;k<dim;k++) {
  a=temp[k+k*dim];
/*verifie la condition "a" different de zero
//sinon on copie une ligne pour resoudre le probleme */
  c=0;
  while(abs(a)<0.000000001) {
   c++;
   for(q=0;q<dim;q++) {
    temp[k+q*dim]=temp[k+q*dim]+temp[k+c+q*dim];
    mat2[k+q*dim]=mat2[k+q*dim]+mat2[k+c+q*dim];
   }
   a=temp[k+k*dim];
  }
/*normalisation la ligne k */
  for(l=0;l<dim;l++) {
   temp[k+l*dim]=temp[k+l*dim]/a;
   mat2[k+l*dim]=mat2[k+l*dim]/a;
  }
/*reduction de gauss-jordan */
  for(i=0;i<dim;i++) {
   b=temp[i+k*dim];
   if(i!=k) {
    for(j=0;j<dim;j++) {
     temp[i+j*dim]=temp[i+j*dim]-b*temp[k+j*dim];
     mat2[i+j*dim]=mat2[i+j*dim]-b*mat2[k+j*dim];
    }
   }
  }
 }
 free(temp);
}

int main() {
/*exemple d'utilisation pour une matrice 2x2 */
 int dim=100;
 double *M;
 double *M_inv;
 int i,j;
 M=(double *)malloc(sizeof(double)*dim*dim);
 M_inv=(double *)malloc(sizeof(double)*dim*dim);

/*remplissage de la matrice M */
/* M[0+0*dim]=2.;
 M[0+1*dim]=4.8;
 M[1+0*dim]=6.2;
 M[1+1*dim]=5.1;
*/
 mat_random(M,dim);
 print_mat(M,dim);
/*inversion de la matrice M */
 gaussR(M,M_inv,dim);
/*affichage des elements de matrice inverse*/
 print_mat(M_inv,dim);
}