From 739ef0b8b4fd2a01188e51d8eb6bb42558b1b0d5 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Tue, 19 May 2020 10:35:31 +0200 Subject: [PATCH] Updated Eigen version. This fixes issue #77 We also now clone Eigen from the eigen git and perform a patch step. This makes updating the eigen versions more easy in the future. --- resources/3rdparty/CMakeLists.txt | 20 +- resources/3rdparty/patches/eigen.patch | 271 +++++++++++++++++++++++++ 2 files changed, 289 insertions(+), 2 deletions(-) create mode 100644 resources/3rdparty/patches/eigen.patch diff --git a/resources/3rdparty/CMakeLists.txt b/resources/3rdparty/CMakeLists.txt index 2a5f13dab..ecc01212e 100644 --- a/resources/3rdparty/CMakeLists.txt +++ b/resources/3rdparty/CMakeLists.txt @@ -44,8 +44,24 @@ list(APPEND STORM_DEP_TARGETS gmm) ## ############################################################# -add_imported_library_interface(StormEigen33 "${PROJECT_SOURCE_DIR}/resources/3rdparty/eigen-3.3-beta1") -list(APPEND STORM_DEP_TARGETS StormEigen33) +# Checkout Eigen version 3.3.7 +ExternalProject_Add( + eigen_src + GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git + GIT_SHALLOW 1 + GIT_TAG 21ae2afd4edaa1b69782c67a54182d34efe43f9c + SOURCE_DIR ${STORM_3RDPARTY_SOURCE_DIR}/eigen + PATCH_COMMAND git apply ${STORM_3RDPARTY_SOURCE_DIR}/patches/eigen.patch + UPDATE_COMMAND "" + CONFIGURE_COMMAND "" + BUILD_COMMAND "" + INSTALL_COMMAND "" + LOG_INSTALL ON +) +add_imported_library_interface(Eigen "${STORM_3RDPARTY_SOURCE_DIR}/eigen") +list(APPEND STORM_DEP_TARGETS Eigen) +add_dependencies(Eigen eigen_src) + ############################################################# diff --git a/resources/3rdparty/patches/eigen.patch b/resources/3rdparty/patches/eigen.patch new file mode 100644 index 000000000..62f8975b6 --- /dev/null +++ b/resources/3rdparty/patches/eigen.patch @@ -0,0 +1,271 @@ +From 9f1a66bef29005a4cbe178241cad563fcb0dfa5f Mon Sep 17 00:00:00 2001 +From: Tim Quatmann +Date: Mon, 18 May 2020 11:58:47 +0200 +Subject: [PATCH 1/2] Adaptions for SparseLU to work with Scalar types that do + not default construct from a double (like CLN numbers) or that do not have an + operator< or std::abs + +--- + Eigen/src/SparseLU/SparseLU.h | 10 +-- + Eigen/src/SparseLU/SparseLU_column_bmod.h | 2 +- + Eigen/src/SparseLU/SparseLU_copy_to_ucol.h | 2 +- + Eigen/src/SparseLU/SparseLU_panel_bmod.h | 6 +- + Eigen/src/SparseLU/SparseLU_pivotL.h | 74 +++++++++++++++------- + 5 files changed, 62 insertions(+), 32 deletions(-) + +diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h +index 7104831c0..9aeceb021 100644 +--- a/Eigen/src/SparseLU/SparseLU.h ++++ b/Eigen/src/SparseLU/SparseLU.h +@@ -97,12 +97,12 @@ class SparseLU : public SparseSolverBase >, + }; + + public: +- SparseLU():m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) ++ SparseLU():m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1),m_detPermR(1) + { + initperfvalues(); + } + explicit SparseLU(const MatrixType& matrix) +- : m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) ++ : m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1),m_detPermR(1) + { + initperfvalues(); + compute(matrix); +@@ -255,7 +255,7 @@ class SparseLU : public SparseSolverBase >, + using std::abs; + eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); + // Initialize with the determinant of the row matrix +- Scalar det = Scalar(1.); ++ Scalar det = Scalar(1); + // Note that the diagonal blocks of U are stored in supernodes, + // which are available in the L part :) + for (Index j = 0; j < this->cols(); ++j) +@@ -286,7 +286,7 @@ class SparseLU : public SparseSolverBase >, + using std::abs; + + eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); +- Scalar det = Scalar(0.); ++ Scalar det = Scalar(0); + for (Index j = 0; j < this->cols(); ++j) + { + for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) +@@ -338,7 +338,7 @@ class SparseLU : public SparseSolverBase >, + { + eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); + // Initialize with the determinant of the row matrix +- Scalar det = Scalar(1.); ++ Scalar det = Scalar(1); + // Note that the diagonal blocks of U are stored in supernodes, + // which are available in the L part :) + for (Index j = 0; j < this->cols(); ++j) +diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h +index b57f06802..be15e5ee6 100644 +--- a/Eigen/src/SparseLU/SparseLU_column_bmod.h ++++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h +@@ -129,7 +129,7 @@ Index SparseLUImpl::column_bmod(const Index jcol, const Ind + { + irow = glu.lsub(isub); + glu.lusup(nextlu) = dense(irow); +- dense(irow) = Scalar(0.0); ++ dense(irow) = Scalar(0); + ++nextlu; + } + +diff --git a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h +index c32d8d8b1..03c8aaf85 100644 +--- a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h ++++ b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h +@@ -87,7 +87,7 @@ Index SparseLUImpl::copy_to_ucol(const Index jcol, const In + irow = glu.lsub(isub); + glu.usub(nextu) = perm_r(irow); // Unlike the L part, the U part is stored in its final order + glu.ucol(nextu) = dense(irow); +- dense(irow) = Scalar(0.0); ++ dense(irow) = Scalar(0); + nextu++; + isub++; + } +diff --git a/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/Eigen/src/SparseLU/SparseLU_panel_bmod.h +index 822cf32c3..5dc819cf7 100644 +--- a/Eigen/src/SparseLU/SparseLU_panel_bmod.h ++++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h +@@ -122,7 +122,7 @@ void SparseLUImpl::panel_bmod(const Index m, const Index w, + + Index isub = lptr + no_zeros; + Index off = u_rows-segsize; +- for (Index i = 0; i < off; i++) U(i,u_col) = 0; ++ for (Index i = 0; i < off; i++) U(i,u_col) = Scalar(0); + for (Index i = 0; i < segsize; i++) + { + Index irow = glu.lsub(isub); +@@ -172,7 +172,7 @@ void SparseLUImpl::panel_bmod(const Index m, const Index w, + { + Index irow = glu.lsub(isub++); + dense_col(irow) = U.coeff(i+off,u_col); +- U.coeffRef(i+off,u_col) = 0; ++ U.coeffRef(i+off,u_col) = Scalar(0); + } + + // Scatter l into SPA dense[] +@@ -180,7 +180,7 @@ void SparseLUImpl::panel_bmod(const Index m, const Index w, + { + Index irow = glu.lsub(isub++); + dense_col(irow) -= L.coeff(i,u_col); +- L.coeffRef(i,u_col) = 0; ++ L.coeffRef(i,u_col) = Scalar(0); + } + u_col++; + } +diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h +index a86dac93f..2ad56531c 100644 +--- a/Eigen/src/SparseLU/SparseLU_pivotL.h ++++ b/Eigen/src/SparseLU/SparseLU_pivotL.h +@@ -32,7 +32,46 @@ + + namespace Eigen { + namespace internal { +- ++ template::value, int>::type = 0> ++ void findLargestAbsolutePivotElement(const Index& nsupc, const Index& nsupr, Scalar const* lu_col_ptr, StorageIndex const* lsub_ptr, Index const& diagind, Index& diag, Index& pivptr, Scalar& pivmax) { ++ double rtemp = 0; ++ for (Index isub = nsupc; isub < nsupr; ++isub) { ++ rtemp = std::abs(lu_col_ptr[isub]); ++ if (rtemp > pivmax) { ++ pivmax = rtemp; ++ pivptr = isub; ++ } ++ if (lsub_ptr[isub] == diagind) diag = isub; ++ } ++ } ++ ++ template::value, int>::type = 0> ++ void findLargestAbsolutePivotElement(const Index& nsupc, const Index& nsupr, Scalar const* lu_col_ptr, StorageIndex const* lsub_ptr, Index const& diagind, Index& diag, Index& pivptr, Scalar& pivmax) { ++ bool foundPivotElement = false; ++ for (Index isub = nsupc; isub < nsupr; ++isub) { ++ Scalar const& rtemp = lu_col_ptr[isub]; ++ if (!foundPivotElement && !storm::utility::isZero(rtemp)) { ++ foundPivotElement = true; ++ pivmax = rtemp; ++ pivptr = isub; ++ } ++ if (lsub_ptr[isub] == diagind) diag = isub; ++ } ++ } ++ ++ template::value, int>::type = 0> ++ bool diagonalElementCanBePivot(Scalar const* lu_col_ptr, Index const& diag, Scalar const& diagpivotthresh, Scalar const& pivmax) { ++ Scalar thresh = diagpivotthresh * pivmax; ++ double rtemp = std::abs(lu_col_ptr[diag]); ++ if (!storm::utility::isZero(rtemp) && rtemp >= thresh) return true; ++ return false; ++ } ++ ++ template::value, int>::type = 0> ++ bool diagonalElementCanBePivot(Scalar const* lu_col_ptr, Index const& diag, Scalar const& diagpivotthresh, Scalar const& pivmax) { ++ if (!storm::utility::isZero(lu_col_ptr[diag])) return true; ++ return false; ++ } + /** + * \brief Performs the numerical pivotin on the current column of L, and the CDIV operation. + * +@@ -70,26 +109,19 @@ Index SparseLUImpl::pivotL(const Index jcol, const RealScal + StorageIndex* lsub_ptr = &(glu.lsub.data()[lptr]); // Start of row indices of the supernode + + // Determine the largest abs numerical value for partial pivoting +- Index diagind = iperm_c(jcol); // diagonal index +- RealScalar pivmax(-1.0); +- Index pivptr = nsupc; +- Index diag = emptyIdxLU; +- RealScalar rtemp; +- Index isub, icol, itemp, k; +- for (isub = nsupc; isub < nsupr; ++isub) { +- using std::abs; +- rtemp = abs(lu_col_ptr[isub]); +- if (rtemp > pivmax) { +- pivmax = rtemp; +- pivptr = isub; +- } +- if (lsub_ptr[isub] == diagind) diag = isub; +- } ++ Index diagind = iperm_c(jcol); // diagonal index ++ RealScalar pivmax(-1); ++ Index pivptr = nsupc; ++ Index diag = emptyIdxLU; ++ Index icol, itemp, k; ++ ++ findLargestAbsolutePivotElement(nsupc, nsupr, lu_col_ptr, lsub_ptr, diagind, diag, pivptr, pivmax); + + // Test for singularity +- if ( pivmax <= RealScalar(0.0) ) { ++ bool columnStructurallyEmpty = nsupr <= nsupc; ++ if ( columnStructurallyEmpty || storm::utility::isZero(pivmax)) { + // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero +- pivrow = pivmax < RealScalar(0.0) ? diagind : lsub_ptr[pivptr]; ++ pivrow = columnStructurallyEmpty ? diagind : lsub_ptr[pivptr]; + perm_r(pivrow) = StorageIndex(jcol); + return (jcol+1); + } +@@ -103,9 +135,7 @@ Index SparseLUImpl::pivotL(const Index jcol, const RealScal + if (diag >= 0 ) + { + // Diagonal element exists +- using std::abs; +- rtemp = abs(lu_col_ptr[diag]); +- if (rtemp != RealScalar(0.0) && rtemp >= thresh) pivptr = diag; ++ if (diagonalElementCanBePivot(lu_col_ptr, diag, diagpivotthresh, pivmax)) pivptr = diag; + } + pivrow = lsub_ptr[pivptr]; + } +@@ -125,7 +155,7 @@ Index SparseLUImpl::pivotL(const Index jcol, const RealScal + } + } + // cdiv operations +- Scalar temp = Scalar(1.0) / lu_col_ptr[nsupc]; ++ Scalar temp = Scalar(1) / lu_col_ptr[nsupc]; + for (k = nsupc+1; k < nsupr; k++) + lu_col_ptr[k] *= temp; + return 0; +-- +2.24.2 (Apple Git-127) + + +From 366d73d04bedd7c58fc6b7e96f32a72422ceea89 Mon Sep 17 00:00:00 2001 +From: Tim Quatmann +Date: Tue, 19 May 2020 08:52:28 +0200 +Subject: [PATCH 2/2] Fixed actually using largest pivot with rational numbers + +--- + Eigen/src/SparseLU/SparseLU_pivotL.h | 8 ++++---- + 1 file changed, 4 insertions(+), 4 deletions(-) + +diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h +index 2ad56531c..cc2dd9292 100644 +--- a/Eigen/src/SparseLU/SparseLU_pivotL.h ++++ b/Eigen/src/SparseLU/SparseLU_pivotL.h +@@ -32,11 +32,11 @@ + + namespace Eigen { + namespace internal { +- template::value, int>::type = 0> ++ template::value, int>::type = 0> + void findLargestAbsolutePivotElement(const Index& nsupc, const Index& nsupr, Scalar const* lu_col_ptr, StorageIndex const* lsub_ptr, Index const& diagind, Index& diag, Index& pivptr, Scalar& pivmax) { +- double rtemp = 0; ++ Scalar rtemp = 0; + for (Index isub = nsupc; isub < nsupr; ++isub) { +- rtemp = std::abs(lu_col_ptr[isub]); ++ rtemp = storm::utility::abs(lu_col_ptr[isub]); + if (rtemp > pivmax) { + pivmax = rtemp; + pivptr = isub; +@@ -45,7 +45,7 @@ namespace internal { + } + } + +- template::value, int>::type = 0> ++ template::value, int>::type = 0> + void findLargestAbsolutePivotElement(const Index& nsupc, const Index& nsupr, Scalar const* lu_col_ptr, StorageIndex const* lsub_ptr, Index const& diagind, Index& diag, Index& pivptr, Scalar& pivmax) { + bool foundPivotElement = false; + for (Index isub = nsupc; isub < nsupr; ++isub) { +-- +2.24.2 (Apple Git-127) +