ピボット選択無し．3次の行列にしか対応してない．
元の行列にLU行列を重ねて上書き．
けっこうだるい．
// 数値解析-LU分解.cpp :
// 3次の行列のLU分解．ピボット選択なし
#include "stdafx.h"
#include "iostream"
#include "cstdio"
using namespace std;
#define N 3 //3次のLU分解しか対応せず
void output(void);
void LUdiv(void);
void Ly_b(void);
void Ux_y(void);
double mat[N][N] = { { 4.0, 1.0, 1.0 }, //簡単のためグローバル変数
{ 1.0, 3.0, 1.0 },
{ 2.0, 1.0, 5.0 }
};
double b[N] = {9.0, 10.0, 19.0};
double y[N],x[N]; //解を格納する
int main(int argc, _TCHAR* argv[])
{
output();
LUdiv();
Ly_b();
Ux_y();
output();
printf("x={ %1.3lf, %1.3lf, %1.3lf }\n",x[0],x[1],x[2]); //解の表示
return 0;
}
void LUdiv(void){ //LU分解
double tmp = 0.0,tmp2 = 0.0;
mat[0][1] = mat[0][1] / mat[0][0];
mat[0][2] = mat[0][2] / mat[0][0];
for(int i = 1; i < N;i++){
for(int k = 0;k <= i-1;k++){
tmp += mat[i][k] * mat[k][i];
}
mat[i][i] = mat[i][i] - tmp;
if(i != N-1){ //最後は実行しない
for(int k = 0;k <= i-1;k++){
tmp2 += mat[i+1][k] * mat[k][i];
}
mat[i+1][i] = mat[i+1][i] - tmp2;
tmp2 = tmp = 0.0;
for(int k = 0;k <= i-1;k++){
tmp += mat[i][k] * mat[k][i+1];
}
mat[i][i+1] = (mat[i][i+1] - tmp) / mat[i][i];
tmp = 0.0;
}
}
}
void Ly_b(void){ //Ly = bを満たすyを求める
double s = 0.0;
for(int i = 0;i < N;i++){
for(int j = 0;j < i;j++){
s += mat[i][j] * y[j];
}
y[i] = (b[i] - s) / mat[i][i];
s = 0.0;
}
}
void Ux_y(void){ //Ux_yを満たすxを求める
double s = 0.0;
for(int i = N-1;i >= 0;i--){
x[i] = y[i];
for(int j = N-1;j > i;j--){
s += mat[0][j] * x[j];
}
x[i] -= s;
s = 0.0;
}
}
void output(void){
for(int i = 0;i < N;i++){
for(int j = 0;j < N;j++){
printf("%1.3lf ",mat[i][j]);
}
printf("\n");
}
printf("\n");
}
