矩阵的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;
}