From 4d2074e54b27877bcd694add91447b569d88e918 Mon Sep 17 00:00:00 2001 From: Tim Quatmann 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 >, }; 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 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::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::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::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 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::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::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,42 +109,30 @@ 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); } - 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::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)