<br>
D言語で作るLU分解 part2 ~PA=LU編~
前回はA = LUとなる一般的なLU分解を行いました。
プログラムも一般的な簡単なものになったかと思います。
前回のA = LUとなるLU分解では、たとえばA11=0のときには正常に計算できません。
(L21はU11で割っている。U11=A11であるのでL21はA11で割っている。)
そこで前回の予告通り、不完全ピボット選択となるPA = LU分解を説明したいと思います。
仕組みは簡単で
A = LU となる行列A,L,Uが存在するとすると、
両辺に置換行列Pを左からかけて
PA = PLU
となります。
PA = A'として
A' = PLU
となります。
つまり、Aを行について適当に置換させた行列A'のLU分解の結果はもちろんPLUとなります。
ではプログラムではどうなるでしょうか？
下は実際の似非フローチャート的なものです。
U[i][i]を生成
↓
if( U[i][i] == 0 ){
//i行とi以降の行を入れ替え
--i;
continue;
}
となります。
if文の中の行の入れ替えは行列AとPについて行うので、
A,Pがstruct Matrix!T であり、その行についての配列演算子が使えるなら
入れ替える行の番号をiとjとすると、
T[] temp = A[i];
A[i] = A[ｊ];
A[ｊ] = temp;
と書くだけです。
Pについても同様にかけます。
では前回のプログラムに上記を組み込んだものを載せておきます。
Matrix!(T)[] pludec(Matrix!T Mat){
int row = Mat.row;
int column = Mat.column;
//LU分解を行う
Matrix!T L = Matrix!T(row,column,0);
Matrix!T U = Matrix!T(row,column,0);
Matrix!T P = Matrix!T(row,column,0);
Matrix!T A = Matrix!T(row,column,0);
//Lの対角成分をすべて1に
for(int i=0;i<row && i<column;i++)
L[i][i] = 1;
if(Mat[0][0] == 0){
//置換を行う
for(int i=1;i<row;i++){
if(Mat[i][0] != 0){//入れ替え
for(int j=0;j<column;j++){
A[0][j] = Mat[i][j];
A[i][j] = Mat[0][j];
}
P[0][i] = 1;
P[i][0] = 1;
}else{
for(int j=0;j<column;j++)
A[i][j] = Mat[i][j];
P[i][i] = 1;
}
}
}else{
for(int i=0;i<row;i++){
P[i][i] = 1;
for(int j=0;j<column;j++)
A[i][j] = Mat[i][j];
}
}
//まずUの一列目の計算を行う
for(int j=0;j<column;j++)
U[0][j] = A[0][j];
//1行目のLについて計算を行う
for(int i=1;i<row;i++)
L[i][0] = A[i][0] / U[0][0];
//1列目(行目)以降のL,Uについて計算する
for(int i=1;i<row;i++){
//Uについて
for(int j=i;j<column;j++){
U[i][j] = A[i][j];
//ここから引いていく
for(int k=0;k<i;k++)
U[i][j] -= L[i][k]*U[k][j];
}
//0なら入れ替え
if(U[i][i] == 0){
for(int l=i+1;l<row;l++)
if(Mat[l][0] != 0){//入れ替え
T[] temp;
temp = A[i];
A[i] = A[l];
A[l] = temp;
temp = P[i];
P[i] = P[l];
P[l] = temp;
}
--i;
continue;
}
//Lについて
for(int j=i+1;j<row;j++){
L[j][i] = A[j][i];
//ここから引いていく
for(int k=0;k<i;k++)
L[j][i] -= L[j][k]*U[k][i];
//最後に割ってやる
L[j][i] /= U[i][i];
}
}
}
ただし、この方法でもまだ不完全で、もし最終行でU[i][i]が0だった場合は正常に動作しません。
その場合はさかのぼってPA=LUで計算しなおしか、PAQ=LU分解をして、こちらで再計算かどちらかになると思います。
次回は上記のどちらかについて書きたいと思います。
