矩阵分解之LU

分享我这几天整理的代码
原理参考如下：https://blog.csdn.net/u010945683/article/details/45803141

直接上代码，注释部分是求解线性方程的

void LU_Decom(double *matrix, double *L,double *U, double *b,int order)
{
	for (int i = 0; i < order; ++i)
	{
		for (int j = 0; j < i; ++j)
		{
			*(U + i*order + j) = 0;
		}//初始化U的下三角为0

		*(L + i*order + i) = 1;//初始化L的对角线为1

		for (int j = i + 1; j < order; ++j)
		{
			*(L + i*order + j) = 0;
		}//初始化L的上三角为0
	}
	for (int i = 0; i < order; ++i)
	{
		*(U + i) = *(matrix + i);
	}//计算U的第一行

	for (int i = 1; i < order; i++)
	{
		*(L + i*order) = *(matrix + i * order) / *U;
		//*(L + i * order) = *(matrix + i * order);
	}/
	float temp;
	for (int i = 1; i < order; ++i)
	{
		for (int j = i; j < order; ++j)
		{   
			temp = 0;
			for (int k = 0; k < i; ++k)
			{
				//cout << *(U + k*order + j) << ' ' << *(L + i*order + k) << endl;
				temp += (*(U + k * order + j) * (*(L + i * order + k)));;
			}

			*(U + i*order + j) = *(matrix + i*order + j) - temp;
		}//计算U的其他行
		for (int j = i + 1; j < order; ++j)
		{
			temp = 0;
			for (int k = 0; k < i; ++k)
			{
				temp += *(U + k*order + i)*(*(L + j*order + k));
			}
			*(L + j*order + i) = (*(matrix + j*order + i) - temp) / (*(U + i*order + i));
		}//计算L的其他行
	}
	//求解线性方程LY=B,UX=Y
	/*double *x = new double[order]();
	double *y = new double[order]();
	double sum;
	*y = *b;
	for (int i = 1; i < order; ++i)
	{
		sum = 0;
		for (int j = 0; j < i; ++j)
		{
			sum += *(L + i*order + j)*(*(y + j));
		}
		*(y + i) = *(b + i) - sum;
	}
	
	*(x + order - 1) = *(y + order - 1) / *(U+order*order-1);
	for (int i = order - 2; i>=0; i--)
	{
		sum = 0;
		for (int j = i + 1; j < order; ++j)
		{
			sum += *(U + i*order + j)*(*(x + j));
		}
		*(x + i) = (*(y + i) - sum)/ (*(U+(i+order)*i));
	}

	for (int i = 0; i < order; ++i)
	{
		cout << *(y + i) << endl;
	}

	for (int i = 0; i < order; ++i)
	{
		cout << *(x + i) << endl;
	}*/

	//delete [] y;
	//delete [] x;
}