LU分解法求逆矩阵 C语言实现

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

最近在网上找了下，没有找到我想要的C语言版本，找到的也是错误的。故自己写了一个，并进行了相关测试，贴出来分享。

具体的LU分解算法就不细说了，随便找本书就知道了，关键是分解的处理流程，细节特别容易出错，一切都在代码里面。

#include <stdio.h>
#include <memory.h>
#include <stdlib.h>#define N 4
#define DEBUG 1//debug label，0即不打印相关结果，非0打印相关输出结果void matrix_inverse_LU(float a[][N])
{
float l[N][N], u[N][N];
float l_inverse[N][N], u_inverse[N][N];
float a_inverse[N][N];
int i, j, k;
float s, t;memset(l, 0, sizeof(l));
memset(u, 0, sizeof(u));
memset(l_inverse, 0, sizeof(l_inverse));
memset(u_inverse, 0, sizeof(u_inverse));
memset(a_inverse, 0, sizeof(u_inverse));for (i = 0; i < N;i++)//计算l矩阵对角线
{
l[i][i] = 1;
}for (i = 0;i < N;i++)
{
for (j = i;j < N;j++)
{
s = 0;
for (k = 0;k < i;k++)
{
s += l[i][k] * u[k][j];
}
u[i][j] = a[i][j] - s;//按行计算u值
}for (j = i + 1;j < N;j++)
{
s = 0;
for (k = 0; k < i; k++)
{
s += l[j][k] * u[k][i];
}
l[j][i] = (a[j][i] - s) / u[i][i];//按列计算l值
}
}for (i = 0;i < N;i++)//按行序，行内从高到低，计算l的逆矩阵
{
l_inverse[i][i] = 1;
}
for (i= 1;i < N;i++)
{
for (j = 0;j < i;j++)
{
s = 0;
for (k = 0;k < i;k++)
{
s += l[i][k] * l_inverse[k][j];
}
l_inverse[i][j] = -s;
}
}#if DEBUG
printf("test l_inverse:\n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
{
s = 0;
for (k = 0; k < N; k++)
{
s += l[i][k] * l_inverse[k][j];
}printf("%f ", s);
}
putchar('\n');
}
#endiffor (i = 0;i < N;i++)//按列序，列内按照从下到上，计算u的逆矩阵
{
u_inverse[i][i] = 1 / u[i][i];
}
for (i = 1;i < N;i++)
{
for (j = i - 1;j >=0;j--)
{
s = 0;
for (k = j + 1;k <= i;k++)
{
s += u[j][k] * u_inverse[k][i];
}
u_inverse[j][i] = -s / u[j][j];
}
}#if DEBUG
printf("test u_inverse:\n");
for (i = 0;i < N;i++)
{
for (j = 0;j < N;j++)
{
s = 0;
for (k = 0;k < N;k++)
{
s += u[i][k] * u_inverse[k][j];
}printf("%f ",s);
}
putchar('\n');
}
#endiffor (i = 0;i < N;i++)//计算矩阵a的逆矩阵
{
for (j = 0;j < N;j++)
{
for (k = 0;k < N;k++)
{
a_inverse[i][j] += u_inverse[i][k] * l_inverse[k][j];
}
}
}#if DEBUG
printf("test a:\n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
{
s = 0;
for (k = 0; k < N; k++)
{
s += a[i][k] * a_inverse[k][j];
}printf("%f ", s);
}
putchar('\n');
}
#endif
}void main()
{
int i, j, k;
float a[N][N];for (i = 0;i < N;i++)
{
for (j = 0;j < N;j++)
{
a[i][j] = rand() % 10;
}
}matrix_inverse_LU(a);
}