From 42fcc6e0954bd2f7f2fded2a1d0eef4f226443d7 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sat, 16 Dec 2023 17:40:06 -0500
Subject: [PATCH] Revert "Start moving towards zero-based indexing with linpack
 (#536)"

This reverts commit 1f0c666f18c0418a3ab0aa817de00f8fd52ebac7.
---
 util/linpack.H | 39 +++++++++++++++++----------------------
 1 file changed, 17 insertions(+), 22 deletions(-)

diff --git a/util/linpack.H b/util/linpack.H
index bdb826a915..5fc8c8baf1 100644
--- a/util/linpack.H
+++ b/util/linpack.H
@@ -8,38 +8,35 @@
 
 template <int num_eqs>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void dgesl (RArray2D& a1, IArray1D& pivot1, RArray1D& b1)
+void dgesl (RArray2D& a, IArray1D& pivot, RArray1D& b)
 {
-    auto const& a = reinterpret_cast<Array2D<Real, 0, num_eqs-1, 0, num_eqs-1>&>(a1);
-    auto const& pivot = reinterpret_cast<Array1D<short, 0, num_eqs-1>&>(pivot1);
-    auto& b = reinterpret_cast<Array1D<Real, 0, num_eqs-1>&>(b1);
 
     int nm1 = num_eqs - 1;
 
     // solve a * x = b
     // first solve l * y = b
     if (nm1 >= 1) {
-        for (int k = 0; k < nm1; ++k) {
-            int l = pivot(k) - 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;
             }
 
-            for (int j = k+1; j < num_eqs; ++j) {
+            for (int j = k+1; j <= num_eqs; ++j) {
                 b(j) += t * a(j,k);
             }
         }
     }
 
     // now solve u * x = y
-    for (int kb = 0; kb < num_eqs; ++kb) {
+    for (int kb = 1; kb <= num_eqs; ++kb) {
 
-        int k = num_eqs - kb - 1;
+        int k = num_eqs + 1 - kb;
         b(k) = b(k) / a(k,k);
         Real t = -b(k);
-        for (int j = 0; j < k; ++j) {
+        for (int j = 1; j <= k-1; ++j) {
             b(j) += t * a(j,k);
         }
     }
@@ -50,10 +47,8 @@ void dgesl (RArray2D& a1, IArray1D& pivot1, RArray1D& b1)
 
 template <int num_eqs>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void dgefa (RArray2D& a1, IArray1D& pivot1, int& info)
+void dgefa (RArray2D& a, IArray1D& pivot, int& info)
 {
-    auto& a = reinterpret_cast<Array2D<Real, 0, num_eqs-1, 0, num_eqs-1>&>(a1);
-    auto& pivot = reinterpret_cast<Array1D<short, 0, num_eqs-1>&>(pivot1);
 
     // dgefa factors a matrix by gaussian elimination.
     // a is returned in the form a = l * u where
@@ -69,19 +64,19 @@ void dgefa (RArray2D& a1, IArray1D& pivot1, int& info)
 
     if (nm1 >= 1) {
 
-        for (int k = 0; k < nm1; ++k) {
+        for (int k = 1; k <= nm1; ++k) {
 
             // find l = pivot index
             int l = k;
             Real dmax = std::abs(a(k,k));
-            for (int i = k+1; i < num_eqs; ++i) {
+            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 + 1;
+            pivot(k) = l;
 
             // zero pivot implies this column already triangularized
             if (a(l,k) != 0.0e0_rt) {
@@ -95,25 +90,25 @@ void dgefa (RArray2D& a1, IArray1D& pivot1, int& info)
 
                 // compute multipliers
                 t = -1.0e0_rt / a(k,k);
-                for (int j = k+1; j < num_eqs; ++j) {
+                for (int j = k+1; j <= num_eqs; ++j) {
                     a(j,k) *= t;
                 }
 
                 // row elimination with column indexing
-                for (int j = k+1; j < num_eqs; ++j) {
+                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;
                     }
-                    for (int i = k+1; i < num_eqs; ++i) {
+                    for (int i = k+1; i <= num_eqs; ++i) {
                         a(i,j) += t * a(i,k);
                     }
                 }
             }
             else {
 
-                info = k+1;
+                info = k;
 
             }
 
@@ -121,9 +116,9 @@ void dgefa (RArray2D& a1, IArray1D& pivot1, int& info)
 
     }
 
-    pivot(num_eqs-1) = num_eqs;
+    pivot(num_eqs) = num_eqs;
 
-    if (a(num_eqs-1,num_eqs-1) == 0.0e0_rt) {
+    if (a(num_eqs,num_eqs) == 0.0e0_rt) {
         info = num_eqs;
     }