devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

merge with default branch

Browse Source
This commit is contained in:
Gael Guennebaud
2013-11-05 10:31:59 +01:00
parent 7c9cdd6030 ddf7753631
commit 1bb1a57ef7
66 changed files with 1152 additions and 290 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/Eigenvalues/EigenSolver.h
View File

@@ -568,7 +568,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
}
// Overflow control
using std::max;
EIGEN_USING_STD_MATH(max);
Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
if ((eps * t) * t > Scalar(1))
m_matT.block(i, n-1, size-i, 2) /= t;

40
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@@ -112,6 +112,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
* Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp
* Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out
*/
EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver()
: m_eivec(),
m_eivalues(),
@@ -131,6 +132,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*
* \sa compute() for an example
*/
EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver(Index size)
: m_eivec(size, size),
m_eivalues(size),
@@ -153,6 +155,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*
* \sa compute(const MatrixType&, int)
*/
EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors)
: m_eivec(matrix.rows(), matrix.cols()),
m_eivalues(matrix.cols()),
@@ -192,6 +195,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*
* \sa SelfAdjointEigenSolver(const MatrixType&, int)
*/
EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
/** \brief Computes eigendecomposition of given matrix using a direct algorithm
@@ -208,6 +212,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*
* \sa compute(const MatrixType&, int options)
*/
EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors);
/**
@@ -241,6 +246,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*
* \sa eigenvalues()
*/
EIGEN_DEVICE_FUNC
const MatrixType& eigenvectors() const
{
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
@@ -263,6 +269,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*
* \sa eigenvectors(), MatrixBase::eigenvalues()
*/
EIGEN_DEVICE_FUNC
const RealVectorType& eigenvalues() const
{
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
@@ -287,6 +294,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
* \sa operatorInverseSqrt(),
* \ref MatrixFunctions_Module "MatrixFunctions Module"
*/
EIGEN_DEVICE_FUNC
MatrixType operatorSqrt() const
{
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
@@ -312,6 +320,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
* \sa operatorSqrt(), MatrixBase::inverse(),
* \ref MatrixFunctions_Module "MatrixFunctions Module"
*/
EIGEN_DEVICE_FUNC
MatrixType operatorInverseSqrt() const
{
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
@@ -323,6 +332,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*
* \returns \c Success if computation was succesful, \c NoConvergence otherwise.
*/
EIGEN_DEVICE_FUNC
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
@@ -337,6 +347,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
static const int m_maxIterations = 30;
#ifdef EIGEN2_SUPPORT
EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors)
: m_eivec(matrix.rows(), matrix.cols()),
m_eivalues(matrix.cols()),
@@ -346,6 +357,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
compute(matrix, computeEigenvectors);
}
EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
: m_eivec(matA.cols(), matA.cols()),
m_eivalues(matA.cols()),
@@ -355,11 +367,13 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
}
EIGEN_DEVICE_FUNC
void compute(const MatrixType& matrix, bool computeEigenvectors)
{
compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
}
EIGEN_DEVICE_FUNC
void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
{
compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
@@ -393,10 +407,12 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*/
namespace internal {
template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
EIGEN_DEVICE_FUNC
static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
}
template<typename MatrixType>
EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
::compute(const MatrixType& matrix, int options)
{
@@ -538,6 +554,7 @@ ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag
template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues
{
EIGEN_DEVICE_FUNC
static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options)
{ eig.compute(A,options); }
};
@@ -548,12 +565,13 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
typedef typename SolverType::RealVectorType VectorType;
typedef typename SolverType::Scalar Scalar;
EIGEN_DEVICE_FUNC
static inline void computeRoots(const MatrixType& m, VectorType& roots)
{
using std::sqrt;
using std::atan2;
using std::cos;
using std::sin;
EIGEN_USING_STD_MATH(sqrt)
EIGEN_USING_STD_MATH(atan2)
EIGEN_USING_STD_MATH(cos)
EIGEN_USING_STD_MATH(sin)
const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
const Scalar s_sqrt3 = sqrt(Scalar(3.0));
@@ -588,15 +606,16 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
// Sort in increasing order.
if (roots(0) >= roots(1))
std::swap(roots(0),roots(1));
internal::swap(roots(0),roots(1));
if (roots(1) >= roots(2))
{
std::swap(roots(1),roots(2));
internal::swap(roots(1),roots(2));
if (roots(0) >= roots(1))
std::swap(roots(0),roots(1));
internal::swap(roots(0),roots(1));
}
}
EIGEN_DEVICE_FUNC
static inline void run(SolverType& solver, const MatrixType& mat, int options)
{
using std::sqrt;
@@ -717,12 +736,14 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
};
// 2x2 direct eigenvalues decomposition, code from Hauke Heibel
template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2,false>
template<typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,2,false>
{
typedef typename SolverType::MatrixType MatrixType;
typedef typename SolverType::RealVectorType VectorType;
typedef typename SolverType::Scalar Scalar;
EIGEN_DEVICE_FUNC
static inline void computeRoots(const MatrixType& m, VectorType& roots)
{
using std::sqrt;
@@ -732,6 +753,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2
roots(1) = t1 + t0;
}
EIGEN_DEVICE_FUNC
static inline void run(SolverType& solver, const MatrixType& mat, int options)
{
using std::sqrt;
@@ -785,6 +807,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2
}
template<typename MatrixType>
EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
::computeDirect(const MatrixType& matrix, int options)
{
@@ -794,6 +817,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
namespace internal {
template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
EIGEN_DEVICE_FUNC
static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
{
using std::abs;