GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-23 05:14:26 +08:00
Code Issues Packages Projects Releases Wiki Activity

compute the rank on the fly rather than at the end, and stop early

Browse Source 
in the case of non-full rank (so big optimization in case the rank is low)
This commit is contained in:
Benoit Jacob 2009-01-26 16:14:20 +00:00
parent 9ea1050281
commit b05c5dd1be
1 changed files with 18 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
Eigen/src/LU/LU.h
View File

@ -86,14 +86,6 @@ template<typename MatrixType> class LU
MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns. MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns.
> ImageResultType; > ImageResultType;
typedef Matrix<typename MatrixType::Scalar,
MatrixType::RowsAtCompileTime,
MatrixType::RowsAtCompileTime,
MatrixType::Options,
MatrixType::MaxRowsAtCompileTime,
MatrixType::MaxRowsAtCompileTime
> MatrixLType;
/** Constructor. /** Constructor.
* *
* \param matrix the matrix of which to compute the LU decomposition. * \param matrix the matrix of which to compute the LU decomposition.
@ -343,16 +335,31 @@ LU<MatrixType>::LU(const MatrixType& matrix)
int number_of_transpositions = 0; int number_of_transpositions = 0;
RealScalar biggest = RealScalar(0); RealScalar biggest = RealScalar(0);
m_rank = size;
for(int k = 0; k < size; ++k) for(int k = 0; k < size; ++k)
{ {
int row_of_biggest_in_corner, col_of_biggest_in_corner; int row_of_biggest_in_corner, col_of_biggest_in_corner;
RealScalar biggest_in_corner; RealScalar biggest_in_corner;
biggest_in_corner = m_lu.corner(Eigen::BottomRight, rows-k, cols-k) biggest_in_corner = m_lu.corner(Eigen::BottomRight, rows-k, cols-k)
.cwise().abs() .cwise().abs()
.maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
row_of_biggest_in_corner += k; row_of_biggest_in_corner += k;
col_of_biggest_in_corner += k; col_of_biggest_in_corner += k;
if(k==0) biggest = biggest_in_corner;
// if the corner is negligible, then we have less than full rank, and we can finish early
if(ei_isMuchSmallerThan(biggest_in_corner, biggest))
{
m_rank = k;
for(int i = k; i < size; i++)
{
rows_transpositions.coeffRef(i) = i;
cols_transpositions.coeffRef(i) = i;
}
break;
}
rows_transpositions.coeffRef(k) = row_of_biggest_in_corner; rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
cols_transpositions.coeffRef(k) = col_of_biggest_in_corner; cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
if(k != row_of_biggest_in_corner) { if(k != row_of_biggest_in_corner) {
@ -363,12 +370,8 @@ LU<MatrixType>::LU(const MatrixType& matrix)
m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner)); m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
++number_of_transpositions; ++number_of_transpositions;
} }
if(k==0) biggest = biggest_in_corner;
const Scalar lu_k_k = m_lu.coeff(k,k);
if(ei_isMuchSmallerThan(lu_k_k, biggest)) continue;
if(k<rows-1) if(k<rows-1)
m_lu.col(k).end(rows-k-1) /= lu_k_k; m_lu.col(k).end(rows-k-1) /= m_lu.coeff(k,k);
if(k<size-1) if(k<size-1)
for(int col = k + 1; col < cols; ++col) for(int col = k + 1; col < cols; ++col)
m_lu.col(col).end(rows-k-1) -= m_lu.col(k).end(rows-k-1) * m_lu.coeff(k,col); m_lu.col(col).end(rows-k-1) -= m_lu.col(k).end(rows-k-1) * m_lu.coeff(k,col);
@ -383,12 +386,6 @@ LU<MatrixType>::LU(const MatrixType& matrix)
std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k))); std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_det_pq = (number_of_transpositions%2) ? -1 : 1;
RealScalar biggest_diagonal_coeff = m_lu.diagonal().cwise().abs().maxCoeff();
m_rank = 0;
for(int k = 0; k < size; ++k)
if(!ei_isMuchSmallerThan(m_lu.diagonal().coeff(k), biggest_diagonal_coeff))
++m_rank;
} }
template<typename MatrixType> template<typename MatrixType>