矩阵的LU分解以及求逆矩阵--c语言
作者:互联网
matrix.h
#ifndef MATRIX_H_
#define MATRIX_H_
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <limits.h>
#include <stdbool.h>
#include <assert.h>
#ifdef N
#define PRINTMATRIX(A, S) \
for (int i = 0; i < N; i++) \
{ \
for (int j = 0; j < N; j++) \
{ \
printf("%lf\t", A[i][j]); \
} \
printf("\n"); \
} \
printf("=================================\n");
bool floatCompare(double f1, double f2)
{
return fabs(f1 - f2) < 1.0e-05;
}
//求对角矩阵的det
double DiagDet(double Matrix[N][N])
{
double res = Matrix[0][0];
for (int i = 0; i < N; i++)
{
res = res * Matrix[i][i];
}
return res;
}
void MatMul(double lhs[N][N], double rhs[N][N], double res[N][N])
{
for (int row = 0; row < N; row++)
{
for (int col = 0; col < N; col++)
{
double ret = 0.0;
for (int j = 0; j < N; j++)
{
ret += (lhs[row][j] * rhs[j][col]);
}
res[row][col] = ret;
}
}
}
double RAND(int start, int end)
{
return rand() % (end - start + 1) + start;
}
//矩阵的预处理
void MatrixPre(double A[N][N])
{
/*
前置的对矩阵进行处理
*/
for (int col = 0; col < N; col++)
{
int index, k;
k = 0;
double temp[N];
for (int j = 0; j < N; j++)
{
//找出目前的col的最大值的下标index 当前的index和对应的row的下标进行交换
temp[k] = A[j][col];
k++;
}
//找出temp数组中最大元的下标index 交换 A[col]和temp的结果
double max = temp[0];
index = 0;
for (int i = 0; i < N; i++)
{
if (temp[i] > max)
{
index = i;
max = temp[i];
}
}
//交换A[index] 和 A[col]
double swap[N];
for (int s = 0; s < N; s++)
{
swap[s] = A[index][s];
A[index][s] = A[col][s];
A[col][s] = swap[s];
}
}
}
void UpperInverse(double Matrix[N][N], double inverse[N][N])
{
for (int i = 0; i < N; i++)
{
inverse[i][i] = 1 / Matrix[i][i];
for (int k = i - 1; k >= 0; k--)
{
double s = 0;
for (int j = k + 1; j <= i; j++)
s = s + Matrix[k][j] * inverse[j][i];
inverse[k][i] = -s / Matrix[k][k];
}
}
}
void LowerInverse(double Matrix[N][N], double inverse[N][N])
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
if (i == j)
{
inverse[i][j] = 1 / Matrix[i][i];
}
else if (i < j)
{
inverse[i][j] = 0;
}
else if (i > j)
{
double s = 0.0;
for (int k = j; k <= i - 1; k++)
{
s = s + (Matrix[i][k] * inverse[k][j]);
}
inverse[i][j] = -s / (1 / Matrix[i][i]);
}
}
}
}
#endif
#endif
main.c
#define N 10
#include "matrix.h"
int main()
{
#ifdef N
/*
根据列最大元对对应得行交换
*/
printf("Notation: Matrix must be square Matrix!!!!\n");
printf("The size of Martix is %d\n", N);
double A[N][N];
//set Matrix
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
A[i][j] = RAND(10, 100);
double B[N][N];
double L[N][N];
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
if (i == j)
L[i][j] = 1;
else
L[i][j] = 0;
printf("=================================================\n");
printf("input Matrix:\n");
PRINTMATRIX(A, N);
MatrixPre(A);
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
B[i][j] = A[i][j];
printf("After swap rows:\n");
PRINTMATRIX(A, N);
for (int i = 0; i < N; i++)
{
for (int k = i + 1; k < N; k++)
{
double weight = (A[k][i] / (double)A[i][i]);
for (int s = 0; s < N; s++)
{
A[k][s] = A[k][s] - (A[i][s] * weight);
}
L[k][i] = weight;
}
}
printf("LU decomposition:\n");
printf("U Matrix as following:\n");
PRINTMATRIX(A, N);
printf("L matrix as following:\n");
PRINTMATRIX(L, N);
if (floatCompare(DiagDet(A), 0.0) || floatCompare(DiagDet(L), 0.0))
{
printf("There is no inverse matrix for this matrix!!!\n");
}
else
{
printf("Inverse Matrix as following:\n");
double L_inv[N][N];
double U_inv[N][N];
LowerInverse(L, L_inv);
UpperInverse(A, U_inv);
double inverse[N][N];
MatMul(U_inv, L_inv, inverse);
PRINTMATRIX(inverse, N);
printf("Check:\n");
double res[N][N];
MatMul(B, inverse, res);
PRINTMATRIX(res, N);
}
#else
printf("Please define Matrix size!!!\n");
#endif
system("pause");
return 0;
}
标签:求逆,Matrix,--,double,矩阵,++,int,printf,col 来源: https://www.cnblogs.com/dfxdd/p/14364424.html