LU分解法のピボッティングについて(C言語/gcc-9)

講義の課題で4変数連立方程式をLU分解法で解くピボッティング機能付きプログラムを作成しており、一旦プログラムは書けたのですが、どうもうまく動作しません。講義中に先生に聞いて対処法を聞いてみたのですが、ざっと見た感じでは先生も見つけられなかったみたいで解決できず...そこでお願いなのですが、以下にプログラミングのコードをはっつけますので、どこがまずいか見つけていただけないでしょうか。また、対処法や変更案等あれば幸いです。どうかお願いします。
（コード）

#include <stdio.h>
#define M 4
int main(){
//変数宣⾔
double a[M][M];
double b[M];
double c[M];
double l[M][M];
double u[M][M];
double x[M];
double d[M];
int i, j, k;
double tmp, max;
//⼊⼒データの受け取り
for(i=0;i<M;i++){
for(j=0;j<M;j++) scanf("%lf",&a[i][j]);
}
for(i=0;i<M;i++){
scanf("%lf",&b[i]);
}
//L⾏列，U⾏列の初期化
for(i = 0; i < M; i++){ /* L行列，U行列を1と0で初期化 */
for(j = 0; j < M; j++){
u[i][j] = 0;
if(i == j)
l[i][j] = 1;
else
l[i][j] = 0;
}
}
//⼊⼒⾏列の出⼒(ピボッティングでa⾏列を操作するため事前に出⼒)

printf("⼊⼒⾏列\n");
for(i=0;i<M;i++){
for(j=0;j<M;j++){
printf("%10.5lf",a[i][j]);
}
printf("%10.5lf\n",b[i]);
}
//結果の出⼒
for(i = 0; i < M; i++){
/*a行列をピボッティング*/
d[0] = a[0][i];
d[1] = a[1][i];
d[2] = a[2][i];
d[3] = a[3][i];
max = d[0];
for(j = 0; j < M; j++){
for(j = 0; j < M; j++){
if(max < d[j])
max = d[j];
}

if(max == d[j])
for(k = 0;k < M; k++){
tmp = a[j][k];
a[j][k] = a[j][i];
a[j][i] = tmp;
}
}
}
printf("⼊⼒⾏列(ピボッティング後)\n");
for(i=0;i<M;i++){
for(j=0;j<M;j++){
printf("%10.5lf",a[i][j]);
}
printf("%10.5lf\n",b[i]);
}
for(i = 0; i < M; i++){
for(j = i; j < M; j++){
/* U行列の生成 */
u[i][j] = a[i][j];
for(k = 0; k < i; k++){
u[i][j] -= u[k][j] * l[i][k];
}
}
for(j = i + 1; j < M; j++){
/* L行列の生成 */
l[j][i] = a[j][i];
for(k = 0; k < i; k++){
l[j][i] -= u[k][i] * l[j][k];
}
l[j][i] /= u[i][i];
}
}
for(i = 0; i < M; i++){
/* c行列の生成 */
c[i] = b[i];
for(j = 0; j < i; j++){
c[i] -= l[i][j] * c[j];
}
}
for(i = M - 1; i >= 0; i--){
/* x行列の生成 */
x[i] = c[i];
for(j = M - 1; j > i; j--){
x[i] -= u[i][j] * x[j];
}
x[i] /= u[i][i];
}
//L⾏列
printf("\nL⾏列\n");
for(i=0;i<M;i++){
for(j=0;j<M;j++){
printf("%10.5lf",l[i][j]);
}
printf("\n");
}
//U⾏列
printf("\nU⾏列\n");
for(i=0;i<M;i++){
for(j=0;j<M;j++){
printf("%10.5lf",u[i][j]);
}
printf("\n");
}


//解の出⼒
printf("\n答え\n");
for(i=0;i<M;i++){
printf("%10.5lf\n",x[i]);
}
}

質問者からの補足コメント

• 

  回答いただきありがとうございます！
  やっぱりそこですよね...色々コードを改造して標準入力ぶちこんで色々試してみたんですが、ピボッティング機能を消したときだけうまくいったので。
  「i列の中で最大のやつ見つけてそれが含まれる行をi行目と交換する～」っていうメカニズム的には考え方はそこまで遠くないと思うので、あとはそのピボッティング機能とそれ以外のコードの組み合わせ方次第だとは思うのですが...

  No.2の回答に寄せられた補足コメントです。 補足日時：2022/07/12 21:31

No.3 回答者： Tacosan 回答日時： 2022/07/13 01:02

#2 で挙げた部分で外の for がまわしているのは

for(j = 0; j < M; j++){
if(max < d[j])
max = d[j];
}

if(max == d[j])
for(k = 0;k < M; k++){
tmp = a[j][k];
a[j][k] = a[j][i];
a[j][i] = tmp;
}
の部分. で, このコードでは上の for が抜けたときには j の値は M になってる. その j の値で次の
if(max == d[j])
に入るわけだけど, それでいいの?

あと「最大値」でいいのか (「絶対値最大」であるべき?) とか, 「どの j で「最大」なのか覚えておけばいいんじゃねとか.

j を変える for の中で j を変える for を使うってところはあきらかにおかしいんだけど, ね.

No.2 回答者： Tacosan 回答日時： 2022/07/12 21:21

あ,

for(j = 0; j < M; j++){
for(j = 0; j < M; j++){
if(max < d[j])
max = d[j];
}

if(max == d[j])
for(k = 0;k < M; k++){
tmp = a[j][k];
a[j][k] = a[j][i];
a[j][i] = tmp;
}
}
の部分がおかしい.

No.1 回答者： Tacosan 回答日時： 2022/07/12 02:23

正直にいってテストする気はないけど, 手計算の過程とこのプログラムでの計算過程とを比較すればいいんじゃないかな.