// ガウスの消去法を行う関数を作成せよ。引数として拡大係数行列を受け取ることとする。

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

#define MN_MAX 100
#define EPS 1.0E-3 //許容誤差

struct matrix_data{
  int row; //行数
  int col; //列数
  double dat[MN_MAX];
  
};

struct matrix_data AugmentedMatrix = {
  3,
  4,
  {
    2.0,  3.0, -1.0,  1.0,
    1.0,  1.0,  2.0,  0.0,
    3.0, -1.0,  1.0,  2.0,
  }
};

struct matrix_data init_matrix(int m, int n)
{
  struct matrix_data mat;

  mat.row = m;
  mat.col = n;

  return mat;
}


void set_matrix(struct matrix_data* mat, int i, int j, double val)
{
  mat->dat[i*mat->col + j] = val;
}

double get_matrix(struct matrix_data mat, int i, int j)
{
  return mat.dat[i*mat.col + j];
}

//p行とq行を入れ替える
struct matrix_data exchange_matrix(struct matrix_data mat,int p, int q)
{
  int i;
  double dp,dq;

  for(i = 0; i < mat.col; ++i){
    dp = get_matrix(mat, p, i);
    dq = get_matrix(mat, q, i);
    set_matrix(&mat, q, i, dp);
    set_matrix(&mat, p, i, dq);
  }
  return mat;
}

//行列の表示
void show_matrix(struct matrix_data mat)
{
  int i, j;

  for (i = 0; i < mat.row; i++){
    for(j = 0; j < mat.col; j++){
      printf("%f ",get_matrix(mat, i, j));
    }
    printf("\n");
  }
}

//ガウスの消去法
struct matrix_data Gaussian_Elimination(struct matrix_data mat)
{
  int i,j,k;
  int n;
  double p,q;
  
  struct matrix_data x; 
  x = init_matrix(mat.col-1,1);
  n = mat.col;  
  
  // 前進消去
  for(k = 0; k < n-1; k++){
    for(i = k+1; i < n-1; i++){
      if(fabs(get_matrix(mat, k, k)) < fabs(get_matrix(mat, i, k))){
	mat = exchange_matrix(mat, i, k);
      }
    }
    if(fabs(get_matrix(mat, k, k)) < EPS){
      exit(2);
    }

    p = get_matrix(mat, k, k);
    for(j = k; j < n; j++){
      set_matrix(&mat, k, j, get_matrix(mat, k, j) / p);
    }
    
    for(i = k+1; i < n-1; i++){
      q = get_matrix(mat, i, k);
      for(j = k; j< n; j++){
	set_matrix(&mat, i, j, get_matrix(mat, i ,j) - get_matrix(mat, k, j)*q);
      }
    }
  }
  printf("途中経過\n");
  show_matrix(mat);

  // 後退代入
  for(k = n-2; k >= 0; k--){
    q = get_matrix(mat, k, n-1);
    for(j = k+1; j < n-1; j++){
      q = q - get_matrix(mat, k ,j) * get_matrix(x, j, 0);
    }
    set_matrix(&x, k, 0, q);
  }
  return x;
}


int main(void)
{
  struct matrix_data x;

  printf("もとの行列を表示します\n");
  show_matrix(AugmentedMatrix);

  x = Gaussian_Elimination(AugmentedMatrix);
  printf("解を表示します。\n");
  show_matrix(x);

  return 0;
}