# LU 分解 （LU Decomposition）

LU分解是指将一个 NxN 矩阵 A 分解为一个上三角矩阵 U 和下三角矩阵 L 的过程， 即： ＬUＡ。


/*
*　Decompose matrix: PA=LU
*
*/
template< int M>
bool LUDecomposition::decompose(const Matrix<M, M> &matrix,
Matrix<M, M> &ll,
Matrix<M, M> &uu,
int pivot[M])
{
uu = matrix;
ll.identity();

for (int iRow=0; iRow<M; ++iRow)
{
// find pivot rows
pivot[iRow] = iRow;
float maxValue = fabs(uu(iRow, iRow));
for (int iNextRow=iRow+1; iNextRow<M; ++iNextRow)
{
if (maxValue + FLOAT_EPS < fabs(uu(iNextRow, iRow)) )
{
pivot[iRow] = iNextRow;
maxValue = fabs(uu(iNextRow, iRow));
}
}

// if the matrix is singular, return error
if (almostZero(maxValue))
return false;

//  if the pivot row differs from the current row, then
//  interchange the two rows.
if (pivot[iRow] != iRow)
{
uu.swapRows(iRow, pivot[iRow]);
}

// update lower matrix and upper matrix
for (int iNextRow = iRow+1; iNextRow<M; ++iNextRow)
{
double factor = uu(iNextRow, iRow) / uu(iRow, iRow);
// lower matrix
ll(iNextRow, iRow) = (float)factor;

// upper matrix
uu(iNextRow, iRow) = 0.f;
for (int iCol= iRow+1; iCol<M;++iCol)
{
uu(iNextRow, iCol) -= uu(iRow, iCol) * factor;
}
}
}
return true;
}

/*
*　Solve LUx=b
*
*/
template<int M>
bool LUDecomposition::solve(const Matrix<M, M> &ll,
const Matrix<M, M> &uu,
float bb[M],
float xx[M])
{
// first, let y=Ux, solve Ly=b
float yy[M];
for (int ii=0; ii< M; ++ii)
{
yy[ii] = bb[ii];
for(int jj=0; jj<ii; ++jj)
yy[ii] -= yy[jj] * ll(ii, jj);
}

// then solve Ux = y
for (int ii=M-1; ii>=0; --ii)
{
if (almostZero(uu(ii, ii)))
return false;

xx[ii] = yy[ii];

for (int jj=ii+1; jj<M; ++jj)
{
xx[ii] -= xx[jj] * uu(ii, jj);
}

xx[ii] /=  uu(ii, ii);
}
return true;
}
posted @ 2010-06-19 00:04  拔剑四顾  阅读(4451)  评论(1编辑  收藏  举报