次のようにある行列Aを下三角行列Lと上三角行列Uに分解することを考えます．
これは例えば行列Aが3×3の行列であれば次のようになります．
このように分解することのメリットは，連立1次方程式を解くときを考えると， 下三角行列と上三角行列に分解していれば後は簡単に解くことができるという事です．
のような連立1次方程式を解くことを考えてみます．LU分解を行うと，
とでき，ベクトルcを用いて，
とおくと，
となります．ここで例えば３次であれば，
のようになり，
を解くことと同じ事になります．するとベクトルcはすぐにわかりますので， 今度は，
より，
を解けばベクトルxが求められます．これは先程と同じで，
のようになりますので，ベクトルxが簡単に求まります．
このように，行列を下三角行列と上三角行列に分解できれば連立方程式が簡単に解けます．
行列Aを下三角行列Lと上三角行列Uに分解すれば連立方程式が簡単に解けることはわかりましたが， どうやったらそのように分解できるのでしょうか．とりあえず，
を計算してみる事にします．
どうやら対角成分より下と，それ以外で場合分けができそうです． もちろん，３次のときにどうなったからといって，数学的にそれが成り立つ事が証明されるわけではありませんが， ある程度自明なところがあります．
一般的に行列のかけ算を考えると，
となりますが，ここで対角成分よりも下側とそれ以外で分けると，
であり，行列Lの対角成分が１であることを考えると，
とできるため，行列L，行列Uの成分は，
となります．
以上で一般的にプログラムが出来そうですが，よく見ると下三角行列の上の部分は無駄ですし， 上三角行列の下の部分も無駄です． しかも行列Lの対角成分は１ですのでここも無くても 問題ありません．そのため行列Lと行列Uの２つの行列をまとめて，
のように値を保持します．それでは３次のときの例で流れを見てみます．
今この行列は元の行列Aを表しています．
このような順で計算すると，
このようになり，この右下２×２の小行列を考えると先程の３×３のときと 同じ状態になっています．そこでここまでを１セットと考え，同じ様に
と計算すると，
となってLU分解が完成します．
上のような流れでLU分解が進んでいきますが，それぞれの１セットの最初に 行の交換をする必要がある場合があります．というのは上の流れを見てもわかるように， そのときの対角成分の左上の値でその列の残りの部分を割っています． そのため割る値がゼロであれば困りますし，ゼロでなくてもゼロに近ければ 誤差も大きくなります．そこで毎回１セットが始まるときに，例えば上の流れの 最初の場合であれば，３行全ての行の内，１列目の要素の絶対値が最大であるものと 交換し，２セット目であれば２行目以降の中から２列目の値の絶対値が最大であるものと 交換します．もちろん交換が行われれば，どことどこが交換されたかが わかるようにどこかにその情報を保持しておかねばなりません．
width = m->width;
height = m->height;
if((order = (int*)malloc(sizeof(int)*height)) == NULL){
fprintf(stderr, "Allocation Error\n");
exit(1);
}
for(i=0; i<height; i++) order[i] = i;
動的に整数の配列を作成し，orderにポインタを保持しておきます． その後0,1,2・・・と初期化していますが， これは行を交換したときにこの値も同じく交換しておくことで， 例えば２行目と３行目を交換すれば，０，２，１・・・となることになります．
d = i;
max = fabs(m->a[i*width + i]);
for(j=i+1; j<height; j++){
if(max < fabs(m->a[j*width + i])){
d = j;
max = fabs(m->a[j*width + i]);
}
}
if(d!=i){
temp = order[i];
order[i] = order[d];
order[d] = temp;
ExchangeRowvector(m, i, d);
}
if(!(m->a[i*width + i])) continue;
for文までは現在注目している列の中から絶対値が最大の行を探しています． そして最大の行が現在の行でなければ現在の行と最大の行を入れ換えます． ExchangeRowvectorはi行目とd行目を交換する関数です． 交換しても対角成分がゼロであればその回はLU分解が終わっているといえるので， continueします．
肝心のLU分解の心臓部ですが，以上のような部分がごちゃごちゃしているだけで， 上の肝じゃない部分を除くと次のようになります．
for(i=0; i<height-1; i++){
for(j=i+1; j<height; j++){
m->a[j*width + i] /= m->a[i*width + i];
for(k=i+1; k<width; k++){
m->a[j*width + k] -= m->a[j*width + i]*m->a[i*width + k];
}
}
}
長々と書いて来ましたが，実際のLU分解の心臓部はたったのこれだけです．
