連立方程式を部分pivot選択付きガウスの消去法をC言語で書いた

C言語によるガウスの消去法

もう数値計算の世界では一億番煎じになってると思われるガウスの消去法を自分で書いた
ので

 //2x 5y 7z =4
//4x 13y 20z =11
//8x 29y 50z =29

# define N 3
# define M N+1
# include<stdio.h>
# include<iostream>
using namespace std;

int main(){
  double mat[N][N+1]={{2,5,7,4},{4,13,20,11},{8,29,50,29}};
  double mat_out[N];
  double pivot;
  int pivot_ren=0;
  int i,j,t,k;


  for(i=0; i < N;i++){
      pivot= mat[i][i];
      for(k=i+1;k < N;k++){
          if(fabs(pivot) < fabs( mat[k][i])) {
              pivot=mat[k][i];
              pivot_ren=k;
          }
      }
      if(pivot==0){
          printf("エラー！");
          return -1;
      }
      if(i!=pivot_ren){
          for(j=0;j<M;j++){
              pivot=mat[i][j];
              mat[i][j]=mat[pivot_ren][j];
              mat[pivot_ren][j]=pivot;
          }
      }
    
      for(t=i+1;t<N;t++){
          pivot = mat[t][i] / mat[i][i];
          for(j=0;j< M; j++){
              mat[t][j] = mat[t][j] - pivot * mat[i][j];
          }
      }
      mat_out[i] = mat[i][N];
    }

    for( i = N; i >= 0; i--) {
       for(t = N-1;t > i;t--){
        mat_out[i] = mat_out[i] -  (mat[i][t] * mat_out[t]);
        }
      mat_out[i] = mat_out[i] / mat[i][i];
    }
    
   for(i=0;i<N;i++){
    for(j=0;j<N+1;j++){
        printf("%lf ",mat[i][j]);
    }
    printf("\n");
}

for( i = 0; i < N; i++)
{
    cout << mat_out[i] << endl;
}

return 0;    
}

出力結果

1
-1
1

$x=1,y=-1,z=1$が解なので正解

cpp使ってるけど最後の解の出力以外はC言語なので実質C言語ですね。
コンパイルはgcc 5.4.0でした
ソースコード内部にはゼロ除算で$∞$になる成分が出ると計算が成り立たないのでそのようなことがないようにエラー処理を行っています。

アルゴリズムの概要

はっちゃけていうとガウスの消去法は掃き出し法と代入法のハーフ。
正確に言うと下三角部分を掃き出し法で掃出して未知変数をなるべく減らした後に
最下行から上三角部分を代入法で解を求めていく方法です

下三角部分の掃き出し

これを数値計算では前進消去と呼称するそうです

 for(t=i+1;t<N;t++){
          pivot = mat[t][i] / mat[i][i];//ここで対角成分を抽出します
          for(j=0;j< M; j++){
              mat[t][j] = mat[t][j] - pivot * mat[i][j];
          }
      }
      mat_out[i] = mat[i][N];
    }

ソースコード中の上記の部分がこれに当たります
中学生の連立方程式を解く方法とやってることは基本同じです。
加減法を下三角部分にだけ行うと考えればいいです。
掃き出し済みの行の拡大係数成分をmat_out[i]に格納します

上三角部分の後退代入（代入法）

さっき前進消去で格納したmat_out[i]を利用して計算します
最下行の解はすでにこの時点でも求まっているので一つ上の行の解を求める為にいわゆる
「代入法」を使います

 for( i = N; i >= 0; i--) {
       for(t = N-1;t > i;t--){
        mat_out[i] = mat_out[i] -  (mat[i][t] * mat_out[t]);
        }
      mat_out[i] = mat_out[i] / mat[i][i];
    }
    

ソース中ではこの部分です

部分pivot選択付きとは

計算誤差を防ぐために前進消去を絶対値の大きい対角成分を持つ行から行っていく処理です。
なんで計算誤差が絶対値が小さいと起こりやすくなるのかはIEE754 floatの仕様から説明する必要が出てきてしまうのであまりここでは深く触れません。

  for(i=0; i < N;i++){
      pivot= mat[i][i];
      for(k=i+1;k < N;k++){
          if(fabs(pivot) < fabs( mat[k][i])) {
              pivot=mat[k][i];
              pivot_ren=k;
          }
      }
      if(pivot==0){
          printf("エラー！");
          return -1;
      }
      if(i!=pivot_ren){
          for(j=0;j<M;j++){
              pivot=mat[i][j];
              mat[i][j]=mat[pivot_ren][j];
              mat[pivot_ren][j]=pivot;
          }
      }

ソース中ではこの部分です
絶対値の大きいpivot成分を決める部分と本来のpivot行と絶対値が大きいpivot成分を持つ
成分を交換する部分で構成されています

この記事について

非常に適当に書かれてています、自分の理解を確認するためのモノでもあるのでもっと詳しい文献が大量に巷には存在しますのでそちらを参考にしてください。

参考文献

[http://nkl.cc.u-tokyo.ac.jp/12s/3D/pivoting.pdf]
[http://nkl.cc.u-tokyo.ac.jp/13n/SolverDirect.pdf]
どちらも素晴らしいテキストなので是非読んでみてください