diff options
Diffstat (limited to 'eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h')
| -rw-r--r-- | eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 273 |
1 files changed, 171 insertions, 102 deletions
diff --git a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 1131c8a..9ddd553 100644 --- a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -20,6 +20,8 @@ class GeneralizedSelfAdjointEigenSolver; namespace internal { template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues; +template<typename MatrixType, typename DiagType, typename SubDiagType> +ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec); } /** \eigenvalues_module \ingroup Eigenvalues_Module @@ -79,7 +81,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; + typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 typedef Matrix<Scalar,Size,Size,ColMajor,MaxColsAtCompileTime,MaxColsAtCompileTime> EigenvectorsType; @@ -100,6 +102,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver */ typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; typedef Tridiagonalization<MatrixType> TridiagonalizationType; + typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType; /** \brief Default constructor for fixed-size matrices. * @@ -111,6 +114,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver() : m_eivec(), m_eivalues(), @@ -130,7 +134,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa compute() for an example */ - SelfAdjointEigenSolver(Index size) + EIGEN_DEVICE_FUNC + explicit SelfAdjointEigenSolver(Index size) : m_eivec(size, size), m_eivalues(size), m_subdiag(size > 1 ? size - 1 : 1), @@ -152,13 +157,15 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa compute(const MatrixType&, int) */ - SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors) + template<typename InputType> + EIGEN_DEVICE_FUNC + explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors) : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()), m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), m_isInitialized(false) { - compute(matrix, options); + compute(matrix.derived(), options); } /** \brief Computes eigendecomposition of given matrix. @@ -191,24 +198,45 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa SelfAdjointEigenSolver(const MatrixType&, int) */ - SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors); + template<typename InputType> + EIGEN_DEVICE_FUNC + SelfAdjointEigenSolver& compute(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors); - /** \brief Computes eigendecomposition of given matrix using a direct algorithm + /** \brief Computes eigendecomposition of given matrix using a closed-form algorithm * * This is a variant of compute(const MatrixType&, int options) which * directly solves the underlying polynomial equation. * - * Currently only 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d). + * Currently only 2x2 and 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d). * - * This method is usually significantly faster than the QR algorithm + * This method is usually significantly faster than the QR iterative algorithm * but it might also be less accurate. It is also worth noting that * for 3x3 matrices it involves trigonometric operations which are * not necessarily available for all scalar types. + * + * For the 3x3 case, we observed the following worst case relative error regarding the eigenvalues: + * - double: 1e-8 + * - float: 1e-3 * * \sa compute(const MatrixType&, int options) */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors); + /** + *\brief Computes the eigen decomposition from a tridiagonal symmetric matrix + * + * \param[in] diag The vector containing the diagonal of the matrix. + * \param[in] subdiag The subdiagonal of the matrix. + * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. + * \returns Reference to \c *this + * + * This function assumes that the matrix has been reduced to tridiagonal form. + * + * \sa compute(const MatrixType&, int) for more information + */ + SelfAdjointEigenSolver& computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options=ComputeEigenvectors); + /** \brief Returns the eigenvectors of given matrix. * * \returns A const reference to the matrix whose columns are the eigenvectors. @@ -227,6 +255,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa eigenvalues() */ + EIGEN_DEVICE_FUNC const EigenvectorsType& eigenvectors() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -249,6 +278,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."); @@ -270,9 +300,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out * - * \sa operatorInverseSqrt(), - * \ref MatrixFunctions_Module "MatrixFunctions Module" + * \sa operatorInverseSqrt(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a> */ + EIGEN_DEVICE_FUNC MatrixType operatorSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -295,9 +325,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out * - * \sa operatorSqrt(), MatrixBase::inverse(), - * \ref MatrixFunctions_Module "MatrixFunctions Module" + * \sa operatorSqrt(), MatrixBase::inverse(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a> */ + EIGEN_DEVICE_FUNC MatrixType operatorInverseSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -309,6 +339,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."); @@ -322,36 +353,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver */ static const int m_maxIterations = 30; - #ifdef EIGEN2_SUPPORT - SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors) - : m_eivec(matrix.rows(), matrix.cols()), - m_eivalues(matrix.cols()), - m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), - m_isInitialized(false) - { - compute(matrix, computeEigenvectors); - } - - SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) - : m_eivec(matA.cols(), matA.cols()), - m_eivalues(matA.cols()), - m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1), - m_isInitialized(false) - { - static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); - } - - void compute(const MatrixType& matrix, bool computeEigenvectors) - { - compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); - } - - void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) - { - compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); - } - #endif // EIGEN2_SUPPORT - protected: static void check_template_parameters() { @@ -366,6 +367,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver bool m_eigenvectorsOk; }; +namespace internal { /** \internal * * \eigenvalues_module \ingroup Eigenvalues_Module @@ -373,8 +375,12 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Performs a QR step on a tridiagonal symmetric matrix represented as a * pair of two vectors \a diag and \a subdiag. * - * \param matA the input selfadjoint matrix - * \param hCoeffs returned Householder coefficients + * \param diag the diagonal part of the input selfadjoint tridiagonal matrix + * \param subdiag the sub-diagonal part of the input selfadjoint tridiagonal matrix + * \param start starting index of the submatrix to work on + * \param end last+1 index of the submatrix to work on + * \param matrixQ pointer to the column-major matrix holding the eigenvectors, can be 0 + * \param n size of the input matrix * * For compilation efficiency reasons, this procedure does not use eigen expression * for its arguments. @@ -382,17 +388,21 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Implemented from Golub's "Matrix Computations", algorithm 8.3.2: * "implicit symmetric QR step with Wilkinson shift" */ -namespace internal { -template<typename RealScalar, typename Scalar, typename Index> +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> +template<typename InputType> +EIGEN_DEVICE_FUNC SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> -::compute(const MatrixType& matrix, int options) +::compute(const EigenBase<InputType>& a_matrix, int options) { check_template_parameters(); + const InputType &matrix(a_matrix.derived()); + using std::abs; eigen_assert(matrix.cols() == matrix.rows()); eigen_assert((options&~(EigVecMask|GenEigMask))==0 @@ -404,7 +414,8 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> if(n==1) { - m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0)); + m_eivec = matrix; + m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0)); if(computeEigenvectors) m_eivec.setOnes(n,n); m_info = Success; @@ -424,19 +435,74 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> mat.template triangularView<Lower>() /= scale; m_subdiag.resize(n-1); internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors); + + m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); + // scale back the eigen values + m_eivalues *= scale; + + m_isInitialized = true; + m_eigenvectorsOk = computeEigenvectors; + return *this; +} + +template<typename MatrixType> +SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> +::computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options) +{ + //TODO : Add an option to scale the values beforehand + bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; + + m_eivalues = diag; + m_subdiag = subdiag; + if (computeEigenvectors) + { + m_eivec.setIdentity(diag.size(), diag.size()); + } + m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); + + m_isInitialized = true; + m_eigenvectorsOk = computeEigenvectors; + return *this; +} + +namespace internal { +/** + * \internal + * \brief Compute the eigendecomposition from a tridiagonal matrix + * + * \param[in,out] diag : On input, the diagonal of the matrix, on output the eigenvalues + * \param[in,out] subdiag : The subdiagonal part of the matrix (entries are modified during the decomposition) + * \param[in] maxIterations : the maximum number of iterations + * \param[in] computeEigenvectors : whether the eigenvectors have to be computed or not + * \param[out] eivec : The matrix to store the eigenvectors if computeEigenvectors==true. Must be allocated on input. + * \returns \c Success or \c NoConvergence + */ +template<typename MatrixType, typename DiagType, typename SubDiagType> +ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec) +{ + using std::abs; + + ComputationInfo info; + typedef typename MatrixType::Scalar Scalar; + + Index n = diag.size(); Index end = n-1; Index start = 0; Index iter = 0; // total number of iterations - + + typedef typename DiagType::RealScalar RealScalar; + const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); + const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon(); + while (end>0) { for (Index i = start; i<end; ++i) - if (internal::isMuchSmallerThan(abs(m_subdiag[i]),(abs(diag[i])+abs(diag[i+1])))) - m_subdiag[i] = 0; + if (internal::isMuchSmallerThan(abs(subdiag[i]),(abs(diag[i])+abs(diag[i+1])),precision) || abs(subdiag[i]) <= considerAsZero) + subdiag[i] = 0; // find the largest unreduced block - while (end>0 && m_subdiag[end-1]==0) + while (end>0 && subdiag[end-1]==RealScalar(0)) { end--; } @@ -445,51 +511,42 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> // if we spent too many iterations, we give up iter++; - if(iter > m_maxIterations * n) break; + if(iter > maxIterations * n) break; start = end - 1; - while (start>0 && m_subdiag[start-1]!=0) + while (start>0 && subdiag[start-1]!=0) start--; - internal::tridiagonal_qr_step(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); + internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n); } - - if (iter <= m_maxIterations * n) - m_info = Success; + if (iter <= maxIterations * n) + info = Success; else - m_info = NoConvergence; + info = NoConvergence; // Sort eigenvalues and corresponding vectors. // TODO make the sort optional ? // TODO use a better sort algorithm !! - if (m_info == Success) + if (info == Success) { for (Index i = 0; i < n-1; ++i) { Index k; - m_eivalues.segment(i,n-i).minCoeff(&k); + diag.segment(i,n-i).minCoeff(&k); if (k > 0) { - std::swap(m_eivalues[i], m_eivalues[k+i]); + std::swap(diag[i], diag[k+i]); if(computeEigenvectors) - m_eivec.col(i).swap(m_eivec.col(k+i)); + eivec.col(i).swap(eivec.col(k+i)); } } } - - // scale back the eigen values - m_eivalues *= scale; - - m_isInitialized = true; - m_eigenvectorsOk = computeEigenvectors; - return *this; + return info; } - - -namespace internal { 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); } }; @@ -499,21 +556,22 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 typedef typename SolverType::MatrixType MatrixType; typedef typename SolverType::RealVectorType VectorType; typedef typename SolverType::Scalar Scalar; - typedef typename MatrixType::Index Index; typedef typename SolverType::EigenvectorsType EigenvectorsType; + /** \internal * Computes the roots of the characteristic polynomial of \a m. * For numerical stability m.trace() should be near zero and to avoid over- or underflow m should be normalized. */ + EIGEN_DEVICE_FUNC static inline void computeRoots(const MatrixType& m, VectorType& roots) { - using std::sqrt; - using std::atan2; - using std::cos; - using std::sin; - const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0); - const Scalar s_sqrt3 = sqrt(Scalar(3.0)); + 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)/Scalar(3); + const Scalar s_sqrt3 = sqrt(Scalar(3)); // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The // eigenvalues are the roots to this equation, all guaranteed to be @@ -526,14 +584,12 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 // and in solving the equation for the roots in closed form. Scalar c2_over_3 = c2*s_inv3; Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3; - if(a_over_3<Scalar(0)) - a_over_3 = Scalar(0); + a_over_3 = numext::maxi(a_over_3, Scalar(0)); Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1)); Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b; - if(q<Scalar(0)) - q = Scalar(0); + q = numext::maxi(q, Scalar(0)); // Compute the eigenvalues by solving for the roots of the polynomial. Scalar rho = sqrt(a_over_3); @@ -546,6 +602,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta; } + EIGEN_DEVICE_FUNC static inline bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative) { using std::abs; @@ -565,6 +622,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 return true; } + EIGEN_DEVICE_FUNC static inline void run(SolverType& solver, const MatrixType& mat, int options) { eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows()); @@ -606,7 +664,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 Index k(0), l(2); if(d0 > d1) { - std::swap(k,l); + numext::swap(k,l); d0 = d1; } @@ -650,13 +708,15 @@ 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; typedef typename SolverType::EigenvectorsType EigenvectorsType; + EIGEN_DEVICE_FUNC static inline void computeRoots(const MatrixType& m, VectorType& roots) { using std::sqrt; @@ -666,11 +726,12 @@ 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; - using std::abs; - + EIGEN_USING_STD_MATH(sqrt); + EIGEN_USING_STD_MATH(abs); + eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows()); eigen_assert((options&~(EigVecMask|GenEigMask))==0 && (options&EigVecMask)!=EigVecMask @@ -680,14 +741,18 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2 EigenvectorsType& eivecs = solver.m_eivec; VectorType& eivals = solver.m_eivalues; - // map the matrix coefficients to [-1:1] to avoid over- and underflow. - Scalar scale = mat.cwiseAbs().maxCoeff(); - scale = (std::max)(scale,Scalar(1)); - MatrixType scaledMat = mat / scale; - + // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow. + Scalar shift = mat.trace() / Scalar(2); + MatrixType scaledMat = mat; + scaledMat.coeffRef(0,1) = mat.coeff(1,0); + scaledMat.diagonal().array() -= shift; + Scalar scale = scaledMat.cwiseAbs().maxCoeff(); + if(scale > Scalar(0)) + scaledMat /= scale; + // Compute the eigenvalues computeRoots(scaledMat,eivals); - + // compute the eigen vectors if(computeEigenvectors) { @@ -715,10 +780,11 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2 eivecs.col(0) << eivecs.col(1).unitOrthogonal(); } } - + // Rescale back to the original size. eivals *= scale; - + eivals.array() += shift; + solver.m_info = Success; solver.m_isInitialized = true; solver.m_eigenvectorsOk = computeEigenvectors; @@ -728,6 +794,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) { @@ -736,7 +803,8 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> } namespace internal { -template<typename RealScalar, typename Scalar, typename Index> +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; @@ -748,14 +816,14 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); // This explain the following, somewhat more complicated, version: RealScalar mu = diag[end]; - if(td==0) + if(td==RealScalar(0)) mu -= abs(e); else { RealScalar e2 = numext::abs2(subdiag[end-1]); RealScalar h = numext::hypot(td,e); - if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h); - else mu -= e2 / (td + (td>0 ? h : -h)); + if(e2==RealScalar(0)) mu -= (e / (td + (td>RealScalar(0) ? RealScalar(1) : RealScalar(-1)))) * (e / h); + else mu -= e2 / (td + (td>RealScalar(0) ? h : -h)); } RealScalar x = diag[start] - mu; @@ -788,7 +856,8 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta // apply the givens rotation to the unit matrix Q = Q * G if (matrixQ) { - Map<Matrix<Scalar,Dynamic,Dynamic,ColMajor> > q(matrixQ,n,n); + // FIXME if StorageOrder == RowMajor this operation is not very efficient + Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n); q.applyOnTheRight(k,k+1,rot); } } |