利用LU分解法求逆矩阵
利用LU分解法求逆矩阵
程序源代码：
void LU()
{
assert(ROWS == COLS);
double A[ROWS][COLS] ={
{-0.617, 53.9627, 2, 0},
{-53.9627, -0.617, 0, 2},
{49588, 0, -0.617, -53.9627},
{0, 49599, 53.9627,-0.617}};
//=====================================================================
//init the L and U
double L[ROWS][COLS], U[ROWS][COLS];
int i, j, k, m;
for (i = 0; i < ROWS; i++)
{
for (j = 0; j < COLS; j++)
{
U[i][j] = A[i][j];
L[i][j] = (i == j ? 1.0 : 0);
}
}
//=====================================================================
printf("the matrix is:\n");
for (i = 0; i < ROWS; i++)
{
for (j = 0; j < COLS; j++)
{
printf("%
}
printf("\n");
}
//=====================================================================
//Dootlittle: LU decompose
/************************************************************************/
/* method One: put L and U in two matrix, L and U,
第一种方法将L，U矩阵分别存放在两个矩阵中，占空间多些
*/
/************************************************************************/
for (j = 0; j < COLS - 1; j++)
{
for(i = j + 1; i < ROWS; i++)
{
if (U[j][j] == 0)
{
int n = i;
while (U[n][j] == 0 && n < ROWS)
n++;
for (k = j; k < COLS; k++)
{
double t = U[j][k];
U[j][k] = U[n][k];
U[n][k] = t;
}
}
L[i][j] = U[i][j] / U[j][j];
for(m = j; m < COLS; m++)
U[i][m] -= U[j][m] * L[i][j];
}
}
//=====================================================================
/************************************************************************/
/* method Two: put L and U in one matrix,
the Upper trianqular matrix is U,
the lower trianqular matrix is L.
第二种方法：将L，U矩阵放在同一个矩阵中
*/
/************************************************************************/
/* double temp = 0;
for (j = 0; j < COLS - 1; j++)
{
for(i = j + 1; i < ROWS; i++)
{
//exchange two rows if U[j][j] = 0;
if (U[j][j] == 0) //对第j列处理的时候，要判断是否为零，为零需要交换两行
{
int n = i;
while (U[n][j] == 0 && n < ROWS)
n++;
for (k = j; k < COLS; k++)
{
double t = U[j][k];
U[j][k] = U[n][k];
U[n][k] = t;
}
}
temp = U[i][j] / U[j][j];
for(m = j; m < COLS; m++)
U[i][m] -= U[j][m] * temp;
U[i][j] = temp;
}
}
*/
//=====================================================================
printf("============matrix L=============\n");
for (i = 0; i < ROWS; i++)
{
for (j = 0; j < COLS; j++)
{
printf("%
}
printf("\n");
}
printf("============matrix U=============\n");
for (i = 0; i < ROWS; i++)
{
for (j = 0; j < COLS; j++)
{
printf("%
}
printf("\n");
}
CvMat mat1 = cvMat(4, 4, CV_64FC1, L);
CvMat mat2 = cvMat(4, 4, CV_64FC1, U);
CvMat *mat3 = cvCreateMat(4, 4, CV_64FC1);
cvGEMM(&mat1, &mat2, 1, NULL, 0, mat3);
printf("the matrix is(L x U):\n");
for (i = 0; i < mat3->rows; i++)
{
for (j = 0; j < mat3->cols; j++)
{
printf("%
}
printf("\n");
}
}
方法一：
运行结果：
the matrix is:
-0.617 53.9627 2 0
-53.9627 -0.617 0 2
49588 0 -0.617 -53.9627
0 49599 53.9627 -0.617
============matrix L=============
1 0 0 0
87.4598 1 0 0
-80369.5 -918.811 1 0
0 -10.5079 -87.4798 1
============matrix U=============
-0.617 53.9627 2 0
0 -4720.18 -174.92 2
0 0 20.394 1783.66
0 0 0 156055
the matrix is(L x U):
-0.617 53.9627 2 0
-53.9627 -0.617 0 2
49588 0 -0.617 -53.9627
0 49599 53.9627 -0.617
方法二：
the matrix is:
-0.617 53.9627 2 0
-53.9627 -0.617 0 2
49588 0 -0.617 -53.9627
0 49599 53.9627 -0.617
============matrix U=============
-0.617 53.9627 2 0
87.4598 -4720.18 -174.92 2
-80369.5 -918.811 20.394 1783.66
0 -10.5079 -87.4798 156055
