Improving Apple Silicon Build

CHANGELOG.md

@@ -9,6 +9,7 @@ Version 1.7.x
 -------------
 
 ## Version 1.7.1
+- Experimental support for compiling on Apple Silicon
 - Added SoPlex as a possible LP solver
 - Upgraded shipped version of sylvan
 - Upgraded repo / version for carl (for polynomials)

CMakeLists.txt

@@ -191,7 +191,7 @@ set(SHIPPED_CARL_USE_GINAC ON)
 # Detect Apple Silicon and adjust settings
 if(APPLE AND ${CMAKE_SYSTEM_PROCESSOR} MATCHES arm64)
 	message(STATUS "Storm - Detected that target system uses Apple Silicon.")
-	message(WARNING "Storm currently does not compile natively on Apple Silicon! Please consider building Storm using Rosetta 2. For more information visit https://www.stormchecker.org/documentation/obtain-storm/apple-silicon.html")
+	message(WARNING "Compiling natively on Apple Silicon is experimental. Pleare report issues to [email protected]. For more information visit https://www.stormchecker.org/documentation/obtain-storm/apple-silicon.html")
 	set(APPLE_SILICON 1)
 	if(STORM_USE_CLN_EA OR STORM_USE_CLN_RF)
 		message(WARNING "CLN and GiNaC are currently not supported on Apple Silicon-based architectures. Disabling Storm and carl usage of the libraries.")

resources/3rdparty/CMakeLists.txt

@@ -46,22 +46,42 @@ list(APPEND STORM_DEP_TARGETS gmm)
 ##
 #############################################################
 
-# Checkout Eigen version 3.3.9
-message (STATUS "Storm - Including Eigen 3.3.9.")
-ExternalProject_Add(
-        eigen_src
-        GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
-		GIT_SHALLOW 1
-        GIT_TAG 0fd6b4f71dd85b2009ee4d1aeb296e2c11fc9d68
-        SOURCE_DIR ${STORM_3RDPARTY_INCLUDE_DIR}/StormEigen
-        PREFIX ${STORM_3RDPARTY_BINARY_DIR}/StormEigen-3.3.9
-		PATCH_COMMAND git apply ${STORM_3RDPARTY_SOURCE_DIR}/patches/eigen.patch
-        UPDATE_COMMAND ""
-        CONFIGURE_COMMAND ""
-        BUILD_COMMAND ""
-        INSTALL_COMMAND ""
-        LOG_INSTALL ON
-)
+# Checkout Eigen 
+# Version 3.3.9 has a known issue when compiled on Apple Silicon (FullPivLU segfaults, relevant for multi-obj model checking).
+# Version 3.4.0 has issues with gcc https://github.com/moves-rwth/storm/issues/162
+# We use this as a workaround:
+if (APPLE_SILICON)
+	message (STATUS "Storm - Including Eigen 3.4.0 commit b0eded878d5d162d61583a286c0d8a45406ad1bc.")
+	ExternalProject_Add(
+	        eigen_src
+	        GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
+	        GIT_TAG b0eded878d5d162d61583a286c0d8a45406ad1bc
+	        SOURCE_DIR ${STORM_3RDPARTY_INCLUDE_DIR}/StormEigen
+	        PREFIX ${STORM_3RDPARTY_BINARY_DIR}/StormEigen-3.4.0
+			PATCH_COMMAND git apply ${STORM_3RDPARTY_SOURCE_DIR}/patches/eigen340.patch
+	        UPDATE_COMMAND ""
+	        CONFIGURE_COMMAND ""
+	        BUILD_COMMAND ""
+	        INSTALL_COMMAND ""
+	        LOG_INSTALL ON
+	)
+else()
+	message (STATUS "Storm - Including Eigen 3.3.9.")
+	ExternalProject_Add(
+	        eigen_src
+	        GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
+			GIT_SHALLOW 1
+	        GIT_TAG 0fd6b4f71dd85b2009ee4d1aeb296e2c11fc9d68
+	        SOURCE_DIR ${STORM_3RDPARTY_INCLUDE_DIR}/StormEigen
+	        PREFIX ${STORM_3RDPARTY_BINARY_DIR}/StormEigen-3.3.9
+			PATCH_COMMAND git apply ${STORM_3RDPARTY_SOURCE_DIR}/patches/eigen339.patch
+	        UPDATE_COMMAND ""
+	        CONFIGURE_COMMAND ""
+	        BUILD_COMMAND ""
+	        INSTALL_COMMAND ""
+	        LOG_INSTALL ON
+	)
+endif()
 add_imported_library_interface(StormEigen "${STORM_3RDPARTY_INCLUDE_DIR}/StormEigen")
 list(APPEND STORM_DEP_TARGETS StormEigen)
 add_dependencies(StormEigen eigen_src)

resources/3rdparty/patches/eigen340.patch

@@ -0,0 +1,226 @@
+From 938211b9e47c1650fb057e7ff88957cb2d25af7e Mon Sep 17 00:00:00 2001
+From: Tim Quatmann <[email protected]>
+Date: Wed, 12 Apr 2023 21:05:54 +0200
+Subject: [PATCH] stormeigen
+
+---
+ Eigen/src/SparseLU/SparseLU.h              |  8 +--
+ 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       | 79 +++++++++++++++-------
+ 5 files changed, 62 insertions(+), 35 deletions(-)
+
+diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
+index f70aab197..33ba7e1d7 100644
+--- a/Eigen/src/SparseLU/SparseLU.h
++++ b/Eigen/src/SparseLU/SparseLU.h
+@@ -158,12 +158,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);
+@@ -355,7 +355,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)
+@@ -438,7 +438,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 d5c29b353..09abdc402 100644
+--- a/Eigen/src/SparseLU/SparseLU_column_bmod.h
++++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h
+@@ -131,7 +131,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 e06b2a03a..cf91f49df 100644
+--- a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h
++++ b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h
+@@ -89,7 +89,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 fd1243d72..a34cbad6d 100644
+--- a/Eigen/src/SparseLU/SparseLU_panel_bmod.h
++++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h
+@@ -124,7 +124,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); 
+@@ -173,7 +173,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[]
+@@ -181,7 +181,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 6daed9133..78da06f62 100644
+--- a/Eigen/src/SparseLU/SparseLU_pivotL.h
++++ b/Eigen/src/SparseLU/SparseLU_pivotL.h
+@@ -34,7 +34,46 @@
+
+ namespace Eigen {
+ namespace internal {
+-  
++    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) {
++        Scalar rtemp = 0;
++        for (Index isub = nsupc; isub < nsupr; ++isub) {
++            rtemp = storm::utility::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, 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) {
++            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.
+  * 
+@@ -72,42 +111,30 @@ 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);
+   }
+
+-  RealScalar thresh = diagpivotthresh * pivmax; 
+-  
+-  // Choose appropriate pivotal element 
+-  
++  // Choose appropriate pivotal element
+   {
+     // Test if the diagonal element can be used as a pivot (given the threshold value)
+     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];
+   }
+@@ -127,7 +154,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.32.1 (Apple Git-133)
+