矩阵分解（LU分解）C语言实现

最新推荐文章于 2022-05-14 21:48:04 发布

黑曜

最新推荐文章于 2022-05-14 21:48:04 发布

阅读量3.9k

点赞数 1

文章标签： lu分解 c语言

本文链接：https://blog.csdn.net/weixin_41466979/article/details/103732172

版权

矩阵分解（LU分解）C语言实现

基于C代码实现的矩阵LU分解。网上找的其他代码都存在一些问题，例如对角元素存在0的处理。
我在某个C++版本上修改而成，经测试没问题，需要用的自取。

#define N 4//矩阵大小
#define DEBUG 1         //debug label，0即不打印相关结果，非0打印相关输出结果
 
 int main (void)
    {
    double a[4][4] = {
    {0.0, 2.0, 0.0, 1.0},
    {2.0, 2.0, 3.0, 2.0},
    {4.0, -3.0, 0.0, 1.0},
    {6.0, 1.0, -6.0, -5.0}    
    };
    double inva[4][4];

	matrix_inverse_LU（a,inva）;
        
    }

void matrix_inverse_LU(double (*a)[N], double (*y)[N])
{
		double b[N];
		int indx[N];
		printf("Origin matrix:");
	  output(a);
	  lucmp(a,indx);
	  int i,j;
	  for (i = 0; i < N; i++)
	  {
		for (j = 0; j < N; j++) 
		b[j] = 0;
		b[i] = 1.0;
		lubksb(a,b,indx);
		for (j = 0; j < N; j++) 
					y[j][i] = b[j];
	       }    
		printf("Its inverse matrix:\n");
		output(y);    
}

void lucmp(double (*a)[N],int *indx) {
		int imax = 0;
		double vv[N];
		int i,j,k;
		for (i = 0; i < N; i++) {
			double big = 0.0;
			for (j = 0; j < N; j++) {
				double temp = abs(a[i][j]);
				if (temp > big) big = temp;
			}
			vv[i] = 1.0 / big;
		}
		for (j = 0; j < N; j++) {
			for (i = 0; i < j; i++) {
				double sum = a[i][j];
				for (k = 0; k < i; k++) sum -= a[i][k] * a[k][j];
				a[i][j] = sum;
			}
			double big = 0.0;
			for (i = j; i < N; i++) {
				double sum = a[i][j];
				for (k = 0; k < j; k++) sum -= a[i][k] * a[k][j];
				a[i][j] = sum;
				double dum = vv[i] * abs(sum);
				if (dum >= big) {
					big = dum;
					imax = i;
				}
			}
			if (j != imax) {
				for (i = 0; i < N; i++) {
					double mid = a[imax][i];
					a[imax][i] = a[j][i];
					a[j][i] = mid;
				}
				double mid = vv[j];
				vv[j] = vv[imax];
				vv[imax] = mid;
			}
			indx[j] = imax;
			if (j != N - 1) {
				double dum = 1.0/a[j][j];
				for (i = j + 1; i < N; i++) a[i][j] *= dum;
			}
		}
	}
 
void lubksb(double (*a)[N],double *b,int *indx) {
		int ii = 0;
		int i,j;
		for (i = 0; i < N; i++) {
			int ip = indx[i];
			double sum = b[ip];
			b[ip] = b[i];
			if (ii != 0)
			 for (j = ii - 1; j < i; j++) 
				 sum -= a[i][j] * b[j]; 
			else
			 ii = i + 1;
			b[i] = sum;
		}
		// x
		for (i = N - 1; i >= 0; i--) {
			double sum = b[i];
			for (j = i + 1; j < N; j++) 
				sum -= a[i][j]*b[j];
			b[i] = sum / a[i][i];
		}
	}


void output(double (*a)[N]) {
#ifdef DEBUG
	int i;
	for (i = 0; i < N; i++) 
	{
		printf(" [%d]| ",i+1);
		int j;
		for (j=0;j<N;j++)
			printf("%.4f",a[i][j]);
		printf(" |\n ");
	}
	printf("-----------------------------------------------\n");
#endif
	}

C++版本源码：
https://blog.csdn.net/BoyMgl/article/details/1914288