C言語 上三角行列の出力

Question

前提・実現したいこと

C言語で与えられた行列に対する上三角行列をを表示するプログラム作っています。
以下のプログラムでは結果が上手く表示されません。
誤っている個所を教えていただきたいです。

該当のソースコード

#include<stdio.h>
	
	int main(){
	int a[4][4]={{2,-1,6,3},{2,-2,4,2},{-1,3,-4,4},{3,-2,12,12}}; 
	int b[4][4];
	int buf, n=4, i,j,k; 

	for(i=0;i<n;i++){ 
		for(j=0;j<n;j++){ 
			if(i<j){ 
			buf=a[j][i]/a[i][i]; 
				for(k=0;k<n;k++){ 
					b[j][k]-=a[i][k]*buf; 
				} 
			} 
		}
	}
	
	for(i=0;i<n;++i){ 
    	for(j=0;j<n;++j){ 
      		printf("%d  ",b[i][j]); 
      			if(j==n-1) 
          			printf("\n"); 
    	} 
 	 } 
return 0;
}

発生している問題・エラーメッセージ

以下が現在の出力結果です
0  9052424  5651992  2
0  9052433  -6  21
18  9052426  9052452  2
16152324  1703764  1979215300  5651999

Answer 1

要はLU分解っしょ? コレ↓の上三角成分を取り出せばいいんよね?

/*
   N次正方行列AをLU分解する。
   結果でAを上書きする。
 */

template<typename T,size_t N>
bool lu_factorize(T (&A)[N][N]) {
  for ( int i = 0; i < N-1; ++i ) {
    for ( int j = i+1; j < N; ++j ) {
      if ( A[i][i] == T(0) ) { return false; }
      A[j][i] /= A[i][i];
      for ( int k = i+1; k < N; ++k ) {
        A[j][k] -= A[j][i] * A[i][k];
      }
    }
  }
  return true;
}

Answer 2

ガウスの消去法

アルゴリズムについては こちらのpdf がわかりやすいので参考にしてください。
ピボット選択は零割り防止の他、丸め誤差が大きくなる計算を防ぐ目的があります。 ピボット選択については、こちら を参考にしてください。

#include <math.h>
#include <stdio.h>

#define n 4

void print_matrix(double A[][n])
{
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j)
            printf("%f ", A[i][j]);
        printf("\n");
    }
}

int main()
{
    // clang-format off
    double A[4][4] = {
        { 2, -1,  6, 3},
        { 2, -2,  4, 2},
        {-1,  3, -4, 4},
        { 3, -2, 12, 12}};
    // clang-format on
    printf("original\n");
    print_matrix(A);

    for (int i = 0; i < n; ++i) {
        printf(
            "step %d"
            "\n----------------------------\n",
            i);

        // 1. <ピボット選択> k行i列の中で最大の要素を調べる。
        // ----------------------------------
        double max_val = fabs(A[i][i]);
        int max_row = i;
        for (int row = i + 1; row < n; ++row) {
            if (fabs(A[row][i]) > max_val) {
                max_val = fabs(A[row][i]);
                max_row = row;
            }
        }

        // 2. <ピボット選択> 現在の行と値が最も大きい行を入れ替える。
        // ----------------------------------
        if (i != max_row) {
            for (int col = i; col < n; ++col) {
                double tmp = A[max_row][col];
                A[max_row][col] = A[i][col];
                A[i][col] = tmp;
            }
            printf("swap rows: %d <--> %d\n", i, max_row);
            print_matrix(A);
        }

        // 3. <前進消去>
        // ----------------------------------
        for (int row = i + 1; row < n; ++row) {
            double tmp = A[row][i] / A[i][i];
            printf("---> A[%d][%d] / A[%d][%d] = %f\n", row, i, i, i, tmp);

            /* 第 i 行を -A[row][i] / A[i][i] 倍して、第 j 行に加える */
            A[row][i] = 0;  // A[row][i] -c * A[i][i] = 0 は確定しているので、
                            // 計算しなくてもよい。
            for (int col = i + 1; col < n; ++col)
                A[row][col] += -tmp * A[i][col];
            printf("---> eliminate: row %d + %f * row %d\n", row, tmp, i);
            print_matrix(A);
        }
    }

    // 結果
    printf("result\n");
    print_matrix(A);
}

original
2.000000 -1.000000 6.000000 3.000000 
2.000000 -2.000000 4.000000 2.000000 
-1.000000 3.000000 -4.000000 4.000000 
3.000000 -2.000000 12.000000 12.000000 
step 0
----------------------------
swap rows: 0 <--> 3
3.000000 -2.000000 12.000000 12.000000 
2.000000 -2.000000 4.000000 2.000000 
-1.000000 3.000000 -4.000000 4.000000 
2.000000 -1.000000 6.000000 3.000000 
---> A[1][0] / A[0][0] = 0.666667
---> eliminate: row 1 + 0.666667 * row 0
3.000000 -2.000000 12.000000 12.000000 
0.000000 -0.666667 -4.000000 -6.000000 
-1.000000 3.000000 -4.000000 4.000000 
2.000000 -1.000000 6.000000 3.000000 
---> A[2][0] / A[0][0] = -0.333333
---> eliminate: row 2 + -0.333333 * row 0
3.000000 -2.000000 12.000000 12.000000 
0.000000 -0.666667 -4.000000 -6.000000 
0.000000 2.333333 0.000000 8.000000 
2.000000 -1.000000 6.000000 3.000000 
---> A[3][0] / A[0][0] = 0.666667
---> eliminate: row 3 + 0.666667 * row 0
3.000000 -2.000000 12.000000 12.000000 
0.000000 -0.666667 -4.000000 -6.000000 
0.000000 2.333333 0.000000 8.000000 
0.000000 0.333333 -2.000000 -5.000000 
step 1
----------------------------
swap rows: 1 <--> 2
3.000000 -2.000000 12.000000 12.000000 
0.000000 2.333333 0.000000 8.000000 
0.000000 -0.666667 -4.000000 -6.000000 
0.000000 0.333333 -2.000000 -5.000000 
---> A[2][1] / A[1][1] = -0.285714
---> eliminate: row 2 + -0.285714 * row 1
3.000000 -2.000000 12.000000 12.000000 
0.000000 2.333333 0.000000 8.000000 
0.000000 0.000000 -4.000000 -3.714286 
0.000000 0.333333 -2.000000 -5.000000 
---> A[3][1] / A[1][1] = 0.142857
---> eliminate: row 3 + 0.142857 * row 1
3.000000 -2.000000 12.000000 12.000000 
0.000000 2.333333 0.000000 8.000000 
0.000000 0.000000 -4.000000 -3.714286 
0.000000 0.000000 -2.000000 -6.142857 
step 2
----------------------------
---> A[3][2] / A[2][2] = 0.500000
---> eliminate: row 3 + 0.500000 * row 2
3.000000 -2.000000 12.000000 12.000000 
0.000000 2.333333 0.000000 8.000000 
0.000000 0.000000 -4.000000 -3.714286 
0.000000 0.000000 0.000000 -4.285714 

result
3.000000 -2.000000 12.000000 12.000000 
0.000000 2.333333 0.000000 8.000000 
0.000000 0.000000 -4.000000 -3.714286 
0.000000 0.000000 0.000000 -4.285714 

Comments

[slushii]

上三角行列は対角成分よりも下側の成分が 0 である行列ではなかったでしょうか。

[tiitoi]

> 与えられた行列に対する上三角行列をを表示する

すいません。おっしゃる通りです。
だとすると、与えられた行列の上三角行列を表示するというのはどういう意味でしょうか？
行列の対角成分より下側に0を代入して上三角行列を作るということですか？

[slushii]

ガウスの消去法で各行の間で計算していって対角成分より下側に0を作っていく感じではないかと思っているのですが。。。

[tiitoi]

つまり、質問内容は「(ガウスの消去法の前進消去で)与えられてた行列を上三角行列にする方法について知りたい」ということでしょうか？

[slushii]

そうです！

[slushii]

説明不足で済みませんでした。

[tiitoi]

ガウスの消去法を追記しました。

[episteme]

そっちかよ...