C++实现Cholesky分解

题目：
编制程序求解矩阵 A 的 Cholesky 分解，并用程序求解方程组 Ax=b，其中

代码实现：

#include <iostream>
#include <math.h>
#include <iomanip>
using namespace std;
//Cholesky 分解
void Cholesky(double a[4][4], double L[4][4], int m, int n) {
	for (int j = 0; j < n; j++)
	for (int i = j; i < m; i++) {
		if (i == j) {
			double sum = 0;
			for (int k = 0; k < j; k++) {
				sum = sum + L[j][k] * L[j][k];
}
			L[i][j] = sqrt(a[i][j] - sum);
}
		else {
				double sum = 0;
				for (int k = 0; k < j; k++) {
					sum = sum + L[i][k] * L[j][k];
}
				L[i][j] = (a[i][j] - sum) / L[j][j];
}
}
}
//初始化矩阵 L
void InitL(double L[4][4], int m, int n) {
	for (int i = 0; i < m; i++)
		for (int j = 0; j < n; j++) {
			if (i < j) {
				L[i][j] = 0;
}
}
}
//显示 L 和 LT 矩阵
void Display(double L[4][4], double LT[4][4], int m, int n) {
	cout << "Cholesky 分解的 L 矩阵为：" << endl;
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			cout << setprecision(8) << setw(12) << L[i][j] << "";
}
		cout << endl;
}
	cout << "Cholesky 分解的 LT 矩阵为：" << endl;
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			cout << setprecision(8) << setw(12) << L[j][i] << "";
			LT[i][j] = L[j][i];
}
		cout << endl;
}
}
//Ax=b 的解答 中间步骤求 y
double* SolveOne(double L[4][4], double b[4], int m, int n) {
	static double y[4];
	y[0] = b[0] / L[0][0];
	for (int i = 1; i < m; i++) {
		double sum = 0;
		for (int k = 0; k < i; k++) {
			sum = sum + L[i][k] * y[k];
}
		y[i] = (b[i] - sum) / L[i][i];
}
	cout << "解答的中间结果 y 为：" << endl;
	for (int i = 0; i < m; i++) {
		cout << setprecision(8) << setw(12) << y[i] << endl;
}
	cout << endl;
	return y;
}
//Ax=b 的解答 最终结果
void SolveTwo(double L[4][4], double y[4], int m, int n) {
	double x[4];
	x[m - 1] = y[m - 1] / L[n - 1][n - 1];
	for (int i = m - 2; i >= 0; i--) {
		double sum = 0;
		for (int k = i + 1; k < m; k++) {
			sum = sum + L[k][i] * x[k];
}
		x[i] = (y[i] - sum) / L[i][i];
}
	cout << "解答的最终结果 x 为：" << endl;
	for (int i = 0; i < m; i++) {
		cout << setprecision(8) << setw(12) << x[i] << endl;
}
	cout << endl;
}
//主函数
int main() {
	double a[4][4] = { {2,1,-1,1},{1,5,2,7},{-1,2,10,1},{1,7,1,11} };
	double b[4] = { 13,-9,6,0 };
	double L[4][4];//L 矩阵
	double LT[4][4];//LT 矩阵
	double* y;//中间矩阵 y
	InitL(L, 4, 4);//初始化矩阵 L
	Cholesky(a, L, 4, 4);//Cholesky 分解矩阵 A
	Display(L, LT, 4, 4);//显示矩阵 L 和 LT
	//求解原方程 Ax=b 中间步骤：L*y = b; LT*x = y;
	y = SolveOne(L, b, 4, 4);
	SolveTwo(L, y, 4, 4);
	return 0;
}

数值结果：

分析：
首先对 A 矩阵进行 Cholesky 分解，得到 L 矩阵，然后通过中间矩阵 y 求解矩阵 x。