Tempest_in_Action/tempest-devel/resources/3rdparty/patches/eigen.patch

From 4d2074e54b27877bcd694add91447b569d88e918 Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Wed, 20 May 2020 15:48:43 +0200
Subject: [PATCH] Eigen Patch

---
 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       | 79 +++++++++++++++-------
 5 files changed, 63 insertions(+), 36 deletions(-)

diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index 7104831..9aeceb0 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 b57f068..be15e5e 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 c32d8d8..03c8aaf 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 822cf32..5dc819c 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 a86dac9..f1ebd3c 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, 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.
  * 
@@ -70,42 +109,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];
   }
@@ -125,7 +152,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.21.0 (Apple Git-122)