[線形代数-ｃ言語]ガウスジョルダン法を使った連立一次方程式のコードを改造してほしい

Question

解決したいこと

以下の式からなる連立一次方程式をガウスジョルダン法を用いてc言語で実装したい

2x + 4y + 1z - 3w = 0
-1x - 2y + 2z + 4w = 10
4x + 2y - 3z + 5w = 2
5x - 4y - 3z + 1w = 6

やったこと

以下のようにソースコードと結果を示します。

#include <stdio.h>
#define N 4


void gauss_jordan(double a[N][N+1], int size);

int main(void)
{

    double coeff[][N+1] = {
                    {2,4,1,-3,0},
                    {-1,-2,2,4,10},
                    {4,2,-3,5,2},
                    {5,-4,-3,1,6},
                            };
    gauss_jordan(coeff, N);

    // show answer
    printf("Answer\n");
    for (int i=0; i < N; i++) {
        printf("x%d = %.15f\n",i, coeff[i][N]);
    }

    return 0;
}


void gauss_jordan(double a[N][N+1], int size) 
{
    // show default
    printf("Default:\n");
    for (int i=0;i<N;i++) {
        for (int j=0;j<=N;j++){
            printf("%10.5f",a[i][j]);
        }printf("\n");
    }printf("\n");

    for (int k=0; k < size; k++){
        double pivot = a[k][k];
        //　make 1 for pivot
        for (int i =k; i < size+1; i++) {
            a[k][i] *= 1.0/pivot;

        }

        // substract 1*s(=a[k][i]*s) and make 0 
        for (int j=0; j< size; j++) {
            if (j == k) continue;
            double s = a[j][k];
            for (int p=k; p < size+1; p++) {
                a[j][p] = a[j][p] - s*a[k][p];
            }
        }
        // show a[][]
        printf("Loop：k=%d\n",k);
        for (int i=0;i<N;i++) {
            for (int j=0; j<=N;j++){
                printf(" %10.5f",a[i][j]);
            }printf("\n");
        }printf("\n");

    }
    return;
}

結果

Default:
   2.00000   4.00000   1.00000  -3.00000   0.00000
  -1.00000  -2.00000   2.00000   4.00000  10.00000
   4.00000   2.00000  -3.00000   5.00000   2.00000
   5.00000  -4.00000  -3.00000   1.00000   6.00000

Loop：k=0
    1.00000    2.00000    0.50000   -1.50000    0.00000
    0.00000    0.00000    2.50000    2.50000   10.00000
    0.00000   -6.00000   -5.00000   11.00000    2.00000
    0.00000  -14.00000   -5.50000    8.50000    6.00000

Loop：k=1
    1.00000       -nan       -inf       -inf       -inf
    0.00000       -nan        inf        inf        inf
    0.00000       -nan        inf        inf        inf
    0.00000       -nan        inf        inf        inf

Loop：k=2
    1.00000       -nan       -nan       -nan       -nan
    0.00000       -nan       -nan       -nan       -nan
    0.00000       -nan       -nan       -nan       -nan
    0.00000       -nan       -nan       -nan       -nan

Loop：k=3
    1.00000       -nan       -nan       -nan       -nan
    0.00000       -nan       -nan       -nan       -nan
    0.00000       -nan       -nan       -nan       -nan
    0.00000       -nan       -nan       -nan       -nan

Answer
x0 = -nan
x1 = -nan
x2 = -nan
x3 = -nan

一応0がない時なら正常に解を導き出してくれます。
これだと、0が解にあるためか、結果が-nanになってしまっておかしくなってしまいます。
どこをどのように改造したら0を考慮する必要がなくなるのかわかりません。

考え

対角化のところでpivot=0が出てきてしまった時に、行数を入れ替える処理がかけたらいいんですが....
わかりません

お願いします。

上記の内容から0があるときにでも解を導けるようにどのようにそーすこーどをいじればよいのかアドバイスをお願いします。
自分の考えとして間違いないのですが、二回目の対角化のところで1ではなく0が出てきてしまったので、その時に、行を入れ替える処理を書きたいんですが、思いつきません....

初心者ですので、至らない点がありましたら指摘してくださいませ。

Answer 1

久しぶりに、大学のころを思い出しました。
部分ピポッテイングとか、完全ピボッティングがあるのですね。
勉強になりました。
以下を参考にして、部分ピポッテイングでコードを入れてみました。
http://nkl.cc.u-tokyo.ac.jp/12s/3D/pivoting.pdf
参考になれば幸いです。
行を検索して、最大値があれば、検索開始行と最大値が見つかった行を入れ替えています。

// 部分ピボッティング
void partial_pivot(double a[N][N + 1], int pos, int size) {
	double max = a[pos][pos];
	int max_line = pos;
	double tmp;

	// search max value in lines
	for (int i = pos + 1; i < size; i++) {
		if (a[i][pos] > max) {
			max = a[i][pos];
			max_line = i;
		}
	}
	if (max_line != pos) {
		// max_lineとposを入れ替える
		for (int i = 0; i < size + 1; i++) {
			tmp = a[pos][i];
			a[pos][i] = a[max_line][i];
			a[max_line][i] = tmp;
		}
	}
}

以下に上記の関数を追加

	for (int k = 0; k < size; k++) {
		partial_pivot(a, k, size);  //★ここ
		double pivot = a[k][k];
		//　make 1 for pivot
		for (int i = k; i < size + 1; i++) {

このソースでのanswer

Answer
x0 = 2.000000000000000
x1 = -1.000000000000000
x2 = 3.000000000000000
x3 = 1.000000000000000