数値解析課題-LU分解

ピボット選択無し．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");
	
}