ガウスの消去法（実装）

ガウスの消去法の実装を行う。

コードは次のようになる。

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

/************************************************************
 * ガウスの消去法による線形方程式の数値計算
 * 
 * (C) 2012- Gofer.
 ************************************************************/

// 方程式の表示
void display_equtation(int N, double *A) {
 int i, j;
 
 for(i=0; i<N; ++i) {
  for(j=0; j<N; ++j) {
   printf("%f x[%d]", A[i*(N+1) + j], j);
   if(j < N-1) printf(" + ");
  }
  printf(" = %f\n", A[i*(N+1) + N]);
 }
 putchar('\n');
 
 return;
}

// n行目とm行目を入れ替える
double* replace_row(int n, int m, int N, double* A) {
 double* line_temp = (double*)malloc(sizeof(double)*(N+1));
 memmove(line_temp,   &A[n*(N+1)], sizeof(double)*(N+1));
 memmove(&A[n*(N+1)], &A[m*(N+1)], sizeof(double)*(N+1));
 memmove(&A[m*(N+1)], line_temp,  sizeof(double)*(N+1));
 free(line_temp);
 return A;
}
 
// i行目の前進消去を行う
double* forward_elimination(int i, int N, double *A) {
 int j, k, J = -1;
 double divide;
 
 // 対角要素が0ならピボット選択を行う
 if(A[i*(N+1) + i] == 0) {
  for(j=i+1; j<N; ++j) {
   if(A[j*(N+1) + i] != 0) {
    J = j;
    break;
   }
  }
  if(J == -1) return NULL;
  replace_row(i, J, N, A);
 }
 
 // 割算する係数を決定
 divide = A[i*(N+1) + i];
 
 // 標準化
 for(j=0; j<(N+1); ++j)
  A[i*(N+1) + j] /= divide;
 
 // i行目以外の行から引く
 for(k=i+1; k<N; ++k) {
  for(j=(N+1)-1; j>=0; --j)
   A[k*(N+1) + j] -= A[k*(N+1) + i] * A[i*(N+1) + j];
 }
 
 return A;
}

// i行目の後退代入
double* backward_substitution(int i, int N, double *A, double *x) {
 // (N-1)-i個の解は決定済み
 int j, k;
 
 // 1ステップ前に求めた解を代入
 k = i+1;
 if(k < N) { 
  for(j=k-1; j>=0; --j) {
   A[j*(N+1) + N] -= x[k] * A[j*(N+1) + k];
   A[j*(N+1) + k] = 0;
  }
 }
 
 // 第i行を用いて1文字消去
 if(A[i*(N+1) + i] != 0) {
  A[i*(N+1) + N] /= A[i*(N+1) + i];
  A[i*(N+1) + i] = 1;
 }
 
 // xベクタに代入
 x[i] = A[i*(N+1) + N];
 
 return A;
}

int main(void) {
 //int N = 2;
 /*
 double A[2][3] = { {1, 2, 5}, 
      {2, 3, 8}};
 */

 //int N = 3;
 /*
 double A[3][4] = { {2, 4, 2, 8}, 
      {4, 10, 3, 17}, 
      {3, 7, 1, 11}};
 */

 //int N = 4;
 /*
 double A[4][5] = { {1, 1, 2, 1, 1}, 
      {1, 1, 3, 2, -2}, 
      {2, -2, 2, -1, 1}, 
      {-1, 1, 0, 1, -1}};
 */
 
 int N = 4;      
 double A[4][5] = { {1, -1, 2, -1, -8}, 
      {2, -2, 3, -3, -20}, 
      {1, 1, 1, 0, -2}, 
      {1, -1, 4, 3, 4}};
 
 int i;
 double* x = (double*)malloc(sizeof(double) * N);
 
 display_equtation(N, (double*)A);
 
 // 前進消去
 for(i=0; i<N; ++i) {
  if(forward_elimination(i, N, (double*)A) == NULL) {
   fprintf(stderr, "Error!\n");
   free(x);
   return -1;
  }
 }
 
 // 後退代入
 for(i=N-1; i>=0; --i)
  backward_substitution(i, N, (double*)A, x);
 
 for(i=0; i<N; ++i)
  printf("x[%d] = %f\n", i, x[i]);
 
 free(x);
 
 return 0;
}