逆行列の計算

[Serdol]

http://kazuhirokazu.org/2014/01/13/c%E8 ... %E3%82%8B/

上記サイトをほぼそのまま使って逆行列の計算を行おうとしたところ、
{2, 4, -2, -4, 1, 0, 0, 0},
{-6, -12, 12, 24, 0, 1, 0, 0},
{4, 2, 2, -4, 0, 0, 1, 0},
{2, -4, -2, 4, 0, 0, 0, 1}};

でうまく逆行列が出力されません。（サイトの数値をそのまま使った場合はうまくいきました）
以下のように出力されます。
1.000000 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00
-1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00
0.000000 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00
-1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00 -1.#IND00

出力を付け足した以外は特に何も変えていないのですが、なぜか入力例が変わった途端に上手くいきません

原因を調べたところpivotの値がループの2周目からおかしいということはわかりましたが、なぜなのかわかりません。
よろしくお願いします。

ほとんどサイトのままですが、一応ソースを載せておきます

コード:

#include <stdio.h>
#define N 4
int main(void)
{
	
	double M[N][N * 2] ={
    {2, 4,  -2, -4,   1, 0, 0, 0},
    {-6, -12, 12, 24,  0, 1, 0, 0},
    {4, 2, 2, -4,   0, 0, 1, 0},
    {2, -4, -2, 4,    0, 0, 0, 1}};
 
double mul;
double pivot;
int i,j,k,n;
// 対角成分が1で正規化された階段行列を作る
for (i = 0; i < N; ++i)
{
    // 対角成分の選択、この値で行成分を正規化
    pivot = M[i][i];
    for (j = 0; j < N * 2; ++j)
    {
        M[i][j] = (1 / pivot) * M[i][j];
    }
 
    // 階段行列を作る為に、現在の行より下の行につ
    // i列目の成分が0になるような基本変形をする
    for (k = i + 1; k < N; ++k)
    {
        mul = M[k][i];
        for (n = i; n < N * 2; ++n)
        {
            M[k][n] = M[k][n] -  mul * M[i][n];
        }
    }
}
 
// 下から上に向かって変数に代入して、左側について
for (i = N - 1; i > 0; --i)
{
    for (k = i - 1; k >= 0; --k)
    {
        mul = M[k][i];
        for (n = i; n < N * 2; ++n)
        {
            M[k][n] = M[k][n] - mul * M[i][n];
        }
    }
}
        for( i = 0; i<N;i++) {
        		for(j=0;j<2*N;j++) {
        			printf(" %f",M[i][j]);
        		}
        		printf("\n");
        	}
	return 0;
}

[みけCAT]

「対角成分が1で正規化された階段行列を作る」の最初のループが終わった時点で、変形の結果M[1][1]が0になってしまい、
それをそのままピポットにしてしまうのでうまくいかないですね。

ピポットを0にしないために、例えば「今見ている行とそれより下の行の中からピポットになる列の数の絶対値が一番大きい行を選び、今の行と入れ替えてから変形を行う」などの工夫をすると改善できるでしょう。

コード:

#include <stdio.h>
#include <math.h> /* fabs()用 */
#define N 4
int main(void)
{

    double M[N][N * 2] ={
        {2, 4,  -2, -4,    1, 0, 0, 0},
        {-6, -12, 12, 24,  0, 1, 0, 0},
        {4, 2, 2, -4,      0, 0, 1, 0},
        {2, -4, -2, 4,     0, 0, 0, 1}};
 
    double mul;
    double pivot;
    int i,j,k,n;
    // 対角成分が1で正規化された階段行列を作る
    for (i = 0; i < N; ++i)
    {
        // ピポットの候補の中で絶対値が一番大きいものを選ぶ
        int k = i;
        for (j = i + 1; j < N; ++j)
        {
            if (fabs(M[j][i]) > fabs(M[k][i])) k = j;
        }
        // 絶対値が一番大きい候補をピポットにするため、行を入れ替える
        if (k != i)
        {
            for (j = 0; j < N * 2; ++j)
            {
                double temp = M[i][j];
                M[i][j] = M[k][j];
                M[k][j] = temp;
            }
        }
        // 対角成分の選択、この値で行成分を正規化
        pivot = M[i][i];
        for (j = 0; j < N * 2; ++j)
        {
            M[i][j] = (1 / pivot) * M[i][j];
        }
     
        // 階段行列を作る為に、現在の行より下の行につ
        // i列目の成分が0になるような基本変形をする
        for (k = i + 1; k < N; ++k)
        {
            mul = M[k][i];
            for (n = i; n < N * 2; ++n)
            {
                M[k][n] = M[k][n] -  mul * M[i][n];
            }
        }
    }
     
    // 下から上に向かって変数に代入して、左側について
    for (i = N - 1; i > 0; --i)
    {
        for (k = i - 1; k >= 0; --k)
        {
            mul = M[k][i];
            for (n = i; n < N * 2; ++n)
            {
                M[k][n] = M[k][n] - mul * M[i][n];
            }
        }
    }
    for( i = 0; i<N;i++) {
        for(j=0;j<2*N;j++) {
            printf(" %f",M[i][j]);
        }
        printf("\n");
    }
    return 0;
}

出力

コード:

 1.000000 0.000000 0.000000 0.000000 0.142857 0.023810 0.142857 0.142857
 -0.000000 1.000000 0.000000 0.000000 0.428571 0.071429 -0.071429 -0.071429
 0.000000 0.000000 1.000000 0.000000 -0.107143 0.023810 0.142857 -0.107143
 0.000000 0.000000 0.000000 1.000000 0.303571 0.071429 -0.071429 0.053571

[Serdol]

なるほど、ピボット選択がきちんと行われていなかったのですね。ありがとうございます。

サイトの方もその辺のことは正直書いておいてほしかったです‥‥これでは困る人がたくさんいると思うので‥‥