diff options
Diffstat (limited to 'extern/Eigen3/Eigen/src/Eigenvalues/RealSchur.h')
| -rw-r--r-- | extern/Eigen3/Eigen/src/Eigenvalues/RealSchur.h | 45 |
1 files changed, 33 insertions, 12 deletions
diff --git a/extern/Eigen3/Eigen/src/Eigenvalues/RealSchur.h b/extern/Eigen3/Eigen/src/Eigenvalues/RealSchur.h index 16d3875372e..17ea903f5f1 100644 --- a/extern/Eigen3/Eigen/src/Eigenvalues/RealSchur.h +++ b/extern/Eigen3/Eigen/src/Eigenvalues/RealSchur.h @@ -64,7 +64,7 @@ template<typename _MatrixType> class RealSchur }; typedef typename MatrixType::Scalar Scalar; typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; - typedef typename MatrixType::Index Index; + typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; @@ -80,7 +80,7 @@ template<typename _MatrixType> class RealSchur * * \sa compute() for an example. */ - RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) + explicit RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : m_matT(size, size), m_matU(size, size), m_workspaceVector(size), @@ -100,7 +100,8 @@ template<typename _MatrixType> class RealSchur * Example: \include RealSchur_RealSchur_MatrixType.cpp * Output: \verbinclude RealSchur_RealSchur_MatrixType.out */ - RealSchur(const MatrixType& matrix, bool computeU = true) + template<typename InputType> + explicit RealSchur(const EigenBase<InputType>& matrix, bool computeU = true) : m_matT(matrix.rows(),matrix.cols()), m_matU(matrix.rows(),matrix.cols()), m_workspaceVector(matrix.rows()), @@ -109,7 +110,7 @@ template<typename _MatrixType> class RealSchur m_matUisUptodate(false), m_maxIters(-1) { - compute(matrix, computeU); + compute(matrix.derived(), computeU); } /** \brief Returns the orthogonal matrix in the Schur decomposition. @@ -165,7 +166,8 @@ template<typename _MatrixType> class RealSchur * * \sa compute(const MatrixType&, bool, Index) */ - RealSchur& compute(const MatrixType& matrix, bool computeU = true); + template<typename InputType> + RealSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true); /** \brief Computes Schur decomposition of a Hessenberg matrix H = Z T Z^T * \param[in] matrixH Matrix in Hessenberg form H @@ -243,26 +245,45 @@ template<typename _MatrixType> class RealSchur template<typename MatrixType> -RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) +template<typename InputType> +RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU) { + const Scalar considerAsZero = (std::numeric_limits<Scalar>::min)(); + eigen_assert(matrix.cols() == matrix.rows()); Index maxIters = m_maxIters; if (maxIters == -1) maxIters = m_maxIterationsPerRow * matrix.rows(); + Scalar scale = matrix.derived().cwiseAbs().maxCoeff(); + if(scale<considerAsZero) + { + m_matT.setZero(matrix.rows(),matrix.cols()); + if(computeU) + m_matU.setIdentity(matrix.rows(),matrix.cols()); + m_info = Success; + m_isInitialized = true; + m_matUisUptodate = computeU; + return *this; + } + // Step 1. Reduce to Hessenberg form - m_hess.compute(matrix); + m_hess.compute(matrix.derived()/scale); // Step 2. Reduce to real Schur form computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU); + + m_matT *= scale; return *this; } template<typename MatrixType> template<typename HessMatrixType, typename OrthMatrixType> RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) -{ - m_matT = matrixH; +{ + using std::abs; + + m_matT = matrixH; if(computeU) m_matU = matrixQ; @@ -282,7 +303,7 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMa Scalar exshift(0); // sum of exceptional shifts Scalar norm = computeNormOfT(); - if(norm!=0) + if(norm!=Scalar(0)) { while (iu >= 0) { @@ -306,7 +327,7 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMa else // No convergence yet { // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG ) - Vector3s firstHouseholderVector(0,0,0), shiftInfo; + Vector3s firstHouseholderVector = Vector3s::Zero(), shiftInfo; computeShift(iu, iter, exshift, shiftInfo); iter = iter + 1; totalIter = totalIter + 1; @@ -343,7 +364,7 @@ inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT() /** \internal Look for single small sub-diagonal element and returns its index */ template<typename MatrixType> -inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu) +inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu) { using std::abs; Index res = iu; |