#include <memory.h>
// LU分解のみを行う 結果を am に統合して出力
void LU44split(double * am/* 4*4 */, double * m/* 4*4 */)
{
　int　i,j,k;
　memcpy(m, am, sizeof(double)*16);
　for(i=0; i<4-1; i++) for(j=i+1; j<4; j++){
　　am[j*4+i]=am[j*4+i]/am[i*4+i];
　　for(k=i+1; k<4; k++) am[j*4+k]=am[j*4+k]-am[j*4+i]*am[i*4+k];
　}
}
//LU分解され統合された行列 bm と右辺列ベクトルである cm から連立一次方程式を解く
void LU44solve(double * am/* 1*4 */, double * bm/* 4*4 */, double * cm/* 1*4 */)
{
　int　i, j;
　memcpy(am, cm, sizeof(double)*4);
　for(i=1; i<4; i++) for(j=0; j<i; j++){// L*y=b: calc y
　　am[i]=am[i]-bm[i*4+j]*am[j];
　}
　for(i=r-1; i>=0; i--){　// U*x=y: calc x
　　for(j=4-1; j>i; j--) am[i]=am[i]-bm[i*4+j]*am[j];
　　am[i]=am[i]/bm[i*4+j];
　}
}
多分これでいい(確認はしてない）。まえ〜に作った任意の行列用のものを適当に4*4にしたから
間違ってる可能性も高い。そこは何とかしてくれ。ちなみに上位三角行列と
下位三角行列(対角要素が1になるようにしている)で後者の対角要素が1なもんで
これを一つにして am に出力している。ってことでL*U=Aの確認とかは適当にやってね。