sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

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.
main
Tim Quatmann 6 years ago
parent
dd7dc4b797
commit
739ef0b8b4
2 changed files with 289 additions and 2 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 20
      resources/3rdparty/CMakeLists.txt
  2. 271
      resources/3rdparty/patches/eigen.patch

20
resources/3rdparty/CMakeLists.txt
View File

@ -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)
############################################################# #############################################################

271
resources/3rdparty/patches/eigen.patch
View File

@ -0,0 +1,271 @@
From 9f1a66bef29005a4cbe178241cad563fcb0dfa5f Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
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<SparseLU<_MatrixType,_OrderingType> >,
};
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<SparseLU<_MatrixType,_OrderingType> >,
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<SparseLU<_MatrixType,_OrderingType> >,
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<SparseLU<_MatrixType,_OrderingType> >,
{
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<Scalar,StorageIndex>::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<Scalar,StorageIndex>::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<Scalar,StorageIndex>::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<Scalar,StorageIndex>::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<Scalar,StorageIndex>::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<typename Scalar, typename StorageIndex, typename std::enable_if<std::is_same<Scalar, double>::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<typename Scalar, typename StorageIndex, typename std::enable_if<!std::is_same<Scalar, double>::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<typename Scalar, typename std::enable_if<std::is_same<Scalar, double>::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<typename Scalar, typename std::enable_if<!std::is_same<Scalar, double>::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<Scalar,StorageIndex>::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<Scalar,StorageIndex>::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<Scalar,StorageIndex>::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 <tim.quatmann@cs.rwth-aachen.de>
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<typename Scalar, typename StorageIndex, typename std::enable_if<std::is_same<Scalar, double>::value, int>::type = 0>
+ template<typename Scalar, typename StorageIndex, typename std::enable_if<!std::is_same<Scalar, storm::RationalFunction>::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<typename Scalar, typename StorageIndex, typename std::enable_if<!std::is_same<Scalar, double>::value, int>::type = 0>
+ template<typename Scalar, typename StorageIndex, typename std::enable_if<std::is_same<Scalar, storm::RationalFunction>::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)
Write Preview
Loading…
Cancel
Save