Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/Eigen/src/Eigenvalues/CMakeLists.txt b/Eigen/src/Eigenvalues/CMakeLists.txt
deleted file mode 100644
index 193e026..0000000
--- a/Eigen/src/Eigenvalues/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_EIGENVALUES_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_EIGENVALUES_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigenvalues COMPONENT Devel
- )
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index 417c729..dc5fae0 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -60,7 +60,7 @@
/** \brief Scalar type for matrices of type #MatrixType. */
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef typename MatrixType::Index Index;
+ typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
/** \brief Complex scalar type for #MatrixType.
*
@@ -104,7 +104,7 @@
* according to the specified problem \a size.
* \sa ComplexEigenSolver()
*/
- ComplexEigenSolver(Index size)
+ explicit ComplexEigenSolver(Index size)
: m_eivec(size, size),
m_eivalues(size),
m_schur(size),
@@ -122,7 +122,8 @@
*
* This constructor calls compute() to compute the eigendecomposition.
*/
- ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
+ template<typename InputType>
+ explicit ComplexEigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
: m_eivec(matrix.rows(),matrix.cols()),
m_eivalues(matrix.cols()),
m_schur(matrix.rows()),
@@ -130,7 +131,7 @@
m_eigenvectorsOk(false),
m_matX(matrix.rows(),matrix.cols())
{
- compute(matrix, computeEigenvectors);
+ compute(matrix.derived(), computeEigenvectors);
}
/** \brief Returns the eigenvectors of given matrix.
@@ -208,7 +209,8 @@
* Example: \include ComplexEigenSolver_compute.cpp
* Output: \verbinclude ComplexEigenSolver_compute.out
*/
- ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
+ template<typename InputType>
+ ComplexEigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
/** \brief Reports whether previous computation was successful.
*
@@ -248,14 +250,15 @@
EigenvectorType m_matX;
private:
- void doComputeEigenvectors(const RealScalar& matrixnorm);
+ void doComputeEigenvectors(RealScalar matrixnorm);
void sortEigenvalues(bool computeEigenvectors);
};
template<typename MatrixType>
+template<typename InputType>
ComplexEigenSolver<MatrixType>&
-ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
+ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors)
{
check_template_parameters();
@@ -264,13 +267,13 @@
// Do a complex Schur decomposition, A = U T U^*
// The eigenvalues are on the diagonal of T.
- m_schur.compute(matrix, computeEigenvectors);
+ m_schur.compute(matrix.derived(), computeEigenvectors);
if(m_schur.info() == Success)
{
m_eivalues = m_schur.matrixT().diagonal();
if(computeEigenvectors)
- doComputeEigenvectors(matrix.norm());
+ doComputeEigenvectors(m_schur.matrixT().norm());
sortEigenvalues(computeEigenvectors);
}
@@ -281,10 +284,12 @@
template<typename MatrixType>
-void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm)
+void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
{
const Index n = m_eivalues.size();
+ matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)());
+
// Compute X such that T = X D X^(-1), where D is the diagonal of T.
// The matrix X is unit triangular.
m_matX = EigenvectorType::Zero(n, n);
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index 89e6cad..7f38919 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -63,7 +63,7 @@
/** \brief Scalar type for matrices of type \p _MatrixType. */
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef typename MatrixType::Index Index;
+ typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
/** \brief Complex scalar type for \p _MatrixType.
*
@@ -91,7 +91,7 @@
*
* \sa compute() for an example.
*/
- ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
+ explicit ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
: m_matT(size,size),
m_matU(size,size),
m_hess(size),
@@ -109,7 +109,8 @@
*
* \sa matrixT() and matrixU() for examples.
*/
- ComplexSchur(const MatrixType& matrix, bool computeU = true)
+ template<typename InputType>
+ explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true)
: m_matT(matrix.rows(),matrix.cols()),
m_matU(matrix.rows(),matrix.cols()),
m_hess(matrix.rows()),
@@ -117,7 +118,7 @@
m_matUisUptodate(false),
m_maxIters(-1)
{
- compute(matrix, computeU);
+ compute(matrix.derived(), computeU);
}
/** \brief Returns the unitary matrix in the Schur decomposition.
@@ -186,7 +187,8 @@
*
* \sa compute(const MatrixType&, bool, Index)
*/
- ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
+ template<typename InputType>
+ ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);
/** \brief Compute Schur decomposition from a given Hessenberg matrix
* \param[in] matrixH Matrix in Hessenberg form H
@@ -313,14 +315,15 @@
template<typename MatrixType>
-ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
+template<typename InputType>
+ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU)
{
m_matUisUptodate = false;
eigen_assert(matrix.cols() == matrix.rows());
if(matrix.cols() == 1)
{
- m_matT = matrix.template cast<ComplexScalar>();
+ m_matT = matrix.derived().template cast<ComplexScalar>();
if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
m_info = Success;
m_isInitialized = true;
@@ -328,7 +331,7 @@
return *this;
}
- internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
+ internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix.derived(), computeU);
computeFromHessenberg(m_matT, m_matU, computeU);
return *this;
}
diff --git a/Eigen/src/Eigenvalues/ComplexSchur_MKL.h b/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h
similarity index 65%
rename from Eigen/src/Eigenvalues/ComplexSchur_MKL.h
rename to Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h
index 91496ae..4980a3e 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur_MKL.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h
@@ -25,27 +25,24 @@
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
********************************************************************************
- * Content : Eigen bindings to Intel(R) MKL
+ * Content : Eigen bindings to LAPACKe
* Complex Schur needed to complex unsymmetrical eigenvalues/eigenvectors.
********************************************************************************
*/
-#ifndef EIGEN_COMPLEX_SCHUR_MKL_H
-#define EIGEN_COMPLEX_SCHUR_MKL_H
-
-#include "Eigen/src/Core/util/MKL_support.h"
+#ifndef EIGEN_COMPLEX_SCHUR_LAPACKE_H
+#define EIGEN_COMPLEX_SCHUR_LAPACKE_H
namespace Eigen {
-/** \internal Specialization for the data types supported by MKL */
+/** \internal Specialization for the data types supported by LAPACKe */
-#define EIGEN_MKL_SCHUR_COMPLEX(EIGTYPE, MKLTYPE, MKLPREFIX, MKLPREFIX_U, EIGCOLROW, MKLCOLROW) \
-template<> inline \
+#define EIGEN_LAPACKE_SCHUR_COMPLEX(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, LAPACKE_PREFIX_U, EIGCOLROW, LAPACKE_COLROW) \
+template<> template<typename InputType> inline \
ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \
-ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW>& matrix, bool computeU) \
+ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBase<InputType>& matrix, bool computeU) \
{ \
typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> MatrixType; \
- typedef MatrixType::Scalar Scalar; \
typedef MatrixType::RealScalar RealScalar; \
typedef std::complex<RealScalar> ComplexScalar; \
\
@@ -54,25 +51,25 @@
m_matUisUptodate = false; \
if(matrix.cols() == 1) \
{ \
- m_matT = matrix.cast<ComplexScalar>(); \
+ m_matT = matrix.derived().template cast<ComplexScalar>(); \
if(computeU) m_matU = ComplexMatrixType::Identity(1,1); \
m_info = Success; \
m_isInitialized = true; \
m_matUisUptodate = computeU; \
return *this; \
} \
- lapack_int n = matrix.cols(), sdim, info; \
- lapack_int lda = matrix.outerStride(); \
- lapack_int matrix_order = MKLCOLROW; \
+ lapack_int n = internal::convert_index<lapack_int>(matrix.cols()), sdim, info; \
+ lapack_int matrix_order = LAPACKE_COLROW; \
char jobvs, sort='N'; \
- LAPACK_##MKLPREFIX_U##_SELECT1 select = 0; \
+ LAPACK_##LAPACKE_PREFIX_U##_SELECT1 select = 0; \
jobvs = (computeU) ? 'V' : 'N'; \
m_matU.resize(n, n); \
- lapack_int ldvs = m_matU.outerStride(); \
+ lapack_int ldvs = internal::convert_index<lapack_int>(m_matU.outerStride()); \
m_matT = matrix; \
+ lapack_int lda = internal::convert_index<lapack_int>(m_matT.outerStride()); \
Matrix<EIGTYPE, Dynamic, Dynamic> w; \
w.resize(n, 1);\
- info = LAPACKE_##MKLPREFIX##gees( matrix_order, jobvs, sort, select, n, (MKLTYPE*)m_matT.data(), lda, &sdim, (MKLTYPE*)w.data(), (MKLTYPE*)m_matU.data(), ldvs ); \
+ info = LAPACKE_##LAPACKE_PREFIX##gees( matrix_order, jobvs, sort, select, n, (LAPACKE_TYPE*)m_matT.data(), lda, &sdim, (LAPACKE_TYPE*)w.data(), (LAPACKE_TYPE*)m_matU.data(), ldvs ); \
if(info == 0) \
m_info = Success; \
else \
@@ -84,11 +81,11 @@
\
}
-EIGEN_MKL_SCHUR_COMPLEX(dcomplex, MKL_Complex16, z, Z, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_SCHUR_COMPLEX(scomplex, MKL_Complex8, c, C, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_SCHUR_COMPLEX(dcomplex, MKL_Complex16, z, Z, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_SCHUR_COMPLEX(scomplex, MKL_Complex8, c, C, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_SCHUR_COMPLEX(dcomplex, lapack_complex_double, z, Z, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_SCHUR_COMPLEX(scomplex, lapack_complex_float, c, C, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_SCHUR_COMPLEX(dcomplex, lapack_complex_double, z, Z, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_SCHUR_COMPLEX(scomplex, lapack_complex_float, c, C, RowMajor, LAPACK_ROW_MAJOR)
} // end namespace Eigen
-#endif // EIGEN_COMPLEX_SCHUR_MKL_H
+#endif // EIGEN_COMPLEX_SCHUR_LAPACKE_H
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index 20c59a7..f205b18 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -79,7 +79,7 @@
/** \brief Scalar type for matrices of type #MatrixType. */
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef typename MatrixType::Index Index;
+ typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
/** \brief Complex scalar type for #MatrixType.
*
@@ -110,7 +110,7 @@
*
* \sa compute() for an example.
*/
- EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
+ EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
/** \brief Default constructor with memory preallocation
*
@@ -118,7 +118,7 @@
* according to the specified problem \a size.
* \sa EigenSolver()
*/
- EigenSolver(Index size)
+ explicit EigenSolver(Index size)
: m_eivec(size, size),
m_eivalues(size),
m_isInitialized(false),
@@ -143,7 +143,8 @@
*
* \sa compute()
*/
- EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
+ template<typename InputType>
+ explicit EigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
: m_eivec(matrix.rows(), matrix.cols()),
m_eivalues(matrix.cols()),
m_isInitialized(false),
@@ -152,7 +153,7 @@
m_matT(matrix.rows(), matrix.cols()),
m_tmp(matrix.cols())
{
- compute(matrix, computeEigenvectors);
+ compute(matrix.derived(), computeEigenvectors);
}
/** \brief Returns the eigenvectors of given matrix.
@@ -273,12 +274,14 @@
* Example: \include EigenSolver_compute.cpp
* Output: \verbinclude EigenSolver_compute.out
*/
- EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
+ template<typename InputType>
+ EigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
+ /** \returns NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise. */
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
- return m_realSchur.info();
+ return m_info;
}
/** \brief Sets the maximum number of iterations allowed. */
@@ -309,6 +312,7 @@
EigenvalueType m_eivalues;
bool m_isInitialized;
bool m_eigenvectorsOk;
+ ComputationInfo m_info;
RealSchur<MatrixType> m_realSchur;
MatrixType m_matT;
@@ -320,11 +324,12 @@
MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
{
eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
+ const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
Index n = m_eivalues.rows();
MatrixType matD = MatrixType::Zero(n,n);
for (Index i=0; i<n; ++i)
{
- if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i))))
+ if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)), precision))
matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
else
{
@@ -341,11 +346,12 @@
{
eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
+ const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
Index n = m_eivec.cols();
EigenvectorsType matV(n,n);
for (Index j=0; j<n; ++j)
{
- if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j))) || j+1==n)
+ if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j)), precision) || j+1==n)
{
// we have a real eigen value
matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
@@ -368,19 +374,23 @@
}
template<typename MatrixType>
+template<typename InputType>
EigenSolver<MatrixType>&
-EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
+EigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors)
{
check_template_parameters();
using std::sqrt;
using std::abs;
+ using numext::isfinite;
eigen_assert(matrix.cols() == matrix.rows());
// Reduce to real Schur form.
- m_realSchur.compute(matrix, computeEigenvectors);
+ m_realSchur.compute(matrix.derived(), computeEigenvectors);
+
+ m_info = m_realSchur.info();
- if (m_realSchur.info() == Success)
+ if (m_info == Success)
{
m_matT = m_realSchur.matrixT();
if (computeEigenvectors)
@@ -394,14 +404,40 @@
if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
{
m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
+ if(!(isfinite)(m_eivalues.coeffRef(i)))
+ {
+ m_isInitialized = true;
+ m_eigenvectorsOk = false;
+ m_info = NumericalIssue;
+ return *this;
+ }
++i;
}
else
{
Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
- Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
+ Scalar z;
+ // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
+ // without overflow
+ {
+ Scalar t0 = m_matT.coeff(i+1, i);
+ Scalar t1 = m_matT.coeff(i, i+1);
+ Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1)));
+ t0 /= maxval;
+ t1 /= maxval;
+ Scalar p0 = p/maxval;
+ z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
+ }
+
m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
+ if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1))))
+ {
+ m_isInitialized = true;
+ m_eigenvectorsOk = false;
+ m_info = NumericalIssue;
+ return *this;
+ }
i += 2;
}
}
@@ -417,26 +453,6 @@
return *this;
}
-// Complex scalar division.
-template<typename Scalar>
-std::complex<Scalar> cdiv(const Scalar& xr, const Scalar& xi, const Scalar& yr, const Scalar& yi)
-{
- using std::abs;
- Scalar r,d;
- if (abs(yr) > abs(yi))
- {
- r = yi/yr;
- d = yr + r*yi;
- return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
- }
- else
- {
- r = yr/yi;
- d = yi + r*yr;
- return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
- }
-}
-
template<typename MatrixType>
void EigenSolver<MatrixType>::doComputeEigenvectors()
@@ -453,7 +469,7 @@
}
// Backsubstitute to find vectors of upper triangular form
- if (norm == 0.0)
+ if (norm == Scalar(0))
{
return;
}
@@ -469,13 +485,13 @@
Scalar lastr(0), lastw(0);
Index l = n;
- m_matT.coeffRef(n,n) = 1.0;
+ m_matT.coeffRef(n,n) = Scalar(1);
for (Index i = n-1; i >= 0; i--)
{
Scalar w = m_matT.coeff(i,i) - p;
Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
- if (m_eivalues.coeff(i).imag() < 0.0)
+ if (m_eivalues.coeff(i).imag() < Scalar(0))
{
lastw = w;
lastr = r;
@@ -483,9 +499,9 @@
else
{
l = i;
- if (m_eivalues.coeff(i).imag() == 0.0)
+ if (m_eivalues.coeff(i).imag() == Scalar(0))
{
- if (w != 0.0)
+ if (w != Scalar(0))
m_matT.coeffRef(i,n) = -r / w;
else
m_matT.coeffRef(i,n) = -r / (eps * norm);
@@ -523,19 +539,19 @@
}
else
{
- std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
+ ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q);
m_matT.coeffRef(n-1,n-1) = numext::real(cc);
m_matT.coeffRef(n-1,n) = numext::imag(cc);
}
- m_matT.coeffRef(n,n-1) = 0.0;
- m_matT.coeffRef(n,n) = 1.0;
+ m_matT.coeffRef(n,n-1) = Scalar(0);
+ m_matT.coeffRef(n,n) = Scalar(1);
for (Index i = n-2; i >= 0; i--)
{
Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
Scalar w = m_matT.coeff(i,i) - p;
- if (m_eivalues.coeff(i).imag() < 0.0)
+ if (m_eivalues.coeff(i).imag() < Scalar(0))
{
lastw = w;
lastra = ra;
@@ -546,7 +562,7 @@
l = i;
if (m_eivalues.coeff(i).imag() == RealScalar(0))
{
- std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
+ ComplexScalar cc = ComplexScalar(-ra,-sa) / ComplexScalar(w,q);
m_matT.coeffRef(i,n-1) = numext::real(cc);
m_matT.coeffRef(i,n) = numext::imag(cc);
}
@@ -557,10 +573,10 @@
Scalar y = m_matT.coeff(i+1,i);
Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
- if ((vr == 0.0) && (vi == 0.0))
+ if ((vr == Scalar(0)) && (vi == Scalar(0)))
vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
- std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
+ ComplexScalar cc = ComplexScalar(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra) / ComplexScalar(vr,vi);
m_matT.coeffRef(i,n-1) = numext::real(cc);
m_matT.coeffRef(i,n) = numext::imag(cc);
if (abs(x) > (abs(lastw) + abs(q)))
@@ -570,15 +586,14 @@
}
else
{
- cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
+ cc = ComplexScalar(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n)) / ComplexScalar(lastw,q);
m_matT.coeffRef(i+1,n-1) = numext::real(cc);
m_matT.coeffRef(i+1,n) = numext::imag(cc);
}
}
// Overflow control
- using std::max;
- Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
+ Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
if ((eps * t) * t > Scalar(1))
m_matT.block(i, n-1, size-i, 2) /= t;
@@ -590,7 +605,7 @@
}
else
{
- eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
+ eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen
}
}
diff --git a/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h b/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h
index 956e80d..87d789b 100644
--- a/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h
+++ b/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h
@@ -1,8 +1,9 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -72,7 +73,7 @@
/** \brief Scalar type for matrices of type #MatrixType. */
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef typename MatrixType::Index Index;
+ typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
/** \brief Complex scalar type for #MatrixType.
*
@@ -89,7 +90,7 @@
*/
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;
- /** \brief Type for vector of complex scalar values eigenvalues as returned by betas().
+ /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas().
*
* This is a column vector with entries of type #ComplexScalar.
* The length of the vector is the size of #MatrixType.
@@ -114,7 +115,14 @@
*
* \sa compute() for an example.
*/
- GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {}
+ GeneralizedEigenSolver()
+ : m_eivec(),
+ m_alphas(),
+ m_betas(),
+ m_valuesOkay(false),
+ m_vectorsOkay(false),
+ m_realQZ()
+ {}
/** \brief Default constructor with memory preallocation
*
@@ -122,14 +130,13 @@
* according to the specified problem \a size.
* \sa GeneralizedEigenSolver()
*/
- GeneralizedEigenSolver(Index size)
+ explicit GeneralizedEigenSolver(Index size)
: m_eivec(size, size),
m_alphas(size),
m_betas(size),
- m_isInitialized(false),
- m_eigenvectorsOk(false),
+ m_valuesOkay(false),
+ m_vectorsOkay(false),
m_realQZ(size),
- m_matS(size, size),
m_tmp(size)
{}
@@ -149,10 +156,9 @@
: m_eivec(A.rows(), A.cols()),
m_alphas(A.cols()),
m_betas(A.cols()),
- m_isInitialized(false),
- m_eigenvectorsOk(false),
+ m_valuesOkay(false),
+ m_vectorsOkay(false),
m_realQZ(A.cols()),
- m_matS(A.rows(), A.cols()),
m_tmp(A.cols())
{
compute(A, B, computeEigenvectors);
@@ -160,22 +166,20 @@
/* \brief Returns the computed generalized eigenvectors.
*
- * \returns %Matrix whose columns are the (possibly complex) eigenvectors.
+ * \returns %Matrix whose columns are the (possibly complex) right eigenvectors.
+ * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
*
* \pre Either the constructor
* GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
* compute(const MatrixType&, const MatrixType& bool) has been called before, and
* \p computeEigenvectors was set to true (the default).
*
- * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
- * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The
- * eigenvectors are normalized to have (Euclidean) norm equal to one. The
- * matrix returned by this function is the matrix \f$ V \f$ in the
- * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists.
- *
* \sa eigenvalues()
*/
-// EigenvectorsType eigenvectors() const;
+ EigenvectorsType eigenvectors() const {
+ eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.");
+ return m_eivec;
+ }
/** \brief Returns an expression of the computed generalized eigenvalues.
*
@@ -197,7 +201,7 @@
*/
EigenvalueType eigenvalues() const
{
- eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
+ eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
return EigenvalueType(m_alphas,m_betas);
}
@@ -208,7 +212,7 @@
* \sa betas(), eigenvalues() */
ComplexVectorType alphas() const
{
- eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
+ eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
return m_alphas;
}
@@ -219,7 +223,7 @@
* \sa alphas(), eigenvalues() */
VectorType betas() const
{
- eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
+ eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
return m_betas;
}
@@ -250,7 +254,7 @@
ComputationInfo info() const
{
- eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
+ eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
return m_realQZ.info();
}
@@ -270,29 +274,14 @@
EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
}
- MatrixType m_eivec;
+ EigenvectorsType m_eivec;
ComplexVectorType m_alphas;
VectorType m_betas;
- bool m_isInitialized;
- bool m_eigenvectorsOk;
+ bool m_valuesOkay, m_vectorsOkay;
RealQZ<MatrixType> m_realQZ;
- MatrixType m_matS;
-
- typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
- ColumnVectorType m_tmp;
+ ComplexVectorType m_tmp;
};
-//template<typename MatrixType>
-//typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
-//{
-// eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
-// eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
-// Index n = m_eivec.cols();
-// EigenvectorsType matV(n,n);
-// // TODO
-// return matV;
-//}
-
template<typename MatrixType>
GeneralizedEigenSolver<MatrixType>&
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
@@ -302,46 +291,125 @@
using std::sqrt;
using std::abs;
eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
-
+ Index size = A.cols();
+ m_valuesOkay = false;
+ m_vectorsOkay = false;
// Reduce to generalized real Schur form:
// A = Q S Z and B = Q T Z
m_realQZ.compute(A, B, computeEigenvectors);
-
if (m_realQZ.info() == Success)
{
- m_matS = m_realQZ.matrixS();
+ // Resize storage
+ m_alphas.resize(size);
+ m_betas.resize(size);
if (computeEigenvectors)
- m_eivec = m_realQZ.matrixZ().transpose();
-
- // Compute eigenvalues from matS
- m_alphas.resize(A.cols());
- m_betas.resize(A.cols());
- Index i = 0;
- while (i < A.cols())
{
- if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
+ m_eivec.resize(size,size);
+ m_tmp.resize(size);
+ }
+
+ // Aliases:
+ Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size);
+ ComplexVectorType &cv = m_tmp;
+ const MatrixType &mS = m_realQZ.matrixS();
+ const MatrixType &mT = m_realQZ.matrixT();
+
+ Index i = 0;
+ while (i < size)
+ {
+ if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0))
{
- m_alphas.coeffRef(i) = m_matS.coeff(i, i);
- m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i);
+ // Real eigenvalue
+ m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
+ m_betas.coeffRef(i) = mT.diagonal().coeff(i);
+ if (computeEigenvectors)
+ {
+ v.setConstant(Scalar(0.0));
+ v.coeffRef(i) = Scalar(1.0);
+ // For singular eigenvalues do nothing more
+ if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
+ {
+ // Non-singular eigenvalue
+ const Scalar alpha = real(m_alphas.coeffRef(i));
+ const Scalar beta = m_betas.coeffRef(i);
+ for (Index j = i-1; j >= 0; j--)
+ {
+ const Index st = j+1;
+ const Index sz = i-j;
+ if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
+ {
+ // 2x2 block
+ Matrix<Scalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( v.segment(st,sz) );
+ Matrix<Scalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
+ v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
+ j--;
+ }
+ else
+ {
+ v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j));
+ }
+ }
+ }
+ m_eivec.col(i).real().noalias() = m_realQZ.matrixZ().transpose() * v;
+ m_eivec.col(i).real().normalize();
+ m_eivec.col(i).imag().setConstant(0);
+ }
++i;
}
else
{
- Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1));
- Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
- m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z);
- m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z);
+ // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
+ // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00):
- m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i);
- m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i);
+ // T = [a 0]
+ // [0 b]
+ RealScalar a = mT.diagonal().coeff(i),
+ b = mT.diagonal().coeff(i+1);
+ const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
+
+ // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
+ Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
+
+ Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
+ Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
+ const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
+ m_alphas.coeffRef(i) = conj(alpha);
+ m_alphas.coeffRef(i+1) = alpha;
+
+ if (computeEigenvectors) {
+ // Compute eigenvector in position (i+1) and then position (i) is just the conjugate
+ cv.setZero();
+ cv.coeffRef(i+1) = Scalar(1.0);
+ // here, the "static_cast" workaound expression template issues.
+ cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
+ / (static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i));
+ for (Index j = i-1; j >= 0; j--)
+ {
+ const Index st = j+1;
+ const Index sz = i+1-j;
+ if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
+ {
+ // 2x2 block
+ Matrix<ComplexScalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( cv.segment(st,sz) );
+ Matrix<ComplexScalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
+ cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
+ j--;
+ } else {
+ cv.coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
+ / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j)));
+ }
+ }
+ m_eivec.col(i+1).noalias() = (m_realQZ.matrixZ().transpose() * cv);
+ m_eivec.col(i+1).normalize();
+ m_eivec.col(i) = m_eivec.col(i+1).conjugate();
+ }
i += 2;
}
}
+
+ m_valuesOkay = true;
+ m_vectorsOkay = computeEigenvectors;
}
-
- m_isInitialized = true;
- m_eigenvectorsOk = false;//computeEigenvectors;
-
return *this;
}
diff --git a/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h
index 07bf1ea..5f6bb82 100644
--- a/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h
@@ -50,7 +50,6 @@
typedef SelfAdjointEigenSolver<_MatrixType> Base;
public:
- typedef typename Base::Index Index;
typedef _MatrixType MatrixType;
/** \brief Default constructor for fixed-size matrices.
@@ -74,7 +73,7 @@
*
* \sa compute() for an example
*/
- GeneralizedSelfAdjointEigenSolver(Index size)
+ explicit GeneralizedSelfAdjointEigenSolver(Index size)
: Base(size)
{}
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
index 3db0c01..f647f69 100644
--- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -71,7 +71,7 @@
/** \brief Scalar type for matrices of type #MatrixType. */
typedef typename MatrixType::Scalar Scalar;
- typedef typename MatrixType::Index Index;
+ typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
/** \brief Type for vector of Householder coefficients.
*
@@ -97,7 +97,7 @@
*
* \sa compute() for an example.
*/
- HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
+ explicit HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
: m_matrix(size,size),
m_temp(size),
m_isInitialized(false)
@@ -115,8 +115,9 @@
*
* \sa matrixH() for an example.
*/
- HessenbergDecomposition(const MatrixType& matrix)
- : m_matrix(matrix),
+ template<typename InputType>
+ explicit HessenbergDecomposition(const EigenBase<InputType>& matrix)
+ : m_matrix(matrix.derived()),
m_temp(matrix.rows()),
m_isInitialized(false)
{
@@ -147,9 +148,10 @@
* Example: \include HessenbergDecomposition_compute.cpp
* Output: \verbinclude HessenbergDecomposition_compute.out
*/
- HessenbergDecomposition& compute(const MatrixType& matrix)
+ template<typename InputType>
+ HessenbergDecomposition& compute(const EigenBase<InputType>& matrix)
{
- m_matrix = matrix;
+ m_matrix = matrix.derived();
if(matrix.rows()<2)
{
m_isInitialized = true;
@@ -337,7 +339,6 @@
template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
: public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
{
- typedef typename MatrixType::Index Index;
public:
/** \brief Constructor.
*
diff --git a/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h b/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
index 4fec8af..e4e4260 100644
--- a/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
+++ b/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
@@ -66,7 +66,6 @@
inline typename MatrixBase<Derived>::EigenvaluesReturnType
MatrixBase<Derived>::eigenvalues() const
{
- typedef typename internal::traits<Derived>::Scalar Scalar;
return internal::eigenvalues_selector<Derived, NumTraits<Scalar>::IsComplex>::run(derived());
}
@@ -88,7 +87,6 @@
inline typename SelfAdjointView<MatrixType, UpLo>::EigenvaluesReturnType
SelfAdjointView<MatrixType, UpLo>::eigenvalues() const
{
- typedef typename SelfAdjointView<MatrixType, UpLo>::PlainObject PlainObject;
PlainObject thisAsMatrix(*this);
return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix, false).eigenvalues();
}
diff --git a/Eigen/src/Eigenvalues/RealQZ.h b/Eigen/src/Eigenvalues/RealQZ.h
index aa3833e..b3a910d 100644
--- a/Eigen/src/Eigenvalues/RealQZ.h
+++ b/Eigen/src/Eigenvalues/RealQZ.h
@@ -67,7 +67,7 @@
};
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;
@@ -83,7 +83,7 @@
*
* \sa compute() for an example.
*/
- RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) :
+ explicit RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) :
m_S(size, size),
m_T(size, size),
m_Q(size, size),
@@ -276,7 +276,7 @@
/** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */
template<typename MatrixType>
- inline typename MatrixType::Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu)
+ inline Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu)
{
using std::abs;
Index res = iu;
@@ -294,7 +294,7 @@
/** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1) */
template<typename MatrixType>
- inline typename MatrixType::Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l)
+ inline Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l)
{
using std::abs;
Index res = l;
@@ -315,8 +315,8 @@
const Index dim=m_S.cols();
if (abs(m_S.coeff(i+1,i))==Scalar(0))
return;
- Index z = findSmallDiagEntry(i,i+1);
- if (z==i-1)
+ Index j = findSmallDiagEntry(i,i+1);
+ if (j==i-1)
{
// block of (S T^{-1})
Matrix2s STi = m_T.template block<2,2>(i,i).template triangularView<Upper>().
@@ -352,7 +352,7 @@
}
else
{
- pushDownZero(z,i,i+1);
+ pushDownZero(j,i,i+1);
}
}
@@ -552,7 +552,6 @@
m_T.coeffRef(l,l-1) = Scalar(0.0);
}
-
template<typename MatrixType>
RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ)
{
@@ -616,6 +615,37 @@
}
// check if we converged before reaching iterations limit
m_info = (local_iter<m_maxIters) ? Success : NoConvergence;
+
+ // For each non triangular 2x2 diagonal block of S,
+ // reduce the respective 2x2 diagonal block of T to positive diagonal form using 2x2 SVD.
+ // This step is not mandatory for QZ, but it does help further extraction of eigenvalues/eigenvectors,
+ // and is in par with Lapack/Matlab QZ.
+ if(m_info==Success)
+ {
+ for(Index i=0; i<dim-1; ++i)
+ {
+ if(m_S.coeff(i+1, i) != Scalar(0))
+ {
+ JacobiRotation<Scalar> j_left, j_right;
+ internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right);
+
+ // Apply resulting Jacobi rotations
+ m_S.applyOnTheLeft(i,i+1,j_left);
+ m_S.applyOnTheRight(i,i+1,j_right);
+ m_T.applyOnTheLeft(i,i+1,j_left);
+ m_T.applyOnTheRight(i,i+1,j_right);
+ m_T(i+1,i) = m_T(i,i+1) = Scalar(0);
+
+ if(m_computeQZ) {
+ m_Q.applyOnTheRight(i,i+1,j_left.transpose());
+ m_Z.applyOnTheLeft(i,i+1,j_right.transpose());
+ }
+
+ i++;
+ }
+ }
+ }
+
return *this;
} // end compute
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index 16d3875..17ea903 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -64,7 +64,7 @@
};
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 @@
*
* \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 @@
* 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 @@
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 @@
*
* \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>
-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 @@
Scalar exshift(0); // sum of exceptional shifts
Scalar norm = computeNormOfT();
- if(norm!=0)
+ if(norm!=Scalar(0))
{
while (iu >= 0)
{
@@ -306,7 +327,7 @@
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 @@
/** \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;
diff --git a/Eigen/src/Eigenvalues/RealSchur_MKL.h b/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h
similarity index 64%
rename from Eigen/src/Eigenvalues/RealSchur_MKL.h
rename to Eigen/src/Eigenvalues/RealSchur_LAPACKE.h
index ad97364..2c22517 100644
--- a/Eigen/src/Eigenvalues/RealSchur_MKL.h
+++ b/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h
@@ -25,43 +25,37 @@
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
********************************************************************************
- * Content : Eigen bindings to Intel(R) MKL
+ * Content : Eigen bindings to LAPACKe
* Real Schur needed to real unsymmetrical eigenvalues/eigenvectors.
********************************************************************************
*/
-#ifndef EIGEN_REAL_SCHUR_MKL_H
-#define EIGEN_REAL_SCHUR_MKL_H
-
-#include "Eigen/src/Core/util/MKL_support.h"
+#ifndef EIGEN_REAL_SCHUR_LAPACKE_H
+#define EIGEN_REAL_SCHUR_LAPACKE_H
namespace Eigen {
-/** \internal Specialization for the data types supported by MKL */
+/** \internal Specialization for the data types supported by LAPACKe */
-#define EIGEN_MKL_SCHUR_REAL(EIGTYPE, MKLTYPE, MKLPREFIX, MKLPREFIX_U, EIGCOLROW, MKLCOLROW) \
-template<> inline \
+#define EIGEN_LAPACKE_SCHUR_REAL(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, LAPACKE_PREFIX_U, EIGCOLROW, LAPACKE_COLROW) \
+template<> template<typename InputType> inline \
RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \
-RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW>& matrix, bool computeU) \
+RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBase<InputType>& matrix, bool computeU) \
{ \
- typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> MatrixType; \
- typedef MatrixType::Scalar Scalar; \
- typedef MatrixType::RealScalar RealScalar; \
-\
eigen_assert(matrix.cols() == matrix.rows()); \
\
- lapack_int n = matrix.cols(), sdim, info; \
- lapack_int lda = matrix.outerStride(); \
- lapack_int matrix_order = MKLCOLROW; \
+ lapack_int n = internal::convert_index<lapack_int>(matrix.cols()), sdim, info; \
+ lapack_int matrix_order = LAPACKE_COLROW; \
char jobvs, sort='N'; \
- LAPACK_##MKLPREFIX_U##_SELECT2 select = 0; \
+ LAPACK_##LAPACKE_PREFIX_U##_SELECT2 select = 0; \
jobvs = (computeU) ? 'V' : 'N'; \
m_matU.resize(n, n); \
- lapack_int ldvs = m_matU.outerStride(); \
+ lapack_int ldvs = internal::convert_index<lapack_int>(m_matU.outerStride()); \
m_matT = matrix; \
+ lapack_int lda = internal::convert_index<lapack_int>(m_matT.outerStride()); \
Matrix<EIGTYPE, Dynamic, Dynamic> wr, wi; \
wr.resize(n, 1); wi.resize(n, 1); \
- info = LAPACKE_##MKLPREFIX##gees( matrix_order, jobvs, sort, select, n, (MKLTYPE*)m_matT.data(), lda, &sdim, (MKLTYPE*)wr.data(), (MKLTYPE*)wi.data(), (MKLTYPE*)m_matU.data(), ldvs ); \
+ info = LAPACKE_##LAPACKE_PREFIX##gees( matrix_order, jobvs, sort, select, n, (LAPACKE_TYPE*)m_matT.data(), lda, &sdim, (LAPACKE_TYPE*)wr.data(), (LAPACKE_TYPE*)wi.data(), (LAPACKE_TYPE*)m_matU.data(), ldvs ); \
if(info == 0) \
m_info = Success; \
else \
@@ -73,11 +67,11 @@
\
}
-EIGEN_MKL_SCHUR_REAL(double, double, d, D, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_SCHUR_REAL(float, float, s, S, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_SCHUR_REAL(double, double, d, D, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_SCHUR_REAL(float, float, s, S, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_SCHUR_REAL(double, double, d, D, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_SCHUR_REAL(float, float, s, S, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_SCHUR_REAL(double, double, d, D, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_SCHUR_REAL(float, float, s, S, RowMajor, LAPACK_ROW_MAJOR)
} // end namespace Eigen
-#endif // EIGEN_REAL_SCHUR_MKL_H
+#endif // EIGEN_REAL_SCHUR_LAPACKE_H
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index c2e76c8..9ddd553 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -20,6 +20,8 @@
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,9 @@
/** \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;
/** \brief Real scalar type for \p _MatrixType.
*
@@ -98,6 +102,7 @@
*/
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.
*
@@ -109,6 +114,7 @@
* Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp
* Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out
*/
+ EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver()
: m_eivec(),
m_eivalues(),
@@ -128,7 +134,8 @@
*
* \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),
@@ -150,13 +157,15 @@
*
* \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.
@@ -189,24 +198,45 @@
*
* \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.
@@ -225,7 +255,8 @@
*
* \sa eigenvalues()
*/
- const MatrixType& eigenvectors() const
+ EIGEN_DEVICE_FUNC
+ const EigenvectorsType& eigenvectors() const
{
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
@@ -247,6 +278,7 @@
*
* \sa eigenvectors(), MatrixBase::eigenvalues()
*/
+ EIGEN_DEVICE_FUNC
const RealVectorType& eigenvalues() const
{
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
@@ -268,9 +300,9 @@
* 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.");
@@ -293,9 +325,9 @@
* 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.");
@@ -307,6 +339,7 @@
*
* \returns \c Success if computation was succesful, \c NoConvergence otherwise.
*/
+ EIGEN_DEVICE_FUNC
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
@@ -320,43 +353,13 @@
*/
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()
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
}
- MatrixType m_eivec;
+ EigenvectorsType m_eivec;
RealVectorType m_eivalues;
typename TridiagonalizationType::SubDiagonalType m_subdiag;
ComputationInfo m_info;
@@ -364,6 +367,7 @@
bool m_eigenvectorsOk;
};
+namespace internal {
/** \internal
*
* \eigenvalues_module \ingroup Eigenvalues_Module
@@ -371,8 +375,12 @@
* 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.
@@ -380,17 +388,21 @@
* Implemented from Golub's "Matrix Computations", algorithm 8.3.2:
* "implicit symmetric QR step with Wilkinson shift"
*/
-namespace internal {
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
@@ -402,7 +414,8 @@
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;
@@ -413,7 +426,7 @@
// declare some aliases
RealVectorType& diag = m_eivalues;
- MatrixType& mat = m_eivec;
+ EigenvectorsType& mat = m_eivec;
// map the matrix coefficients to [-1:1] to avoid over- and underflow.
mat = matrix.template triangularView<Lower>();
@@ -422,58 +435,8 @@
mat.template triangularView<Lower>() /= scale;
m_subdiag.resize(n-1);
internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
-
- Index end = n-1;
- Index start = 0;
- Index iter = 0; // total number of iterations
- 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;
-
- // find the largest unreduced block
- while (end>0 && m_subdiag[end-1]==0)
- {
- end--;
- }
- if (end<=0)
- break;
-
- // if we spent too many iterations, we give up
- iter++;
- if(iter > m_maxIterations * n) break;
-
- start = end - 1;
- while (start>0 && m_subdiag[start-1]!=0)
- start--;
-
- internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
- }
-
- if (iter <= m_maxIterations * n)
- m_info = Success;
- else
- m_info = NoConvergence;
-
- // Sort eigenvalues and corresponding vectors.
- // TODO make the sort optional ?
- // TODO use a better sort algorithm !!
- if (m_info == Success)
- {
- for (Index i = 0; i < n-1; ++i)
- {
- Index k;
- m_eivalues.segment(i,n-i).minCoeff(&k);
- if (k > 0)
- {
- std::swap(m_eivalues[i], m_eivalues[k+i]);
- if(computeEigenvectors)
- m_eivec.col(i).swap(m_eivec.col(k+i));
- }
- }
- }
+ m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
// scale back the eigen values
m_eivalues *= scale;
@@ -483,11 +446,107 @@
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(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 && subdiag[end-1]==RealScalar(0))
+ {
+ end--;
+ }
+ if (end<=0)
+ break;
+
+ // if we spent too many iterations, we give up
+ iter++;
+ if(iter > maxIterations * n) break;
+
+ start = end - 1;
+ while (start>0 && subdiag[start-1]!=0)
+ start--;
+
+ internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
+ }
+ if (iter <= maxIterations * n)
+ info = Success;
+ else
+ info = NoConvergence;
+
+ // Sort eigenvalues and corresponding vectors.
+ // TODO make the sort optional ?
+ // TODO use a better sort algorithm !!
+ if (info == Success)
+ {
+ for (Index i = 0; i < n-1; ++i)
+ {
+ Index k;
+ diag.segment(i,n-i).minCoeff(&k);
+ if (k > 0)
+ {
+ std::swap(diag[i], diag[k+i]);
+ if(computeEigenvectors)
+ eivec.col(i).swap(eivec.col(k+i));
+ }
+ }
+ }
+ return info;
+}
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); }
};
@@ -497,20 +556,22 @@
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
@@ -523,14 +584,12 @@
// 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);
@@ -543,6 +602,7 @@
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;
@@ -562,6 +622,7 @@
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());
@@ -570,7 +631,7 @@
&& "invalid option parameter");
bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
- MatrixType& eivecs = solver.m_eivec;
+ EigenvectorsType& eivecs = solver.m_eivec;
VectorType& eivals = solver.m_eivalues;
// Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
@@ -603,7 +664,7 @@
Index k(0), l(2);
if(d0 > d1)
{
- std::swap(k,l);
+ numext::swap(k,l);
d0 = d1;
}
@@ -647,12 +708,15 @@
};
// 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;
@@ -662,28 +726,33 @@
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
&& "invalid option parameter");
bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
- MatrixType& eivecs = solver.m_eivec;
+ 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)
{
@@ -711,10 +780,11 @@
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;
@@ -724,6 +794,7 @@
}
template<typename MatrixType>
+EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
::computeDirect(const MatrixType& matrix, int options)
{
@@ -733,6 +804,7 @@
namespace internal {
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;
@@ -744,14 +816,14 @@
// 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;
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h
similarity index 65%
rename from Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h
rename to Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h
index 17c0dad..b0c947d 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h
@@ -25,38 +25,36 @@
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
********************************************************************************
- * Content : Eigen bindings to Intel(R) MKL
+ * Content : Eigen bindings to LAPACKe
* Self-adjoint eigenvalues/eigenvectors.
********************************************************************************
*/
-#ifndef EIGEN_SAEIGENSOLVER_MKL_H
-#define EIGEN_SAEIGENSOLVER_MKL_H
-
-#include "Eigen/src/Core/util/MKL_support.h"
+#ifndef EIGEN_SAEIGENSOLVER_LAPACKE_H
+#define EIGEN_SAEIGENSOLVER_LAPACKE_H
namespace Eigen {
-/** \internal Specialization for the data types supported by MKL */
+/** \internal Specialization for the data types supported by LAPACKe */
-#define EIGEN_MKL_EIG_SELFADJ(EIGTYPE, MKLTYPE, MKLRTYPE, MKLNAME, EIGCOLROW, MKLCOLROW ) \
-template<> inline \
+#define EIGEN_LAPACKE_EIG_SELFADJ_2(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME, EIGCOLROW ) \
+template<> template<typename InputType> inline \
SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \
-SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW>& matrix, int options) \
+SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBase<InputType>& matrix, int options) \
{ \
eigen_assert(matrix.cols() == matrix.rows()); \
eigen_assert((options&~(EigVecMask|GenEigMask))==0 \
&& (options&EigVecMask)!=EigVecMask \
&& "invalid option parameter"); \
bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; \
- lapack_int n = matrix.cols(), lda, matrix_order, info; \
+ lapack_int n = internal::convert_index<lapack_int>(matrix.cols()), lda, info; \
m_eivalues.resize(n,1); \
m_subdiag.resize(n-1); \
m_eivec = matrix; \
\
if(n==1) \
{ \
- m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0)); \
+ m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0)); \
if(computeEigenvectors) m_eivec.setOnes(n,n); \
m_info = Success; \
m_isInitialized = true; \
@@ -64,28 +62,25 @@
return *this; \
} \
\
- lda = matrix.outerStride(); \
- matrix_order=MKLCOLROW; \
+ lda = internal::convert_index<lapack_int>(m_eivec.outerStride()); \
char jobz, uplo='L'/*, range='A'*/; \
jobz = computeEigenvectors ? 'V' : 'N'; \
\
- info = LAPACKE_##MKLNAME( matrix_order, jobz, uplo, n, (MKLTYPE*)m_eivec.data(), lda, (MKLRTYPE*)m_eivalues.data() ); \
+ info = LAPACKE_##LAPACKE_NAME( LAPACK_COL_MAJOR, jobz, uplo, n, (LAPACKE_TYPE*)m_eivec.data(), lda, (LAPACKE_RTYPE*)m_eivalues.data() ); \
m_info = (info==0) ? Success : NoConvergence; \
m_isInitialized = true; \
m_eigenvectorsOk = computeEigenvectors; \
return *this; \
}
+#define EIGEN_LAPACKE_EIG_SELFADJ(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME ) \
+ EIGEN_LAPACKE_EIG_SELFADJ_2(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME, ColMajor ) \
+ EIGEN_LAPACKE_EIG_SELFADJ_2(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME, RowMajor )
-EIGEN_MKL_EIG_SELFADJ(double, double, double, dsyev, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_EIG_SELFADJ(float, float, float, ssyev, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_EIG_SELFADJ(dcomplex, MKL_Complex16, double, zheev, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_EIG_SELFADJ(scomplex, MKL_Complex8, float, cheev, ColMajor, LAPACK_COL_MAJOR)
-
-EIGEN_MKL_EIG_SELFADJ(double, double, double, dsyev, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_EIG_SELFADJ(float, float, float, ssyev, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_EIG_SELFADJ(dcomplex, MKL_Complex16, double, zheev, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_EIG_SELFADJ(scomplex, MKL_Complex8, float, cheev, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_EIG_SELFADJ(double, double, double, dsyev)
+EIGEN_LAPACKE_EIG_SELFADJ(float, float, float, ssyev)
+EIGEN_LAPACKE_EIG_SELFADJ(dcomplex, lapack_complex_double, double, zheev)
+EIGEN_LAPACKE_EIG_SELFADJ(scomplex, lapack_complex_float, float, cheev)
} // end namespace Eigen
diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h
index 192278d..1d102c1 100644
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h
@@ -18,8 +18,10 @@
template<typename MatrixType> struct TridiagonalizationMatrixTReturnType;
template<typename MatrixType>
struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
+ : public traits<typename MatrixType::PlainObject>
{
- typedef typename MatrixType::PlainObject ReturnType;
+ typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix?
+ enum { Flags = 0 };
};
template<typename MatrixType, typename CoeffVectorType>
@@ -67,7 +69,7 @@
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef typename MatrixType::Index Index;
+ typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
enum {
Size = MatrixType::RowsAtCompileTime,
@@ -89,10 +91,8 @@
>::type DiagonalReturnType;
typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
- typename internal::add_const_on_value_type<typename Diagonal<
- Block<const MatrixType,SizeMinusOne,SizeMinusOne> >::RealReturnType>::type,
- const Diagonal<
- Block<const MatrixType,SizeMinusOne,SizeMinusOne> >
+ typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type,
+ const Diagonal<const MatrixType, -1>
>::type SubDiagonalReturnType;
/** \brief Return type of matrixQ() */
@@ -110,7 +110,7 @@
*
* \sa compute() for an example.
*/
- Tridiagonalization(Index size = Size==Dynamic ? 2 : Size)
+ explicit Tridiagonalization(Index size = Size==Dynamic ? 2 : Size)
: m_matrix(size,size),
m_hCoeffs(size > 1 ? size-1 : 1),
m_isInitialized(false)
@@ -126,8 +126,9 @@
* Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp
* Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
*/
- Tridiagonalization(const MatrixType& matrix)
- : m_matrix(matrix),
+ template<typename InputType>
+ explicit Tridiagonalization(const EigenBase<InputType>& matrix)
+ : m_matrix(matrix.derived()),
m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
m_isInitialized(false)
{
@@ -152,9 +153,10 @@
* Example: \include Tridiagonalization_compute.cpp
* Output: \verbinclude Tridiagonalization_compute.out
*/
- Tridiagonalization& compute(const MatrixType& matrix)
+ template<typename InputType>
+ Tridiagonalization& compute(const EigenBase<InputType>& matrix)
{
- m_matrix = matrix;
+ m_matrix = matrix.derived();
m_hCoeffs.resize(matrix.rows()-1, 1);
internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
m_isInitialized = true;
@@ -305,7 +307,7 @@
Tridiagonalization<MatrixType>::diagonal() const
{
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
- return m_matrix.diagonal();
+ return m_matrix.diagonal().real();
}
template<typename MatrixType>
@@ -313,8 +315,7 @@
Tridiagonalization<MatrixType>::subDiagonal() const
{
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
- Index n = m_matrix.rows();
- return Block<const MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
+ return m_matrix.template diagonal<-1>().real();
}
namespace internal {
@@ -346,7 +347,6 @@
void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
{
using numext::conj;
- typedef typename MatrixType::Index Index;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
Index n = matA.rows();
@@ -367,10 +367,10 @@
hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
* (conj(h) * matA.col(i).tail(remainingSize)));
- 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);
+ hCoeffs.tail(n-i-1) += (conj(h)*RealScalar(-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);
+ .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1));
matA.col(i).coeffRef(i+1) = beta;
hCoeffs.coeffRef(i) = h;
@@ -438,7 +438,6 @@
{
typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
- typedef typename MatrixType::Index Index;
template<typename DiagonalType, typename SubDiagonalType>
static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
{
@@ -467,9 +466,10 @@
static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
{
using std::sqrt;
+ const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
diag[0] = mat(0,0);
RealScalar v1norm2 = numext::abs2(mat(2,0));
- if(v1norm2 == RealScalar(0))
+ if(v1norm2 <= tol)
{
diag[1] = mat(1,1);
diag[2] = mat(2,2);
@@ -526,7 +526,6 @@
template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
: public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
{
- typedef typename MatrixType::Index Index;
public:
/** \brief Constructor.
*