##連立一次方程式のGuass-Jordan法
Guass-Jordan法は、連立一次方程式の数値解を求める代表的な方法のひとつで、前進のみの消去法で後進部がないという特徴があります。
この方法はプログラムが簡単で、一般に制度もいいので広く用いられています。
今回は以下の行列の$x_1,x_2,\dots x_n$の値を計算します。

\begin{pmatrix} 
  a_{11}x_1 &+& a_{12}x_2 &+&a_{13}x_3&+& \dots &+& a_{1n}x_n &= & b_1\\
   a_{21}x_1 &+& a_{22}x_2 &+&a_{23}x_3&+& \dots &+& a_{2n}x_n &= & b_2\\
   &  & \vdots &  &  &&&&&& \\
   a_{n1}x_1&+& a_{n2}^{(1)}x_2 &+&a_{n3}^{(1)}x_3&+& \dots &+& a_{nn}^{(1)}x_n &= & b_n^{(1)}\\
   \tag{1}
\end{pmatrix}

これを計算するために、まずは(1)の第1式の両辺を$a_{11}$で割り、得られた第1式を用いて第$i$式($i$=2,3,...,n)の$x_1$の係数$a_{i1}$を削除します。

\begin{pmatrix} 
  x_1 &+& a_{12}^{(1)}x_2 &+&a_{13}^{(1)}x_3&+& \dots &+& a_{1n}^{(1)}x_n &= & b_1^{(1)}\\
   && a_{22}^{(1)}x_2 &+&a_{23}^{(1)}x_3&+& \dots &+& a_{2n}^{(1)}x_n &= & b_2^{(1)}\\
   &  & \dots &  &  &&&&&& \\
   && a_{n2}^{(1)}x_2 &+&a_{n3}^{(1)}x_3&+& \dots &+& a_{nn}^{(1)}x_n &= & b_n^{(1)}\\
   \tag{2}
\end{pmatrix}

ここでの$a_{ij}^{(1)}$および$b_{ij}^{(1)}$は新しい係数を表します。
第二段階では、先ほどの式(2)の第2式を$a_{22}^{(1)}$で割り、得られた第2式を用いて第1式及び第$i$式$(i=3,4,...,n)$の$x_2$の係数を削除すると、

\begin{pmatrix} 
  x_1 &&&+&a_{13}^{(2)}x_3&+& \dots &+& a_{1n}^{(2)}x_n &= & b_1^{(2)}\\
   &&x_2 &+&a_{23}^{(2)}x_3&+& \dots &+& a_{2n}^{(2)}x_n &= & b_2^{(2)}\\
   &  &  & \dots &  &&&&&& \\
   &&&+&a_{n3}^{(2)}x_3&+& \dots &+& a_{nn}^{(2)}x_n &= & b_n^{(2)}\\
   \tag{3}
\end{pmatrix}

となります。
$x_3,...x_n$の列でも同様に計算していくと、

\begin{pmatrix} 
  x_1 &= & b_1^{(n)},x_2 &= & b_2^{(n)}\\
   &  & \dots,  \\
   &&,x_n&=&b_n^{(n)}\\
   \tag{4}
\end{pmatrix}

となり、解が求まります。
ところで、解を求めるためにはそれぞれの係数の計算方法を公式化する必要があります。

a_{12}^{(1)}=a_{12}/a_{11},a_{13}^{(1)}=a_{13}/a_{11},...,a_{1j}^{(1)}=a_{1j}/a_{11} \ (j=2,3,...,n)

a_{23}^{(2)}=a_{23}^{(1)}/a_{22},a_{24}^{(2)}=a_{24}^{(1)}/a_{22},...,a_{2j}^{(2)}=a_{2j}^{(1)}/a_{22} \ (j=3,4,...,n)

a_{34}^{(3)}=a_{34}^{(2)}/a_{33},a_{34}^{(3)}=a_{34}^{(3)}/a_{33},...,a_{3j}^{(3)}=a_{3j}^{(2)}/a_{33} \ (j=4,5,...,n)

となり、また、

a_{13}^{(2)}=a_{13}^{(1)}-a_{23}^{(2)}\cdot a_{12}^{(1)},
a_{14}^{(2)}=a_{14}^{(1)}-a_{24}^{(2)}\cdot a_{13}^{(1)}
,...,a_{1j}^{(2)}=a_{1j}^{(1)}\cdot a_{2j} \ (j=3,...,n)

a_{24}^{(3)}=a_{24}^{(2)}-a_{34}^{(3)}\cdot a_{23}^{(2)},
a_{25}^{(3)}=a_{25}^{(2)}-a_{35}^{(3)}\cdot a_{24}^{(2)}
,...,a_{2j}^{(3)}=a_{2j}^{(2)}\cdot a_{3j} \ (j=4,...,n)

$b$も同様に計算すると、最終的には以下のような形で要素を表せます。

\begin{pmatrix} 
a_{kj}^{(k)}=a_{kj}^{(k-1)}/a_{kk}^{(k-1)} \\
b_{k}^{(k)}=b_{k}^{(k-1)}/a_{kk}^{(k-1)}
\end{pmatrix}
\begin{pmatrix}
j=k+1,...,n\\
k=1,2,...,n
\end{pmatrix}

\begin{pmatrix}
a_{ij}=a_{ij}^{(k-1)}-a_{ik}^{(k-1)}*a_{kj}^{(kj)} \\
b_{i}=b_{i}^{(k-1)}-a_{ik}^{(k-1)}*b_k^{(k)}\\
\end{pmatrix}
\begin{pmatrix}
i=1,...,k-1,k+1,...,n\\
j=k+1,...,n\\
k=1,2,...,n-1
\end{pmatrix}

$なお、b_nの計算時に限り、k=nの値を取るため、jが登場するa_{kj}^{(k)}の計算時、kはnの値を取らないものとします。$
上記の式をプログラムコードに使用し、解を求めます。
ところがこの計算をする際には以下の2点の問題があります。
(1)軸要素(今回の場合は$a_{ii}$($i$=1,2,...n))が0に近いと、精度が落ちる可能性がある。
(2)係数の絶対値が極端に違っていると、情報落ち誤差が発生しやすく、演算が不安定になる場合がある。
なお、情報落ち誤差というのは、絶対値が大きく異なる値を足し引きすることで、小さい値の情報が無視されてしまうことを表します。

プログラムコードでは、これらの問題を解決しなければなりません。
まず軸要素が0に近いと精度が落ちる、という問題を解決する方法は、方程式の中から最大値を係数として$a_{kk}$に持ってきて計算することで解決を図ります。
なお、最大値も係数として0に近すぎる場合は計算できないので、例外処理をプログラム上に作成します。
次に、(2)の問題は、行列を正規化(スケーリング)することで解決を図ります。
ここでいう正規化(スケーリング)とは、各行列の最大値でそれぞれの値を割り、係数の絶対値がほぼ同じ値(0~1)になるように調整することをあらわしています。
つまり、係数の絶対値の差を減らせるので、(2)の情報落ち誤差の発生率を減らすことができます。
下記の手順で正規化を行います。

①与えられた連立一次方程式
\left(\sum_{j=1}^n a_{ij}x_j\right) =
b_i\left(i=1,2,...,n\right)　を

\sum_{j=1}^n \frac{a_{ij}}{s_i}x_{j}=\frac{b_i}{s_i} \
\left(ただし、s_i=\max|a_{ij}|\ \ 1\leq j \leq n \right)

と変換する
②①で得た方程式を

\sum_{j=1}^n \frac{a_{ij}}{s_i t_j}y_{j}=\frac{b_i}{s_i} \
\left(ただし、t_j=\max|\frac{a_{ij}}{s_i}|\ \ 1\leq j \leq n \right)

と変換する
③②の式を解けば,解$x_j$は、

x_{j}=\frac{y_j}{t_j} \

となる。
以上のことをふまえてプログラムを書くと、

#include <stdio.h>
#include <math.h>
#define eps 1.e-5 //マクロの定義。epsとコードに書くと、1.0*10^(-5)に変換されるということをあらわしている。行列の最大値がこれを下回った場合は例外処理する.
void sweep(double *,double *,int,int,int *,int *); 

int main(void)
{
    double  a[4][5],t[4]; //aは行列で、tは正規化用の配列.
    int iwork[4]; //作業領域iworkを用意する.
    int m,n,ILL=0; //ILLは例外処理に使用.
    scanf("%d %d",&n,&m); //nは行数,mは列数を表している.
    for(int i=0;i<n;i++){
        for(int j=0;j<m;j++){
            scanf("%lf",&a[i][j]);
        }
    }

    sweep((double *)a,t,n,m,iwork,&ILL); //正規化と解の計算.なお2次元配列aにおいて、今回はシングルポインタ型で値を渡す.
    if(ILL!=1) //例外処理.
    {
        for(int i=0;i<n;i++)
        {
            for(int j=n;j<m;j++)
            {
                printf("x(%d) =%8.4f\n",i+1,a[i][j]); //解の出力.
            }
        }
    }
 return 0;
}

void sweep(double *a,double *t,int n,int m, int *iwork,int *ILL)
{
double max,w;
int i,j,k,iw,p,q;

for(i=0;i<n;i++){  //ここから正規化の処理を行う.
    max=fabs(*(a+m*i));
    for(j=0;j<n;j++){
        if(max<fabs(*(a+m*i+j))){
            max=fabs(*(a+m*i+j));
        }    
    }
    for(j=0;j<m;j++) *(a+m*i+j)=*(a+m*i+j)/max; //ここでは各行の最大値で各行を割り算している.正規化第一段階.
}
    for(j=0;j<n;j++)
    {
        max=fabs(*(a+j));
        for(i=0;i<n;i++)
        {
        if(max<fabs(*(a+m*i+j))){
            max=fabs(*(a+m*i+j));
        }
        }
        *(t+j)=max;
        for(i=0;i<n;i++) *(a+m*i+j)=*(a+m*i+j)/max; //ここでは各列の最大値で、各列の割り算を行っている.正規化第二段階.
    }
    for(i=0;i<n;i++)  //作業領域iworkの初期化.
    {
        *(iwork+i)=i;
    }
    for(k=0;k<n;k++)  //掃き出し法の第k段階.
    {
        max=fabs(*(a+m*k+k)); //最大要素を持つ行列と、現在計算している行列の場所を入れ替える.
        p=k;
        q=k;
        for(j=k;j<n;j++)
        {
            for(i=k;i<n;i++)
            {
                if(max<fabs(*(a+m*i+j)))
                {
                    max=fabs(*(a+m*i+j));
                    p=i;
                    q=j;
                }
            }
        }
       if(max<=eps) //例外処理.
       {
            *ILL=1;
            printf("Matrix is ill.\n");
            return;
        }
        for(i=0;i<n;i++)  //列の入れ替え.
        {
            w=*(a+m*i+k);
            *(a+m*i+k)=*(a+m*i+q);
            *(a+m*i+q)=w;
        }
        for(j=k;j<m;j++) //行の入れ替え.
        {
            w=*(a+m*k+j);
            *(a+m*k+j)=*(a+m*p+j);
            *(a+m*p+j)=w;
        }
        i=*(iwork+k);  //入れ替え場所を記録.
        *(iwork+k)=*(iwork+q);
        *(iwork+q)=i;

        for(j=k+1;j<m;j++)
        {
            *(a+m*k+j)=*(a+m*k+j)/(*(a+m*k+k)); //a_{kj}^{(k)}の計算.
        }
        for(i=0;i<n;i++)
        {
            if(i!=k)
            {
                for(j=k+1;j<m;j++)
                {
                    *(a+m*i+j)=*(a+m*i+j)-*(a+m*i+k)*(*(a+m*k+j)); //a_{ij}^{(k)}の計算.
                }
            }
        }
    }
    for(j=n;j<m;j++)         //解の整列.
    {
        for(i=0;i<n;i++)
        {
            iw=*(iwork + i);
            *(a+m*iw+n-1)=*(a+m*i+j); //n-1列を作業領域として確保.
        }
        for(i=0;i<n;i++)
        {
            *(a+m*i+j)=*(a+m*i+n-1)/(*(t+i)); //作業領域から戻し、x_jの計算.
        }
    }
    return ;
}

となります。
試しに、以下の方程式を計算します。

\begin{pmatrix} 
x_1&+&3x_2&+&7x_2&+&2x_3&=&4\\
4x_1&+&x_2&+&9x_3&+&x_4&=&16\\
2x_1&+&2x_2&+&3x_3&+&x_4&=&3 \\
x_1&+&x_2&+&x_3&+&5x_4&=&6 
\end{pmatrix}

入力

実行結果

x(1) =  2.0000
x(2) = -3.0000
x(3) =  1.0000
x(4) =  2.0000

このように、Guass-jordan法を使用すれば、連立一次関数の計算をプログラムで行うことができます。
参考文献
鈴木淳子 (2001) 『C言語と数値解析法』株式会社培風館