C/C++ 语言实现矩阵LU分解

应用计算方法C语言程序：01

本博客作为后续博客使用C/C++ 语言实现使用LU分解求解线性方程组 的一部分，讲解使用C/C++ 语言实现矩阵LU分解。方程组中的LU分解，也叫做Doolittle分解。

后接博文：应用计算方法C语言程序：02：C/C++ 语言实现使用LU分解求解线性方程组 。

以下图例子进行说明：Ax=b；
double A[3][3] = { {1,3,3},{2,1,1},{2,3,4} };
double b[3] = { 1,2,1 };

Ax=b,A为方阵时，对矩阵A的LU分解矩阵L、U公式如下图：

对矩阵A的LU分解C代码如下：

#include <iostream>
#include "math.h"

using namespace std;

double L[3][3] = { 0 }, U[3][3] = { 0 };
void Doolittle(double a[3][3])
{
	for (int i = 0; i < 3; i++)
	{
		//更新L矩阵  下三角矩阵
		for (int j = 0; j <= i; j++)
		{
			if (j == i) L[i][i] = 1;
			else
			{
				L[i][j] = a[i][j];
				for (int k = 0; k < j; k++)
				{
					L[i][j] -= L[i][k] * U[k][j];
				}
				L[i][j] /= U[j][j];
			}
		}

		//更新U矩阵  上三角矩阵
		for (int j = i; j < 3; j++)
		{
			U[i][j] = a[i][j];
			for (int k = 0; k < i; k++)
			{
				U[i][j] -= L[i][k] * U[k][j];
			}
		}
	}
}

int main()
{
	double A[3][3] = { {1,3,3},{2,1,1},{2,3,4} };
	double b[3] = { 1,2,1 };

	Doolittle(A);

	cout << "LU分解" << endl << endl;
	cout << "U = " << endl;
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			cout << U[i][j] << " ";
		}
		cout << endl;
	}
	cout << endl;
	cout << "L = " << endl;
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			cout << L[i][j] << " ";
		}
		cout << endl;
	}
	cout << endl;

	return 0;
}

运行结果见下图：