/*
  gauss.c  -  Viktor Leijon(leijon@ludd.luth.se) 1997-99

  Feel free to use this program in any way you wish, with the condition
  that you first notify me of it. Thank you.

  Calculation of inverse matrixes by Gauss-Jordan elimination  

  (Matris means Matrix in Swedish, thus some of the variable names)

*/

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

#define MAXMATRIS 20

typedef float matristyp[MAXMATRIS][MAXMATRIS]; /* x, y*/

/* This functions prints the matrix matrix in screen */
void printmatrix(matristyp matrix,int size) {
 int i,j;

 for(i=0;i<size;i++) {
   for (j=0;j<size;j++) {
     printf("%f ",matrix[j][i]);
   }
   printf("\n");
 }
}

/* This function exchanges row1 and row2 in the matrix matrix*/
void chrow(matristyp matrix,int row1,int row2,int size) {
  int counter;
  float temp;

  for (counter=0;counter<size;counter++) {
    temp = matrix[counter][row1];
    matrix[counter][row1] = matrix[counter][row2];
    matrix[counter][row2] = temp;
  }
  
}

/* Reduces row "rowtoreduce" with factor*"reducer" */
void reduce(matristyp matrix,int rowtoreduce,
	    int reducer,float factor,int size)
{
  int i;

  for (i=0;i<size;i++){
    matrix[i][rowtoreduce] -= factor * matrix[i][reducer];
  }
  
}

/* Divdes row "row" by factor.. (equal to multiplying with 1/factor,
   a valid row operation) */
void divide(matristyp matrix,int row,float factor,int size){
  int counter;

  for(counter=0;counter<size;counter++) {
    matrix[counter][row] /= factor;
  }

}

/* Calculates the inverse of matrix matrix, prints it.. */
int inverse(matristyp matrix,int size)
{
  matristyp identitet; /* "identity" matrix */
  int x,y;
  int gausspos,counter;
  float koeff,fact;
  
  printf("\nOriginal matrix : \n");
  printmatrix(matrix,size);

  /* Skapa identitets matris .. */
  for (x=0;x<size;x++)
    for (y=0;y<size;y++) {
      if (x == y)
	identitet[x][y] = 1;
      else
	identitet[x][y] = 0;
    }

  /* Gauss elimination */
  for (gausspos=0;gausspos<size;gausspos++) {
    koeff = matrix[gausspos][gausspos];

    if (koeff == 0) {
      for(counter=(gausspos+1);counter <size;counter++) {
	if (matrix[gausspos][counter] != 0)
	  {
	    chrow(matrix,gausspos,counter,size);
	    chrow(identitet,gausspos,counter,size);
	    koeff = matrix[gausspos][gausspos];
	    break;
          }
      }
      if (koeff == 0) {
	printf("Cannot invert matrix.\n");
	printmatrix(matrix,size);
	return(0);
      }
    } /* if koeff== 0 */

    divide(matrix,gausspos,koeff,size);
    divide(identitet,gausspos,koeff,size);
    
    for (counter=(gausspos+1);counter <size;counter++) {
      fact = matrix[gausspos][counter];
      reduce(matrix,counter,gausspos,fact,size);
      reduce(identitet,counter,gausspos,fact,size);
    }
  } /* Eliminiation */
    
  /* Resubstitution.. */
  for (gausspos=(size-1);gausspos>=0;gausspos--) {
    for (counter=(gausspos-1);counter >=0;counter--) {
      fact = matrix[gausspos][counter];
      reduce(matrix,counter,gausspos,fact,size);
      reduce(identitet,counter,gausspos,fact,size);
    }
  }
    
  printf("Inverted matrix : \n");
  printmatrix(identitet,size);
}


int main()
{
  matristyp matrix;
  int m,i,j;
  
  printf("Creating matrix (mxm)\nEnter m : ");
  scanf("%d",&m);
  printf("Enter the matrix..\n");
  for (i=0;i<m;i++) {
    for (j=0;j<m;j++) {
      printf("Enter value for row %d, col %d : ",i+1,j+1);
      scanf("%f",&matrix[j][i]);
    }
  }
  inverse(matrix,m);
  return 0;
}