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
*/
如果对您有用请点赞哦,欢迎评论交流