diff options
Diffstat (limited to 'Eigen/src/Eigenvalues/RealSchur.h')
-rw-r--r-- | Eigen/src/Eigenvalues/RealSchur.h | 52 |
1 files changed, 27 insertions, 25 deletions
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h index f9d49c6b7..92ff448ed 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -77,6 +77,8 @@ template<typename _MatrixType> class RealSchur }; typedef typename MatrixType::Scalar Scalar; typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; + typedef typename MatrixType::Index Index; + typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; @@ -91,7 +93,7 @@ template<typename _MatrixType> class RealSchur * * \sa compute() for an example. */ - RealSchur(int size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) + RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : m_matT(size, size), m_matU(size, size), m_workspaceVector(size), @@ -177,11 +179,11 @@ template<typename _MatrixType> class RealSchur typedef Matrix<Scalar,3,1> Vector3s; Scalar computeNormOfT(); - int findSmallSubdiagEntry(int iu, Scalar norm); - void splitOffTwoRows(int iu, Scalar exshift); - void computeShift(int iu, int iter, Scalar& exshift, Vector3s& shiftInfo); - void initFrancisQRStep(int il, int iu, const Vector3s& shiftInfo, int& im, Vector3s& firstHouseholderVector); - void performFrancisQRStep(int il, int im, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace); + Index findSmallSubdiagEntry(Index iu, Scalar norm); + void splitOffTwoRows(Index iu, Scalar exshift); + void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); + void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); + void performFrancisQRStep(Index il, Index im, Index iu, const Vector3s& firstHouseholderVector, Scalar* workspace); }; @@ -204,14 +206,14 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix) // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. // Rows il,...,iu is the part we are working on (the active window). // Rows iu+1,...,end are already brought in triangular form. - int iu = m_matU.cols() - 1; - int iter = 0; // iteration count + Index iu = m_matU.cols() - 1; + Index iter = 0; // iteration count Scalar exshift = 0.0; // sum of exceptional shifts Scalar norm = computeNormOfT(); while (iu >= 0) { - int il = findSmallSubdiagEntry(iu, norm); + Index il = findSmallSubdiagEntry(iu, norm); // Check for convergence if (il == iu) // One root found @@ -233,7 +235,7 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix) Vector3s firstHouseholderVector, shiftInfo; computeShift(iu, iter, exshift, shiftInfo); iter = iter + 1; // (Could check iteration count here.) - int im; + Index im; initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector); performFrancisQRStep(il, im, iu, firstHouseholderVector, workspace); } @@ -246,21 +248,21 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix) template<typename MatrixType> inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT() { - const int size = m_matU.cols(); + const Index size = m_matU.cols(); // FIXME to be efficient the following would requires a triangular reduxion code // Scalar norm = m_matT.upper().cwiseAbs().sum() // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum(); Scalar norm = 0.0; - for (int j = 0; j < size; ++j) - norm += m_matT.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum(); + for (Index j = 0; j < size; ++j) + norm += m_matT.row(j).segment(std::max(j-1,Index(0)), size-std::max(j-1,Index(0))).cwiseAbs().sum(); return norm; } /** \internal Look for single small sub-diagonal element and returns its index */ template<typename MatrixType> -inline int RealSchur<MatrixType>::findSmallSubdiagEntry(int iu, Scalar norm) +inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm) { - int res = iu; + Index res = iu; while (res > 0) { Scalar s = ei_abs(m_matT.coeff(res-1,res-1)) + ei_abs(m_matT.coeff(res,res)); @@ -275,9 +277,9 @@ inline int RealSchur<MatrixType>::findSmallSubdiagEntry(int iu, Scalar norm) /** \internal Update T given that rows iu-1 and iu decouple from the rest. */ template<typename MatrixType> -inline void RealSchur<MatrixType>::splitOffTwoRows(int iu, Scalar exshift) +inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, Scalar exshift) { - const int size = m_matU.cols(); + const Index size = m_matU.cols(); // The eigenvalues of the 2x2 matrix [a b; c d] are // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc @@ -307,7 +309,7 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(int iu, Scalar exshift) /** \internal Form shift in shiftInfo, and update exshift if an exceptional shift is performed. */ template<typename MatrixType> -inline void RealSchur<MatrixType>::computeShift(int iu, int iter, Scalar& exshift, Vector3s& shiftInfo) +inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo) { shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu); shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1); @@ -317,7 +319,7 @@ inline void RealSchur<MatrixType>::computeShift(int iu, int iter, Scalar& exshif if (iter == 10) { exshift += shiftInfo.coeff(0); - for (int i = 0; i <= iu; ++i) + for (Index i = 0; i <= iu; ++i) m_matT.coeffRef(i,i) -= shiftInfo.coeff(0); Scalar s = ei_abs(m_matT.coeff(iu,iu-1)) + ei_abs(m_matT.coeff(iu-1,iu-2)); shiftInfo.coeffRef(0) = Scalar(0.75) * s; @@ -338,7 +340,7 @@ inline void RealSchur<MatrixType>::computeShift(int iu, int iter, Scalar& exshif s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0); s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s; exshift += s; - for (int i = 0; i <= iu; ++i) + for (Index i = 0; i <= iu; ++i) m_matT.coeffRef(i,i) -= s; shiftInfo.setConstant(Scalar(0.964)); } @@ -347,7 +349,7 @@ inline void RealSchur<MatrixType>::computeShift(int iu, int iter, Scalar& exshif /** \internal Compute index im at which Francis QR step starts and the first Householder vector. */ template<typename MatrixType> -inline void RealSchur<MatrixType>::initFrancisQRStep(int il, int iu, const Vector3s& shiftInfo, int& im, Vector3s& firstHouseholderVector) +inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector) { Vector3s& v = firstHouseholderVector; // alias to save typing @@ -373,14 +375,14 @@ inline void RealSchur<MatrixType>::initFrancisQRStep(int il, int iu, const Vecto /** \internal Perform a Francis QR step involving rows il:iu and columns im:iu. */ template<typename MatrixType> -inline void RealSchur<MatrixType>::performFrancisQRStep(int il, int im, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace) +inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, const Vector3s& firstHouseholderVector, Scalar* workspace) { assert(im >= il); assert(im <= iu-2); - const int size = m_matU.cols(); + const Index size = m_matU.cols(); - for (int k = im; k <= iu-2; ++k) + for (Index k = im; k <= iu-2; ++k) { bool firstIteration = (k == im); @@ -422,7 +424,7 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(int il, int im, int iu, } // clean up pollution due to round-off errors - for (int i = im+2; i <= iu; ++i) + for (Index i = im+2; i <= iu; ++i) { m_matT.coeffRef(i,i-2) = Scalar(0); if (i > im+2) |