diff options
Diffstat (limited to 'Eigen/src/Eigenvalues/Tridiagonalization.h')
-rw-r--r-- | Eigen/src/Eigenvalues/Tridiagonalization.h | 64 |
1 files changed, 37 insertions, 27 deletions
diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index 23ae748d4..fac79cf96 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -26,6 +26,11 @@ #ifndef EIGEN_TRIDIAGONALIZATION_H #define EIGEN_TRIDIAGONALIZATION_H +namespace internal { +template<typename MatrixType, typename CoeffVectorType> +void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs); +} + /** \eigenvalues_module \ingroup Eigenvalues_Module * * @@ -78,15 +83,15 @@ template<typename _MatrixType> class Tridiagonalization }; typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType; - typedef typename ei_plain_col_type<MatrixType, RealScalar>::type DiagonalType; + typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType; typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType; - typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex, + typedef typename internal::meta_if<NumTraits<Scalar>::IsComplex, typename Diagonal<MatrixType,0>::RealReturnType, Diagonal<MatrixType,0> >::ret DiagonalReturnType; - typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex, + typedef typename internal::meta_if<NumTraits<Scalar>::IsComplex, typename Diagonal< Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 >::RealReturnType, Diagonal< @@ -129,7 +134,7 @@ template<typename _MatrixType> class Tridiagonalization m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1), m_isInitialized(false) { - ei_tridiagonalization_inplace(m_matrix, m_hCoeffs); + internal::tridiagonalization_inplace(m_matrix, m_hCoeffs); m_isInitialized = true; } @@ -154,7 +159,7 @@ template<typename _MatrixType> class Tridiagonalization { m_matrix = matrix; m_hCoeffs.resize(matrix.rows()-1, 1); - ei_tridiagonalization_inplace(m_matrix, m_hCoeffs); + internal::tridiagonalization_inplace(m_matrix, m_hCoeffs); m_isInitialized = true; return *this; } @@ -177,7 +182,7 @@ template<typename _MatrixType> class Tridiagonalization */ inline CoeffVectorType householderCoefficients() const { - ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); + eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); return m_hCoeffs; } @@ -214,7 +219,7 @@ template<typename _MatrixType> class Tridiagonalization */ inline const MatrixType& packedMatrix() const { - ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); + eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); return m_matrix; } @@ -235,7 +240,7 @@ template<typename _MatrixType> class Tridiagonalization */ HouseholderSequenceType matrixQ() const { - ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); + eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1); } @@ -296,7 +301,7 @@ template<typename MatrixType> const typename Tridiagonalization<MatrixType>::DiagonalReturnType Tridiagonalization<MatrixType>::diagonal() const { - ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); + eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); return m_matrix.diagonal(); } @@ -304,7 +309,7 @@ template<typename MatrixType> const typename Tridiagonalization<MatrixType>::SubDiagonalReturnType Tridiagonalization<MatrixType>::subDiagonal() const { - ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); + eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); Index n = m_matrix.rows(); return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal(); } @@ -315,7 +320,7 @@ Tridiagonalization<MatrixType>::matrixT() const { // FIXME should this function (and other similar ones) rather take a matrix as argument // and fill it ? (to avoid temporaries) - ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); + eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); Index n = m_matrix.rows(); MatrixType matT = m_matrix; matT.topRightCorner(n-1, n-1).diagonal() = subDiagonal().template cast<Scalar>().conjugate(); @@ -327,6 +332,8 @@ Tridiagonalization<MatrixType>::matrixT() const return matT; } +namespace internal { + /** \internal * Performs a tridiagonal decomposition of the selfadjoint matrix \a matA in-place. * @@ -351,10 +358,10 @@ Tridiagonalization<MatrixType>::matrixT() const * \sa Tridiagonalization::packedMatrix() */ template<typename MatrixType, typename CoeffVectorType> -void ei_tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs) +void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs) { - ei_assert(matA.rows()==matA.cols()); - ei_assert(matA.rows()==hCoeffs.size()+1); + eigen_assert(matA.rows()==matA.cols()); + eigen_assert(matA.rows()==hCoeffs.size()+1); typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; @@ -371,9 +378,9 @@ void ei_tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs) matA.col(i).coeffRef(i+1) = 1; hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>() - * (ei_conj(h) * matA.col(i).tail(remainingSize))); + * (conj(h) * matA.col(i).tail(remainingSize))); - hCoeffs.tail(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1); + hCoeffs.tail(n-i-1) += (conj(h)*Scalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1); matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>() .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1); @@ -387,7 +394,7 @@ void ei_tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs) template<typename MatrixType, int Size=MatrixType::ColsAtCompileTime, bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex> -struct ei_tridiagonalization_inplace_selector; +struct tridiagonalization_inplace_selector; /** \brief Performs a full tridiagonalization in place * @@ -430,19 +437,19 @@ struct ei_tridiagonalization_inplace_selector; * \sa class Tridiagonalization */ template<typename MatrixType, typename DiagonalType, typename SubDiagonalType> -void ei_tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) +void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) { typedef typename MatrixType::Index Index; //Index n = mat.rows(); - ei_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1); - ei_tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ); + eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1); + tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ); } /** \internal * General full tridiagonalization */ template<typename MatrixType, int Size, bool IsComplex> -struct ei_tridiagonalization_inplace_selector +struct tridiagonalization_inplace_selector { typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType; typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType; @@ -451,7 +458,7 @@ struct ei_tridiagonalization_inplace_selector static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) { CoeffVectorType hCoeffs(mat.cols()-1); - ei_tridiagonalization_inplace(mat,hCoeffs); + tridiagonalization_inplace(mat,hCoeffs); diag = mat.diagonal().real(); subdiag = mat.template diagonal<-1>().real(); if(extractQ) @@ -464,7 +471,7 @@ struct ei_tridiagonalization_inplace_selector * Especially useful for plane fitting. */ template<typename MatrixType> -struct ei_tridiagonalization_inplace_selector<MatrixType,3,false> +struct tridiagonalization_inplace_selector<MatrixType,3,false> { typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; @@ -473,7 +480,7 @@ struct ei_tridiagonalization_inplace_selector<MatrixType,3,false> static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) { diag[0] = mat(0,0); - RealScalar v1norm2 = ei_abs2(mat(2,0)); + RealScalar v1norm2 = abs2(mat(2,0)); if(v1norm2 == RealScalar(0)) { diag[1] = mat(1,1); @@ -485,7 +492,7 @@ struct ei_tridiagonalization_inplace_selector<MatrixType,3,false> } else { - RealScalar beta = ei_sqrt(ei_abs2(mat(1,0)) + v1norm2); + RealScalar beta = sqrt(abs2(mat(1,0)) + v1norm2); RealScalar invBeta = RealScalar(1)/beta; Scalar m01 = mat(1,0) * invBeta; Scalar m02 = mat(2,0) * invBeta; @@ -508,16 +515,19 @@ struct ei_tridiagonalization_inplace_selector<MatrixType,3,false> * Trivial specialization for 1x1 matrices */ template<typename MatrixType, bool IsComplex> -struct ei_tridiagonalization_inplace_selector<MatrixType,1,IsComplex> +struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex> { typedef typename MatrixType::Scalar Scalar; template<typename DiagonalType, typename SubDiagonalType> static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ) { - diag(0,0) = ei_real(mat(0,0)); + diag(0,0) = real(mat(0,0)); if(extractQ) mat(0,0) = Scalar(1); } }; + +} // end namespace internal + #endif // EIGEN_TRIDIAGONALIZATION_H |