diff options
Diffstat (limited to 'Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h')
| -rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 29 |
1 files changed, 15 insertions, 14 deletions
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 76343640d..2c53655d1 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -82,6 +82,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver /** \brief Scalar type for matrices of type \p _MatrixType. */ typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; /** \brief Real scalar type for \p _MatrixType. * @@ -105,7 +106,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * perform decompositions via compute(const MatrixType&, bool) or * compute(const MatrixType&, const MatrixType&, bool). This constructor * can only be used if \p _MatrixType is a fixed-size matrix; use - * SelfAdjointEigenSolver(int) for dynamic-size matrices. + * SelfAdjointEigenSolver(Index) for dynamic-size matrices. * * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out @@ -132,7 +133,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa compute(const MatrixType&, bool) for an example */ - SelfAdjointEigenSolver(int size) + SelfAdjointEigenSolver(Index size) : m_eivec(size, size), m_eivalues(size), m_tridiag(size), @@ -379,8 +380,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Implemented from Golub's "Matrix Computations", algorithm 8.3.2: * "implicit symmetric QR step with Wilkinson shift" */ -template<typename RealScalar, typename Scalar> -static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n); +template<typename RealScalar, typename Scalar, typename Index> +static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n); template<typename MatrixType> SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) @@ -389,7 +390,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute( m_eigenvectorsOk = computeEigenvectors; #endif assert(matrix.cols() == matrix.rows()); - int n = matrix.cols(); + Index n = matrix.cols(); m_eivalues.resize(n,1); m_eivec.resize(n,n); @@ -407,11 +408,11 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute( if (computeEigenvectors) m_eivec = m_tridiag.matrixQ(); - int end = n-1; - int start = 0; + Index end = n-1; + Index start = 0; while (end>0) { - for (int i = start; i<end; ++i) + for (Index i = start; i<end; ++i) if (ei_isMuchSmallerThan(ei_abs(m_subdiag[i]),(ei_abs(diag[i])+ei_abs(diag[i+1])))) m_subdiag[i] = 0; @@ -430,9 +431,9 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute( // Sort eigenvalues and corresponding vectors. // TODO make the sort optional ? // TODO use a better sort algorithm !! - for (int i = 0; i < n-1; ++i) + for (Index i = 0; i < n-1; ++i) { - int k; + Index k; m_eivalues.segment(i,n-i).minCoeff(&k); if (k > 0) { @@ -473,7 +474,7 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors { // transform back the eigen vectors: evecs = inv(U) * evecs cholB.matrixU().solveInPlace(m_eivec); - for (int i=0; i<m_eivec.cols(); ++i) + for (Index i=0; i<m_eivec.cols(); ++i) m_eivec.col(i) = m_eivec.col(i).normalized(); } return *this; @@ -482,8 +483,8 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors #endif // EIGEN_HIDE_HEAVY_CODE #ifndef EIGEN_EXTERN_INSTANTIATIONS -template<typename RealScalar, typename Scalar> -static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n) +template<typename RealScalar, typename Scalar, typename Index> +static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) { RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5); RealScalar e2 = ei_abs2(subdiag[end-1]); @@ -491,7 +492,7 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st RealScalar x = diag[start] - mu; RealScalar z = subdiag[start]; - for (int k = start; k < end; ++k) + for (Index k = start; k < end; ++k) { PlanarRotation<RealScalar> rot; rot.makeGivens(x, z); |