LU分解

Posted on | In

LU分解的详解，请见下面维基百科的超链接。

维基百科LU详解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#include <iostream>
#include <iomanip>

using namespace std;

template<typename T>
void print_mat(T* matrix, int row, int col)
{
	for (int i = 0; i < row; ++i)
	{
		for (int j = 0; j < col; ++j)
			cout << std::left << setw(9) << matrix[i * col + j] << "      ";

		cout << endl;
	}
	cout << endl;
}

template<typename T>
void tri_decomposition(T* mat_A, T* mat_L, T* mat_U, int row, int col)
{
	/*
	*	外层循环代表循环次数
	*	内层循环是分别处理L和U矩阵
	*/
	col--;					//统一从下标为0开始
	row--;					//统一从下标为0开始
	double cum_val = 0.0;
	for (int r = 0; r <= row; ++r)    //统一从下标0开始,与课本公式上的变量名保持一致
	{
		for (int i = r; i <= row; ++i) // U矩阵
		{
			if (0 == r)
				mat_U[r * (col+1) + i] = mat_A[r * (col+1) + i];
			else
			{
				cum_val = 0.0;						//置为0
				for (int k = 0; k <= r - 1; ++k)  //累加
				{
					cum_val += mat_L[r * (col + 1) + k] * mat_U[k * (col + 1) + i];
					//cout << "debug ";
					//cout << "step 1 cum_val = " << mat_L[r * (col + 1) + k] << "  " << r * (col + 1) + k <<"  "<< mat_U[k * (col + 1) + i]<<"  "<<cum_val << endl;
				}
				mat_U[r * (col + 1) + i] = mat_A[r * (col + 1) + i] - cum_val;
		
			}
		}
		
		for (int i = r; i <= row; ++i)	// L矩阵
		{
			if (0 == r)
				mat_L[i * (col + 1) + r] = (mat_A[i * (col + 1) + r]) / (mat_U[0 * (col + 1) + r]);

			else
			{
				cum_val = 0.0;
				for (int k = 0; k <= r - 1; ++k)   //累加
				{
					cum_val += mat_L[i * (col + 1) + k] * mat_U[k * (col + 1) + r];
				}
					//cum_val += mat_L[i * (col + 1) + k] * mat_U[k * (col + 1) + r];
				mat_L[i * (col + 1) + r] = (mat_A[i * (col + 1) + r] - cum_val) / mat_U[r * (col + 1) + r];
			}
				
		}
	}
}



int main()
{
	double mat_A [][3] = { {1,2,3}, {4,5,6}, {3,-3,5} };
	double mat_L [][3] = { {0,0,0}, {0,0,0} ,{0,0,0} };
	double mat_U [][3] = { {0,0,0}, {0,0,0}, {0,0,0} };

	print_mat((double*)mat_A, 3, 3);
	print_mat((double*)mat_L, 3, 3);
	print_mat((double*)mat_U, 3, 3);
	 
	tri_decomposition((double*)mat_A,(double*)mat_L,(double*)mat_U,3,3);
	print_mat((double*)mat_A, 3, 3);
	print_mat((double*)mat_L, 3, 3);
	print_mat((double*)mat_U, 3, 3);
}