From 59ba53194d9b7ff87c1c7eb8f7d5b5af284cc352 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 27 Jan 2024 11:52:58 -0500 Subject: [PATCH] template pivoting in the linear algebra now we can do `integrator.linalg_do_pivoting=0` to disable pivoting --- integration/BackwardEuler/be_integrator.H | 17 +++++- integration/VODE/vode_dvjac.H | 8 ++- integration/VODE/vode_dvnlsd.H | 8 ++- integration/_parameters | 3 + util/linpack.H | 67 ++++++++++++++--------- 5 files changed, 73 insertions(+), 30 deletions(-) diff --git a/integration/BackwardEuler/be_integrator.H b/integration/BackwardEuler/be_integrator.H index 9ef6532110..6e972f7e29 100644 --- a/integration/BackwardEuler/be_integrator.H +++ b/integration/BackwardEuler/be_integrator.H @@ -96,14 +96,27 @@ int single_step (BurnT& state, BeT& be, const Real dt) int ierr_linpack; IArray1D pivot; - dgefa(be.jac, pivot, ierr_linpack); + + if (integrator_rp::linalg_do_pivoting == 1) { + constexpr bool allow_pivot{true}; + dgefa(be.jac, pivot, ierr_linpack); + } else { + constexpr bool allow_pivot{false}; + dgefa(be.jac, pivot, ierr_linpack); + } if (ierr_linpack != 0) { ierr = IERR_LU_DECOMPOSITION_ERROR; break; } - dgesl(be.jac, pivot, b); + if (integrator_rp::linalg_do_pivoting == 1) { + constexpr bool allow_pivot{true}; + dgesl(be.jac, pivot, b); + } else { + constexpr bool allow_pivot{false}; + dgesl(be.jac, pivot, b); + } // update our current guess for the solution diff --git a/integration/VODE/vode_dvjac.H b/integration/VODE/vode_dvjac.H index 35113e8c5c..2975cc5c55 100644 --- a/integration/VODE/vode_dvjac.H +++ b/integration/VODE/vode_dvjac.H @@ -173,7 +173,13 @@ void dvjac (IArray1D& pivot, int& IERPJ, BurnT& state, DvodeT& vstate) RHS::dgefa(vstate.jac); IER = 0; #else - dgefa(vstate.jac, pivot, IER); + if (integrator_rp::linalg_do_pivoting == 1) { + constexpr bool allow_pivot{true}; + dgefa(vstate.jac, pivot, IER); + } else { + constexpr bool allow_pivot{false}; + dgefa(vstate.jac, pivot, IER); + } #endif if (IER != 0) { diff --git a/integration/VODE/vode_dvnlsd.H b/integration/VODE/vode_dvnlsd.H index 5237b96fdf..bf6499490a 100644 --- a/integration/VODE/vode_dvnlsd.H +++ b/integration/VODE/vode_dvnlsd.H @@ -114,7 +114,13 @@ Real dvnlsd (IArray1D& pivot, int& NFLAG, BurnT& state, DvodeT& vstate) #ifdef NEW_NETWORK_IMPLEMENTATION RHS::dgesl(vstate.jac, vstate.y); #else - dgesl(vstate.jac, pivot, vstate.y); + if (integrator_rp::linalg_do_pivoting == 1) { + constexpr bool allow_pivot{true}; + dgesl(vstate.jac, pivot, vstate.y); + } else { + constexpr bool allow_pivot{false}; + dgesl(vstate.jac, pivot, vstate.y); + } #endif if (vstate.RC != 1.0_rt) { diff --git a/integration/_parameters b/integration/_parameters index 6bae04c821..98b60c657f 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -88,3 +88,6 @@ nse_deriv_dt_factor real 0.05 # for NSE update, do we include the weak rate neutrino losses? nse_include_enu_weak integer 1 + +# for the linear algebra, do we allow pivoting? +linalg_do_pivoting integer 1 diff --git a/util/linpack.H b/util/linpack.H index 5fc8c8baf1..d86f91ef10 100644 --- a/util/linpack.H +++ b/util/linpack.H @@ -6,7 +6,7 @@ #include -template +template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void dgesl (RArray2D& a, IArray1D& pivot, RArray1D& b) { @@ -17,11 +17,17 @@ void dgesl (RArray2D& a, IArray1D& pivot, RArray1D& b) // first solve l * y = b if (nm1 >= 1) { for (int k = 1; k <= nm1; ++k) { - int l = pivot(k); - Real t = b(l); - if (l != k) { - b(l) = b(k); - b(k) = t; + + Real t{}; + if constexpr (allow_pivot) { + int l = pivot(k); + t = b(l); + if (l != k) { + b(l) = b(k); + b(k) = t; + } + } else { + t = b(k); } for (int j = k+1; j <= num_eqs; ++j) { @@ -45,7 +51,7 @@ void dgesl (RArray2D& a, IArray1D& pivot, RArray1D& b) -template +template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void dgefa (RArray2D& a, IArray1D& pivot, int& info) { @@ -68,24 +74,29 @@ void dgefa (RArray2D& a, IArray1D& pivot, int& info) // find l = pivot index int l = k; - Real dmax = std::abs(a(k,k)); - for (int i = k+1; i <= num_eqs; ++i) { - if (std::abs(a(i,k)) > dmax) { - l = i; - dmax = std::abs(a(i,k)); + + if constexpr (allow_pivot) { + Real dmax = std::abs(a(k,k)); + for (int i = k+1; i <= num_eqs; ++i) { + if (std::abs(a(i,k)) > dmax) { + l = i; + dmax = std::abs(a(i,k)); + } } - } - pivot(k) = l; + pivot(k) = l; + } // zero pivot implies this column already triangularized if (a(l,k) != 0.0e0_rt) { - // interchange if necessary - if (l != k) { - t = a(l,k); - a(l,k) = a(k,k); - a(k,k) = t; + if constexpr (allow_pivot) { + // interchange if necessary + if (l != k) { + t = a(l,k); + a(l,k) = a(k,k); + a(k,k) = t; + } } // compute multipliers @@ -97,26 +108,30 @@ void dgefa (RArray2D& a, IArray1D& pivot, int& info) // row elimination with column indexing for (int j = k+1; j <= num_eqs; ++j) { t = a(l,j); - if (l != k) { - a(l,j) = a(k,j); - a(k,j) = t; + + if constexpr (allow_pivot) { + if (l != k) { + a(l,j) = a(k,j); + a(k,j) = t; + } } + for (int i = k+1; i <= num_eqs; ++i) { a(i,j) += t * a(i,k); } } - } - else { + } else { info = k; - } } } - pivot(num_eqs) = num_eqs; + if constexpr (allow_pivot) { + pivot(num_eqs) = num_eqs; + } if (a(num_eqs,num_eqs) == 0.0e0_rt) { info = num_eqs;