ガウスの消去法

連立一次方程式を解くというような話になった時に一番最初に出てくるであろうガウスの消去法から取り上げていきたいと思います。おそらくこれは中高生のころに習った方程式の解き方みたいな感じなので理解はしやすいです(筆者がそうでした)。
何をしているのかを簡単に説明すると、連立なので複数の式があります。これらの式を足したり、引いたり、掛けたうえで加減したすることで変数を消していくというようなものです。最終的に一つの変数の値がわかれば、芋づる式に他の変数の値が求められるようになります。前半部分の変数を消していく作業は数式の左側の変数から順次やっていって、最終的に方程式の係数を表す行列である係数行列が上三角にするのが目的です。
具体的に例を使って説明してみます。

$3 x_0 + 5 x_1 - 7 x_2 = 0 \\ x_0 + x_1 - x_2 = 2 \\ 2 x_0 - 3 x_1 + x_2 = 5$
この方程式を解くことを考えます。
まず $x_0$ に着目します。
すると、2行目 - 1行目 * 1/3を行うことによって、2行目の $x_0$ は消え、3行目 - 1行目 * 2/3を行うことで3行目の $x_0$ は消えます。
このように1つずつ変数を消していく作業を行うのですが、次の $x_1$ については、3行目にだけ行ってあげ、上の1行目には行いません。これは後で簡単に値が求められるようになるので、わざわざ無駄な計算をする必要は無いというものです。ちなみに逐一上の行の変数も消していく手法はガウス・ジョルダン法といいます。
最終的に
$3 x_0 + 5 x_1 - 7 x_2 = 0 \\ -2/3 x_1 + 4/3 x_2 = 2 \\ -7 x_2 = 14$
となります。
ここまでの操作を前進消去(Forward elimination)といいます。
ここから1つずつxの値を求めていく作業になりまして、これを後進代入(Backward substitution)といいます。
下の方の数式から順番に値を求めていくだけという簡単な話で
$x_2 = -2, x_1 = 1, x_0 = 3$ と計算できます。

このような形で求めていくのがガウスの消去法で非常に簡単に理解できるものです。
ただ、係数行列の対角成分がゼロであるような場合にはうまく消去が出来ません。そういった際には対象とする数式を変更するというようなことを行います。ピボット選択と言われる処理になります。ここでは面倒なので割愛します。

ここからが本番でプログラムを書いてみます。
ちなみに行列の表記についてですが、私は2次元配列を使うのがあまり好きではないので、いつも1次元で表しています。この際にメモリの連続性をどう取るかという話が出てきます。行方向に連続に取るか、列方向に連続に取るかというものです。ここでは行方向に連続になるように配置を考えます。つまり、メモリ上には $a_{1,1}, a_{1,2}, a_{1,3},...,a_{1,n}, a_{2,1},...$ と並んでいる状態になります。
解法部分のプログラムが下のものになります。
ちなみにrealは浮動小数点数を表す型でこちらで作成したものです。また、行列はMxNのものを考えています。

void solve_gauss(real *dmat, real *x, real *y, int M, int N) {

  int i, j, k;
  real d;
  real *eq;
  
  eq = (real *)malloc(sizeof(real) * M * (N + 1));
  for (i = 0; i < M; i++) {
    for (j = 0; j < N; j++) {
      eq[i * (N + 1) + j] =dmat->val[i * N + j];
    }
    eq[i * (N + 1) + N] = y[i];
  }
 
  /* Forward elimination */
  for (i = 0; i < M; i++) {
    for (j = i + 1; j < M; j++) {
      d = eq[j * (N + 1) + i] / eq[i * (N + 1) + i];
      for (k = i; k <= M; k++) {
        eq[j * (N + 1) + k] -= d * eq[i * (N + 1) + k];
      }
    }
  }

  /* Backward substitution */
  for (i = M - 1; i >= 0; i--) {
    d = eq[i * (N + 1) + N];
    for (j = i + 1; j < M; j++) {
      d -= eq[i * (N + 1) + j] * eq[j * (N + 1) + M];
    }
    eq[i * (N + 1) + M] = d / eq[i * (N + 1) + i];
  }
  
  /* Copy to x vector */
  for (i = 0; i < M; i++) {
    x[i] = eq[i * (N + 1) + N];
  }
}

順を追って最初から説明していきますと、まずはじめに行列と右辺ベクトルから数式eqを作成します。その後、前進消去を行い、上三角行列を作ったら、後進代入を行っていきます。最終的に右辺ベクトルが合った場所に求めるxの解がありますので、それをxにコピーしてあげることで終了となります。

配列のアドレス計算に気をつけてあげれば、難しくはないプログラムです。

ガウスの消去法では前進消去のところに関しては、並列化が簡単に出来そうだなあと思いました。2番目のjループを並列化すればある程度性能は上がりそうですね。iが大きくなると並列度は下がりますが。

というわけでざっくりした感じでしたが、ガウスの消去法についてはここまでということで。

全然関係ないけど、数式うまく書けるようにしときたいなあ