LU分解法のピボット選択機能実装について（C言語・gcc-9）

C言語の質問です。
講義の最終テストが終わったので講義で飛ばされた「LU分解法にピボット選択機能をつけたプログラム」を実装しようと思い試行錯誤したのですが、どうにもうまく実装できません。とはいえまだ他の講義のテストはまだ全然終わっていないのであまり時間をかけたくないのが正直なところです。どなたか、説明なしの答えだけでも良いので教えていただけませんか？
参考までに以下に、ピボット機能なしのプログラムコードを貼っておきます。
*連立方程式を
「aw+bx+cy+dz=e
fw+gx+hy+iz=j
kw+lx+my+nz=o
pw+qx+ry+sz=t」
と置くと、標準入力として
a b c d
f g h i
k l m n
p q r s
e j o t
と入力し、答えを表示させるコードになります。(見にくくてすみません。講義で取り扱ったピボット選択機能なし版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];
int i, j, k;
//⼊⼒データの受け取り
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++){
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]);
}

No.1 ベストアンサー 回答者： Ogre7077 回答日時： 2022/07/24 01:48

行だけ入れ替える、部分ピボット選択だけならば

LU 分解前に以下の処理追加

for (i=0; i<M; i++) {
_ pivot = 絶対数が最大である行を探す(a, M, i);
_ // 検索対象 : a[i][i] から a[i][M-1] の範囲
_
_ if (i != pivot) {
_ _ 行を入れ替える(a, M, i, pivot); // i行とpivot行を入れ替え
_ _ 値を入れ替える(b, M, i, pivot); // ベクトル B も同様
_ }
}

行と列を入れ替える、完全ピボット選択を目指すなら

for (i=0; i<M; i++) {
_ pivot = 絶対数が最大である行と列を探す(a, M, i);
_ // 検索対象 : a[i][i] から a[M-1][M-1] の範囲
_
_ if (i != pivot.row) {
_ _ 行を入れ替える(a, M, i, pivot.row); // i行とpivot行を入れ替え
_ _ 値を入れ替える(b, M, i, pivot.row); // ベクトル B も同様
_ }
_ if (i != pivot.col) {
_ _ 列を入れ替える(a, M, i, pivot.col);
_ }
_ ベクトルxの並び順[i] = pivot.col; // 計算後にxを並び替える
}

この回答へのお礼

ありがとうございます！おかげで解決できました。

お礼日時：2022/07/25 10:20