LUP分解的C++实现

LUP分解的C++实现

简单粗暴上代码

在这里插入代码片

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
using namespace std;

void Swap(double* a, double* b) {//交换函数 
	double temp = *a;
	*a = *b;
	*b = temp;
}
void Swap1(int* a, int* b) {//交换函数 
	int temp = *a;
	*a = *b;
	*b = temp;
}


double Round(double dVal, short iPlaces) {//四舍五入保留小数位数 
	double dRetval;
	double dMod = 0.0000001;
	if (dVal < 0.0) dMod = -0.0000001;
	dRetval = dVal;
	dRetval += (5.0 / pow(10.0, iPlaces + 1.0));
	dRetval *= pow(10.0, iPlaces);
	dRetval = floor(dRetval + dMod);
	dRetval /= pow(10.0, iPlaces);
	return(dRetval);
}


void Input(double** A,  int n) {
	int i, j;
	cout << "请输入A矩阵:" << endl;
	for (i = 1; i <= n; i++) {
		for (j = 1; j <= n; j++) {
			cin >> A[i][j];
		}
	}

	/*cout << "请输入b数组:" << endl;
	for (i = 1; i <= n; i++) {
		cin >> b[i];
	}*/
}

void LUP_DECOMPOSITION(double** A, double** L, double** U, int* P, int n) {
	int i, j, k, k1;
	double temp, temp1;;
	for (i = 1; i <= n; i++) {
		P[i] = i;
	}

	for (k = 1; k <= n; k++) {
		int p = 0;
		for (i = k; i <= n; i++) {
			if (fabs(A[i][k]) > p) {
				p = fabs(A[i][k]);
				k1 = i;
			}
		}
		if (p == 0) {
			cout << "singular matrix" << endl;
		}
		else {
			Swap1(&P[k], &P[k1]);
		}

		for (i = 1; i <= n; i++) {
			Swap(&A[k][i], &A[k1][i]);
		}

		for (i = k + 1; i <= n; i++) {
			temp = A[i][k] / A[k][k];
			A[i][k] = Round(temp, 2);
			for (j = k + 1; j <= n; j++) {
				temp = A[i][k] * A[k][j];
				temp1 = Round(temp, 2);
				A[i][j] = A[i][j] - temp1;
			}
		}
	}
	for (i = 1; i <= n; i++) {
		for (j = 1; j <= n; j++) {
			if (i == j) {
				L[i][j] = 1;
			}
			else {
				if (i < j) {
					L[i][j] = 0;
				}
				else {
					L[i][j] = A[i][j];
				}
			}
		}
	}

	for (i = 1; i <= n; i++) {
		for (j = 1; j <= n; j++) {
			if (i <= j) {
				U[i][j] = A[i][j];
			}
			else {
				U[i][j] = 0;
			}
		}
	}

	for (i = 1; i <= n; i++) {
		P[i] -= 1;
	}

}

void Print(int* x, int n) {
	int i, j;
	for (i = 1; i <= n; i++) {
		cout << x[i] << "\t";
	}
}

void PrintL(double** L, int n) {
	int i, j;
	for (i = 1; i <= n; i++) {
		for (j = 1; j <= n; j++) {
			cout << L[i][j] << "\t";
		}
		cout << endl;
	}
}

int main() {
	int n;
	double* x = NULL;
	double* y = NULL;
	double** L = NULL;
	double** U = NULL;
	int* P = NULL;
	double** A = NULL;
	double* b = NULL;
	cin >> n;
	int i = 0;
	L = new double* [n + 1];
	for (i = 1; i <= n; i++) {
		L[i] = new double[n + 1];
	}

	U = new double* [n + 1];
	for (i = 1; i <= n; i++) {
		U[i] = new double[n + 1];
	}

	A = new double* [n + 1];
	for (i = 1; i <= n; i++) {
		A[i] = new double[n + 1];
	}

	P = new int[n ];
	x = new double[n + 1];
	y = new double[n + 1];
	b = new double[n + 1];

	Input(A,  n);
	LUP_DECOMPOSITION(A, L, U, P, n);
	//LUP_SOLVE(A, L, U, P, x, y, b, n);
	//Print(x, n);
	cout << "L:" << endl;
	PrintL(L, n);
	cout << "U:" << endl;
	PrintL(U, n);
	//PrintL(A, n);
	cout << "P:" << endl;
	Print(P, n);
	return 0;
}

/*

1 0 0 0 1 0 1 0
0 0 0 0 1 0 1 1
0 1 0 0 0 1 0 0
0 0 0 0 1 1 1 1
0 0 1 0 1 0 0 1
0 0 0 0 0 1 1 0
0 0 0 1 0 1 0 0
0 0 0 0 0 0 1 1

1 0 0 0 1 1 0 0
1 1 0 0 1 0 1 0
0 1 1 0 0 1 0 0
1 0 1 1 1 1 1 1
1 1 0 1 0 0 1 0
0 1 1 0 1 0 0 0
0 0 1 1 0 1 0 0
0 0 0 1 1 0 1 0
*/

如果对您有用请点赞哦，欢迎评论交流