aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen')
-rw-r--r--Eigen/Cholesky6
-rw-r--r--Eigen/src/Cholesky/LDLT.h6
-rw-r--r--Eigen/src/Cholesky/LLT.h6
-rw-r--r--Eigen/src/Core/MathFunctions.h9
-rw-r--r--Eigen/src/Core/arch/NEON/PacketMath.h1
-rw-r--r--Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h2
-rw-r--r--Eigen/src/Core/util/MKL_support.h19
-rw-r--r--Eigen/src/Core/util/Memory.h2
-rw-r--r--Eigen/src/Eigenvalues/EigenSolver.h40
-rw-r--r--Eigen/src/SVD/JacobiSVD.h9
10 files changed, 85 insertions, 15 deletions
diff --git a/Eigen/Cholesky b/Eigen/Cholesky
index f727f5d89..7314d326c 100644
--- a/Eigen/Cholesky
+++ b/Eigen/Cholesky
@@ -10,9 +10,11 @@
*
*
* This module provides two variants of the Cholesky decomposition for selfadjoint (hermitian) matrices.
- * Those decompositions are accessible via the following MatrixBase methods:
- * - MatrixBase::llt(),
+ * Those decompositions are also accessible via the following methods:
+ * - MatrixBase::llt()
* - MatrixBase::ldlt()
+ * - SelfAdjointView::llt()
+ * - SelfAdjointView::ldlt()
*
* \code
* #include <Eigen/Cholesky>
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index b43e85e7f..efac7fe40 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -43,7 +43,7 @@ namespace internal {
* Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky
* decomposition to determine whether a system of equations has a solution.
*
- * \sa MatrixBase::ldlt(), class LLT
+ * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt(), class LLT
*/
template<typename _MatrixType, int _UpLo> class LDLT
{
@@ -179,7 +179,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
* least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function
* computes the least-square solution of \f$ A x = b \f$ is \f$ A \f$ is singular.
*
- * \sa MatrixBase::ldlt()
+ * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt()
*/
template<typename Rhs>
inline const internal::solve_retval<LDLT, Rhs>
@@ -582,6 +582,7 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
#ifndef __CUDACC__
/** \cholesky_module
* \returns the Cholesky decomposition with full pivoting without square root of \c *this
+ * \sa MatrixBase::ldlt()
*/
template<typename MatrixType, unsigned int UpLo>
inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
@@ -592,6 +593,7 @@ SelfAdjointView<MatrixType, UpLo>::ldlt() const
/** \cholesky_module
* \returns the Cholesky decomposition with full pivoting without square root of \c *this
+ * \sa SelfAdjointView::ldlt()
*/
template<typename Derived>
inline const LDLT<typename MatrixBase<Derived>::PlainObject>
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 2201c641e..45ed8438f 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -41,7 +41,7 @@ template<typename MatrixType, int UpLo> struct LLT_Traits;
* Example: \include LLT_example.cpp
* Output: \verbinclude LLT_example.out
*
- * \sa MatrixBase::llt(), class LDLT
+ * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT
*/
/* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
* Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
@@ -115,7 +115,7 @@ template<typename _MatrixType, int _UpLo> class LLT
* Example: \include LLT_solve.cpp
* Output: \verbinclude LLT_solve.out
*
- * \sa solveInPlace(), MatrixBase::llt()
+ * \sa solveInPlace(), MatrixBase::llt(), SelfAdjointView::llt()
*/
template<typename Rhs>
inline const internal::solve_retval<LLT, Rhs>
@@ -468,6 +468,7 @@ MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
#ifndef __CUDACC__
/** \cholesky_module
* \returns the LLT decomposition of \c *this
+ * \sa SelfAdjointView::llt()
*/
template<typename Derived>
inline const LLT<typename MatrixBase<Derived>::PlainObject>
@@ -478,6 +479,7 @@ MatrixBase<Derived>::llt() const
/** \cholesky_module
* \returns the LLT decomposition of \c *this
+ * \sa SelfAdjointView::llt()
*/
template<typename MatrixType, unsigned int UpLo>
inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index 63fb92b75..20fc2be74 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -669,6 +669,15 @@ bool (isfinite)(const T& x)
return x<NumTraits<T>::highest() && x>NumTraits<T>::lowest();
}
+template<typename T>
+EIGEN_DEVICE_FUNC
+bool (isfinite)(const std::complex<T>& x)
+{
+ using std::real;
+ using std::imag;
+ return isfinite(real(x)) && isfinite(imag(x));
+}
+
} // end namespace numext
namespace internal {
diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h
index 7f3301f51..c9b9e5e9b 100644
--- a/Eigen/src/Core/arch/NEON/PacketMath.h
+++ b/Eigen/src/Core/arch/NEON/PacketMath.h
@@ -97,6 +97,7 @@ template<> struct packet_traits<int> : default_packet_traits
// workaround gcc 4.2, 4.3 and 4.4 compilatin issue
EIGEN_STRONG_INLINE float32x4_t vld1q_f32(const float* x) { return ::vld1q_f32((const float32_t*)x); }
EIGEN_STRONG_INLINE float32x2_t vld1_f32 (const float* x) { return ::vld1_f32 ((const float32_t*)x); }
+EIGEN_STRONG_INLINE float32x2_t vld1_dup_f32 (const float* x) { return ::vld1_dup_f32 ((const float32_t*)x); }
EIGEN_STRONG_INLINE void vst1q_f32(float* to, float32x4_t from) { ::vst1q_f32((float32_t*)to,from); }
EIGEN_STRONG_INLINE void vst1_f32 (float* to, float32x2_t from) { ::vst1_f32 ((float32_t*)to,from); }
#endif
diff --git a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h
index ffa871cae..225b994d1 100644
--- a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h
+++ b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h
@@ -265,6 +265,8 @@ template<typename MatrixType, unsigned int UpLo>
template<typename ProductDerived, typename _Lhs, typename _Rhs>
TriangularView<MatrixType,UpLo>& TriangularView<MatrixType,UpLo>::assignProduct(const ProductBase<ProductDerived, _Lhs,_Rhs>& prod, const Scalar& alpha)
{
+ eigen_assert(m_matrix.rows() == prod.rows() && m_matrix.cols() == prod.cols());
+
general_product_to_triangular_selector<MatrixType, ProductDerived, UpLo, (_Lhs::ColsAtCompileTime==1) || (_Rhs::RowsAtCompileTime==1)>::run(m_matrix.const_cast_derived(), prod.derived(), alpha);
return *this;
diff --git a/Eigen/src/Core/util/MKL_support.h b/Eigen/src/Core/util/MKL_support.h
index 1e6e355d6..8acca9c8c 100644
--- a/Eigen/src/Core/util/MKL_support.h
+++ b/Eigen/src/Core/util/MKL_support.h
@@ -54,8 +54,25 @@
#endif
#if defined EIGEN_USE_MKL
+# include <mkl.h>
+/*Check IMKL version for compatibility: < 10.3 is not usable with Eigen*/
+# ifndef INTEL_MKL_VERSION
+# undef EIGEN_USE_MKL /* INTEL_MKL_VERSION is not even defined on older versions */
+# elif INTEL_MKL_VERSION < 100305 /* the intel-mkl-103-release-notes say this was when the lapacke.h interface was added*/
+# undef EIGEN_USE_MKL
+# endif
+# ifndef EIGEN_USE_MKL
+ /*If the MKL version is too old, undef everything*/
+# undef EIGEN_USE_MKL_ALL
+# undef EIGEN_USE_BLAS
+# undef EIGEN_USE_LAPACKE
+# undef EIGEN_USE_MKL_VML
+# undef EIGEN_USE_LAPACKE_STRICT
+# undef EIGEN_USE_LAPACKE
+# endif
+#endif
-#include <mkl.h>
+#if defined EIGEN_USE_MKL
#include <mkl_lapacke.h>
#define EIGEN_MKL_VML_THRESHOLD 128
diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h
index 34b387ded..e1a12aef1 100644
--- a/Eigen/src/Core/util/Memory.h
+++ b/Eigen/src/Core/util/Memory.h
@@ -89,7 +89,7 @@ inline void throw_std_bad_alloc()
#ifdef EIGEN_EXCEPTIONS
throw std::bad_alloc();
#else
- std::size_t huge = -1;
+ std::size_t huge = static_cast<std::size_t>(-1);
new int[huge];
#endif
}
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index bf20e03ef..d2563d470 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -275,10 +275,11 @@ template<typename _MatrixType> class EigenSolver
*/
EigenSolver& compute(const MatrixType& 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. */
@@ -302,6 +303,7 @@ template<typename _MatrixType> class EigenSolver
EigenvalueType m_eivalues;
bool m_isInitialized;
bool m_eigenvectorsOk;
+ ComputationInfo m_info;
RealSchur<MatrixType> m_realSchur;
MatrixType m_matT;
@@ -366,12 +368,16 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect
{
using std::sqrt;
using std::abs;
+ using std::max;
+ using numext::isfinite;
eigen_assert(matrix.cols() == matrix.rows());
// Reduce to real Schur form.
m_realSchur.compute(matrix, computeEigenvectors);
+
+ m_info = m_realSchur.info();
- if (m_realSchur.info() == Success)
+ if (m_info == Success)
{
m_matT = m_realSchur.matrixT();
if (computeEigenvectors)
@@ -385,14 +391,40 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect
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 = (max)(abs(p),(max)(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;
}
}
@@ -581,7 +613,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
}
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/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index eee31ca97..439eb5d29 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -415,6 +415,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
JacobiRotation<RealScalar> *j_right)
{
using std::sqrt;
+ using std::abs;
Matrix<RealScalar,2,2> m;
m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
@@ -428,9 +429,11 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
}
else
{
- RealScalar u = d / t;
- rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u));
- rot1.s() = rot1.c() * u;
+ RealScalar t2d2 = numext::hypot(t,d);
+ rot1.c() = abs(t)/t2d2;
+ rot1.s() = d/t2d2;
+ if(t<RealScalar(0))
+ rot1.s() = -rot1.s();
}
m.applyOnTheLeft(0,1,rot1);
j_right->makeJacobi(m,0,1);
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-04 19:21:44 +0000