aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen/src/Cholesky
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen/src/Cholesky')
-rw-r--r--Eigen/src/Cholesky/CMakeLists.txt6
-rw-r--r--Eigen/src/Cholesky/LDLT.h54
-rw-r--r--Eigen/src/Cholesky/LLT.h20
-rw-r--r--Eigen/src/Cholesky/LLT_LAPACKE.h (renamed from Eigen/src/Cholesky/LLT_MKL.h)35
4 files changed, 79 insertions, 36 deletions
diff --git a/Eigen/src/Cholesky/CMakeLists.txt b/Eigen/src/Cholesky/CMakeLists.txt
deleted file mode 100644
index d01488b41..000000000
--- a/Eigen/src/Cholesky/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_Cholesky_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_Cholesky_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Cholesky COMPONENT Devel
- )
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 538aff956..fcee7b2e3 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -43,6 +43,8 @@ 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.
*
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
* \sa MatrixBase::ldlt(), SelfAdjointView::ldlt(), class LLT
*/
template<typename _MatrixType, int _UpLo> class LDLT
@@ -52,7 +54,6 @@ template<typename _MatrixType, int _UpLo> class LDLT
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
- Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here!
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
UpLo = _UpLo
@@ -61,7 +62,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
typedef typename MatrixType::StorageIndex StorageIndex;
- typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType;
+ typedef Matrix<Scalar, RowsAtCompileTime, 1, 0, MaxRowsAtCompileTime, 1> TmpMatrixType;
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
@@ -97,6 +98,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
/** \brief Constructor with decomposition
*
* This calculates the decomposition for the input \a matrix.
+ *
* \sa LDLT(Index size)
*/
template<typename InputType>
@@ -110,6 +112,23 @@ template<typename _MatrixType, int _UpLo> class LDLT
compute(matrix.derived());
}
+ /** \brief Constructs a LDLT factorization from a given matrix
+ *
+ * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
+ *
+ * \sa LDLT(const EigenBase&)
+ */
+ template<typename InputType>
+ explicit LDLT(EigenBase<InputType>& matrix)
+ : m_matrix(matrix.derived()),
+ m_transpositions(matrix.rows()),
+ m_temporary(matrix.rows()),
+ m_sign(internal::ZeroSign),
+ m_isInitialized(false)
+ {
+ compute(matrix.derived());
+ }
+
/** Clear any existing decomposition
* \sa rankUpdate(w,sigma)
*/
@@ -234,7 +253,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "LDLT is not initialized.");
- return Success;
+ return m_info;
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
@@ -262,6 +281,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
TmpMatrixType m_temporary;
internal::SignMatrix m_sign;
bool m_isInitialized;
+ ComputationInfo m_info;
};
namespace internal {
@@ -279,6 +299,8 @@ template<> struct ldlt_inplace<Lower>
typedef typename TranspositionType::StorageIndex IndexType;
eigen_assert(mat.rows()==mat.cols());
const Index size = mat.rows();
+ bool found_zero_pivot = false;
+ bool ret = true;
if (size <= 1)
{
@@ -337,9 +359,27 @@ template<> struct ldlt_inplace<Lower>
// we should only make sure that we do not introduce INF or NaN values.
// Remark that LAPACK also uses 0 as the cutoff value.
RealScalar realAkk = numext::real(mat.coeffRef(k,k));
- if((rs>0) && (abs(realAkk) > RealScalar(0)))
+ bool pivot_is_valid = (abs(realAkk) > RealScalar(0));
+
+ if(k==0 && !pivot_is_valid)
+ {
+ // The entire diagonal is zero, there is nothing more to do
+ // except filling the transpositions, and checking whether the matrix is zero.
+ sign = ZeroSign;
+ for(Index j = 0; j<size; ++j)
+ {
+ transpositions.coeffRef(j) = IndexType(j);
+ ret = ret && (mat.col(j).tail(size-j-1).array()==Scalar(0)).all();
+ }
+ return ret;
+ }
+
+ if((rs>0) && pivot_is_valid)
A21 /= realAkk;
+ if(found_zero_pivot && pivot_is_valid) ret = false; // factorization failed
+ else if(!pivot_is_valid) found_zero_pivot = true;
+
if (sign == PositiveSemiDef) {
if (realAkk < static_cast<RealScalar>(0)) sign = Indefinite;
} else if (sign == NegativeSemiDef) {
@@ -350,7 +390,7 @@ template<> struct ldlt_inplace<Lower>
}
}
- return true;
+ return ret;
}
// Reference for the algorithm: Davis and Hager, "Multiple Rank
@@ -474,7 +514,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const EigenBase<InputTyp
m_temporary.resize(size);
m_sign = internal::ZeroSign;
- internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign);
+ m_info = internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign) ? Success : NumericalIssue;
m_isInitialized = true;
return *this;
@@ -602,7 +642,6 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
return res;
}
-#ifndef __CUDACC__
/** \cholesky_module
* \returns the Cholesky decomposition with full pivoting without square root of \c *this
* \sa MatrixBase::ldlt()
@@ -624,7 +663,6 @@ MatrixBase<Derived>::ldlt() const
{
return LDLT<PlainObject>(derived());
}
-#endif // __CUDACC__
} // end namespace Eigen
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 19578b216..ddf4875ab 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -41,6 +41,8 @@ template<typename MatrixType, int UpLo> struct LLT_Traits;
* Example: \include LLT_example.cpp
* Output: \verbinclude LLT_example.out
*
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
* \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)
@@ -54,7 +56,6 @@ template<typename _MatrixType, int _UpLo> class LLT
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
- Options = MatrixType::Options,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
@@ -95,6 +96,21 @@ template<typename _MatrixType, int _UpLo> class LLT
compute(matrix.derived());
}
+ /** \brief Constructs a LDLT factorization from a given matrix
+ *
+ * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when
+ * \c MatrixType is a Eigen::Ref.
+ *
+ * \sa LLT(const EigenBase&)
+ */
+ template<typename InputType>
+ explicit LLT(EigenBase<InputType>& matrix)
+ : m_matrix(matrix.derived()),
+ m_isInitialized(false)
+ {
+ compute(matrix.derived());
+ }
+
/** \returns a view of the upper triangular matrix U */
inline typename Traits::MatrixU matrixU() const
{
@@ -491,7 +507,6 @@ MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
return matrixL() * matrixL().adjoint().toDenseMatrix();
}
-#ifndef __CUDACC__
/** \cholesky_module
* \returns the LLT decomposition of \c *this
* \sa SelfAdjointView::llt()
@@ -513,7 +528,6 @@ SelfAdjointView<MatrixType, UpLo>::llt() const
{
return LLT<PlainObject,UpLo>(m_matrix);
}
-#endif // __CUDACC__
} // end namespace Eigen
diff --git a/Eigen/src/Cholesky/LLT_MKL.h b/Eigen/src/Cholesky/LLT_LAPACKE.h
index 0d42cb5bc..bc6489e69 100644
--- a/Eigen/src/Cholesky/LLT_MKL.h
+++ b/Eigen/src/Cholesky/LLT_LAPACKE.h
@@ -25,25 +25,22 @@
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
********************************************************************************
- * Content : Eigen bindings to Intel(R) MKL
+ * Content : Eigen bindings to LAPACKe
* LLt decomposition based on LAPACKE_?potrf function.
********************************************************************************
*/
-#ifndef EIGEN_LLT_MKL_H
-#define EIGEN_LLT_MKL_H
-
-#include "Eigen/src/Core/util/MKL_support.h"
-#include <iostream>
+#ifndef EIGEN_LLT_LAPACKE_H
+#define EIGEN_LLT_LAPACKE_H
namespace Eigen {
namespace internal {
-template<typename Scalar> struct mkl_llt;
+template<typename Scalar> struct lapacke_llt;
-#define EIGEN_MKL_LLT(EIGTYPE, MKLTYPE, MKLPREFIX) \
-template<> struct mkl_llt<EIGTYPE> \
+#define EIGEN_LAPACKE_LLT(EIGTYPE, BLASTYPE, LAPACKE_PREFIX) \
+template<> struct lapacke_llt<EIGTYPE> \
{ \
template<typename MatrixType> \
static inline Index potrf(MatrixType& m, char uplo) \
@@ -53,13 +50,13 @@ template<> struct mkl_llt<EIGTYPE> \
EIGTYPE* a; \
eigen_assert(m.rows()==m.cols()); \
/* Set up parameters for ?potrf */ \
- size = m.rows(); \
+ size = convert_index<lapack_int>(m.rows()); \
StorageOrder = MatrixType::Flags&RowMajorBit?RowMajor:ColMajor; \
matrix_order = StorageOrder==RowMajor ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \
a = &(m.coeffRef(0,0)); \
- lda = m.outerStride(); \
+ lda = convert_index<lapack_int>(m.outerStride()); \
\
- info = LAPACKE_##MKLPREFIX##potrf( matrix_order, uplo, size, (MKLTYPE*)a, lda ); \
+ info = LAPACKE_##LAPACKE_PREFIX##potrf( matrix_order, uplo, size, (BLASTYPE*)a, lda ); \
info = (info==0) ? -1 : info>0 ? info-1 : size; \
return info; \
} \
@@ -69,7 +66,7 @@ template<> struct llt_inplace<EIGTYPE, Lower> \
template<typename MatrixType> \
static Index blocked(MatrixType& m) \
{ \
- return mkl_llt<EIGTYPE>::potrf(m, 'L'); \
+ return lapacke_llt<EIGTYPE>::potrf(m, 'L'); \
} \
template<typename MatrixType, typename VectorType> \
static Index rankUpdate(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) \
@@ -80,7 +77,7 @@ template<> struct llt_inplace<EIGTYPE, Upper> \
template<typename MatrixType> \
static Index blocked(MatrixType& m) \
{ \
- return mkl_llt<EIGTYPE>::potrf(m, 'U'); \
+ return lapacke_llt<EIGTYPE>::potrf(m, 'U'); \
} \
template<typename MatrixType, typename VectorType> \
static Index rankUpdate(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) \
@@ -90,13 +87,13 @@ template<> struct llt_inplace<EIGTYPE, Upper> \
} \
};
-EIGEN_MKL_LLT(double, double, d)
-EIGEN_MKL_LLT(float, float, s)
-EIGEN_MKL_LLT(dcomplex, MKL_Complex16, z)
-EIGEN_MKL_LLT(scomplex, MKL_Complex8, c)
+EIGEN_LAPACKE_LLT(double, double, d)
+EIGEN_LAPACKE_LLT(float, float, s)
+EIGEN_LAPACKE_LLT(dcomplex, lapack_complex_double, z)
+EIGEN_LAPACKE_LLT(scomplex, lapack_complex_float, c)
} // end namespace internal
} // end namespace Eigen
-#endif // EIGEN_LLT_MKL_H
+#endif // EIGEN_LLT_LAPACKE_H
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-18 17:56:08 +0000