mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

- add kernel computation using the triangular solver

Browse Source 
- take advantage of the fact that our LU dec sorts the eigenvalues of U
in decreasing order
- add meta selector in determinant
This commit is contained in:
Benoit Jacob 2008-08-07 21:48:21 +00:00
parent 58ba9ca72f
commit 5f350448e1
2 changed files with 123 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
Eigen/src/LU/Determinant.h
View File

@ -67,6 +67,73 @@ const typename Derived::Scalar ei_bruteforce_det(const MatrixBase<Derived>& m)
}
}
const int TriangularDeterminant = 0;
template<typename Derived,
int DeterminantType =
(Derived::Flags & (UpperTriangularBit | LowerTriangularBit))
? TriangularDeterminant : Derived::RowsAtCompileTime
> struct ei_determinant_impl
{
static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
{
return m.lu().determinant();
}
};
template<typename Derived> struct ei_determinant_impl<Derived, TriangularDeterminant>
{
static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
{
if (Derived::Flags & UnitDiagBit)
return 1;
else if (Derived::Flags & ZeroDiagBit)
return 0;
else
return m.diagonal().redux(ei_scalar_product_op<typename ei_traits<Derived>::Scalar>());
}
};
template<typename Derived> struct ei_determinant_impl<Derived, 1>
{
static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
{
return m.coeff(0,0);
}
};
template<typename Derived> struct ei_determinant_impl<Derived, 2>
{
static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
{
return m.coeff(0,0) * m.coeff(1,1) - m.coeff(1,0) * m.coeff(0,1);
}
};
template<typename Derived> struct ei_determinant_impl<Derived, 3>
{
static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
{
return ei_bruteforce_det3_helper(m,0,1,2)
- ei_bruteforce_det3_helper(m,1,0,2)
+ ei_bruteforce_det3_helper(m,2,0,1);
}
};
template<typename Derived> struct ei_determinant_impl<Derived, 4>
{
static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
{
// trick by Martin Costabel to compute 4x4 det with only 30 muls
return ei_bruteforce_det4_helper(m,0,1,2,3)
- ei_bruteforce_det4_helper(m,0,2,1,3)
+ ei_bruteforce_det4_helper(m,0,3,1,2)
+ ei_bruteforce_det4_helper(m,1,2,0,3)
- ei_bruteforce_det4_helper(m,1,3,0,2)
+ ei_bruteforce_det4_helper(m,2,3,0,1);
}
};
/** \lu_module
*
* \returns the determinant of this matrix
@ -75,17 +142,7 @@ template<typename Derived>
typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
{
assert(rows() == cols());
if (Derived::Flags & (UpperTriangularBit | LowerTriangularBit))
{
if (Derived::Flags & UnitDiagBit)
return 1;
else if (Derived::Flags & ZeroDiagBit)
return 0;
else
return derived().diagonal().redux(ei_scalar_product_op<Scalar>());
}
else if(rows() <= 4) return ei_bruteforce_det(derived());
else return lu().determinant();
return ei_determinant_impl<Derived>::run(derived());
}
#endif // EIGEN_DETERMINANT_H

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

@ -51,8 +51,15 @@ template<typename MatrixType> class LU
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Matrix<int, MatrixType::ColsAtCompileTime, 1> IntRowVectorType;
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
enum { MaxKerDimAtCompileTime = EIGEN_ENUM_MIN(
MatrixType::MaxColsAtCompileTime,
MatrixType::MaxRowsAtCompileTime)
};
LU(const MatrixType& matrix);
@ -81,9 +88,9 @@ template<typename MatrixType> class LU
return m_q;
}
inline const Matrix<Scalar, MatrixType::RowsAtCompileTime, Dynamic,
MatrixType::MaxRowsAtCompileTime,
MatrixType::MaxColsAtCompileTime> kernel() const;
inline const Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, Dynamic,
MatrixType::MaxColsAtCompileTime,
LU<MatrixType>::MaxKerDimAtCompileTime> kernel() const;
template<typename OtherDerived>
typename ProductReturnType<Transpose<MatrixType>, OtherDerived>::Type::Eval
@ -131,7 +138,6 @@ template<typename MatrixType> class LU
IntColVectorType m_p;
IntRowVectorType m_q;
int m_det_pq;
Scalar m_biggest_eigenvalue_of_u;
int m_rank;
};
@ -173,8 +179,7 @@ LU<MatrixType>::LU(const MatrixType& matrix)
if(k==0) biggest = biggest_in_corner;
const Scalar lu_k_k = m_lu.coeff(k,k);
std::cout << lu_k_k << " " << biggest << std::endl;
if(ei_isMuchSmallerThan(lu_k_k, biggest)) { std::cout << "hello" << std::endl; continue; }
if(ei_isMuchSmallerThan(lu_k_k, biggest)) continue;
if(k<rows-1)
m_lu.col(k).end(rows-k-1) /= lu_k_k;
if(k<size-1)
@ -192,35 +197,61 @@ LU<MatrixType>::LU(const MatrixType& matrix)
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
int index_of_biggest_in_diagonal;
m_lu.diagonal().cwise().abs().maxCoeff(&index_of_biggest_in_diagonal);
m_biggest_eigenvalue_of_u = m_lu.diagonal().coeff(index_of_biggest_in_diagonal);
m_rank = 0;
for(int k = 0; k < size; k++)
m_rank += !ei_isMuchSmallerThan(m_lu.diagonal().coeff(k), m_biggest_eigenvalue_of_u);
m_rank += !ei_isMuchSmallerThan(m_lu.diagonal().coeff(k),
m_lu.diagonal().coeff(0));
}
template<typename MatrixType>
typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
{
if(!isInvertible()) return Scalar(0);
Scalar res = m_det_pq;
for(int k = 0; k < m_lu.diagonal().size(); k++) res *= m_lu.diagonal().coeff(k);
return res;
return m_lu.diagonal().redux(ei_scalar_product_op<Scalar>()) * Scalar(m_det_pq);
}
#if 0
template<typename MatrixType>
inline const Matrix<Scalar, RowsAtCompileTime, Dynamic,
MaxRowsAtCompileTime, MaxColsAtCompileTime> kernel() const
inline const Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, Dynamic,
MatrixType::MaxColsAtCompileTime,
LU<MatrixType>::MaxKerDimAtCompileTime>
LU<MatrixType>::kernel() const
{
Matrix<Scalar, RowsAtCompileTime, Dynamic,
MaxRowsAtCompileTime, MaxColsAtCompileTime>
result(m_lu.rows(), dimensionOfKernel());
ei_assert(!isInvertible());
const int dimker = dimensionOfKernel(), rows = m_lu.rows(), cols = m_lu.cols();
Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic,
MatrixType::MaxColsAtCompileTime,
LU<MatrixType>::MaxKerDimAtCompileTime>
result(cols, dimker);
/* Let us use the following lemma:
*
* Lemma: If the matrix A has the LU decomposition PAQ = LU,
* then Ker A = Q( Ker U ).
*
* Proof: trivial: just keep in mind that P, Q, L are invertible.
*/
/* Thus, all we need to do is to compute Ker U, and then apply Q.
*
* U is upper triangular, with eigenvalues sorted in decreasing order of
* absolute value. Thus, the diagonal of U ends with exactly
* m_dimKer zero's. Let us use that to construct m_dimKer linearly
* independent vectors in Ker U.
*/
Matrix<Scalar, Dynamic, Dynamic, MatrixType::MaxColsAtCompileTime, MaxKerDimAtCompileTime>
y(-m_lu.corner(TopRight, m_rank, dimker));
m_lu.corner(TopLeft, m_rank, m_rank)
.template marked<Upper>()
.inverseProductInPlace(y);
for(int i = 0; i < m_rank; i++)
result.row(m_q.coeff(i)) = y.row(i);
for(int i = m_rank; i < cols; i++) result.row(m_q.coeff(i)).setZero();
for(int k = 0; k < dimker; k++) result.coeffRef(m_q.coeff(m_rank+k), k) = Scalar(1);
return result;
}
#endif
/** \return the LU decomposition of \c *this.
*