miscellaneous formatting changes + adjustments to memory management

ORNL · Jun 11, 2024 · 377537a · 377537a
1 parent fdf2af5
commit 377537a
Show file tree

Hide file tree

Showing 12 changed files with 76 additions and 83 deletions.
diff --git a/resolve/LinSolverDirectLUSOL.cpp b/resolve/LinSolverDirectLUSOL.cpp
@@ -56,8 +56,7 @@ namespace ReSolve
     parmlu_[2] = std::pow(std::numeric_limits<real_type>::epsilon(), 0.8);
 
     // TODO: figure out where this exponent comes from :)
-    parmlu_[4] = parmlu_[3] =
-        std::pow(std::numeric_limits<real_type>::epsilon(), 0.67);
+    parmlu_[4] = parmlu_[3] = std::pow(std::numeric_limits<real_type>::epsilon(), 0.67);
 
     parmlu_[5] = 3.0;
     parmlu_[6] = 0.3;
@@ -103,6 +102,11 @@ namespace ReSolve
 
   int LinSolverDirectLUSOL::analyze()
   {
+    // TODO: replace this with something better
+    if (a_ != nullptr || indc_ != nullptr || indr_ != nullptr) {
+        return -1;
+    }
+
     // NOTE: LUSOL does not come with any discrete analysis operation. it is
     //       possible to break apart bits of lu1fac into that, but for now,
     //       we don't bother and shunt it all into ::factorize()
@@ -128,49 +132,47 @@ namespace ReSolve
     index_type* indc_in = A_->getRowData(memory::HOST);
     index_type* indr_in = A_->getColData(memory::HOST);
 
-    a_ = static_cast<real_type*>(std::realloc(a_, lena_ * sizeof(real_type)));
-    indc_ =
-      static_cast<index_type*>(std::realloc(indc_, lena_ * sizeof(index_type)));
-    indr_ =
-      static_cast<index_type*>(std::realloc(indr_, lena_ * sizeof(index_type)));
+    a_ = new real_type[lena_];
+    indc_ = new index_type[lena_];
+    indr_ = new index_type[lena_];
 
     for (index_type i = 0; i < nelem_; i++) {
       a_[i] = a_in[i];
       indc_[i] = indc_in[i] + 1;
       indr_[i] = indr_in[i] + 1;
     }
 
-    p_ = static_cast<index_type*>(std::realloc(p_, m_ * sizeof(index_type)));
+    p_ = new index_type[m_];
     std::fill_n(p_, m_, 0);
 
-    q_ = static_cast<index_type*>(std::realloc(q_, n_ * sizeof(index_type)));
+    q_ = new index_type[n_];
     std::fill_n(q_, n_, 0);
 
-    lenc_ = static_cast<index_type*>(std::realloc(lenc_, n_ * sizeof(index_type)));
+    lenc_ = new index_type[n_];
     std::fill_n(lenc_, n_, 0);
 
-    lenr_ = static_cast<index_type*>(std::realloc(lenr_, m_ * sizeof(index_type)));
+    lenr_ = new index_type[m_];
     std::fill_n(lenr_, m_, 0);
 
-    locc_ = static_cast<index_type*>(std::realloc(locc_, n_ * sizeof(index_type)));
+    locc_ = new index_type[n_];
     std::fill_n(locc_, n_, 0);
 
-    locr_ = static_cast<index_type*>(std::realloc(locr_, m_ * sizeof(index_type)));
+    locr_ = new index_type[m_];
     std::fill_n(locr_, m_, 0);
 
-    iploc_ = static_cast<index_type*>(std::realloc(iploc_, n_ * sizeof(index_type)));
+    iploc_ = new index_type[n_];
     std::fill_n(iploc_, n_, 0);
 
-    iqloc_ = static_cast<index_type*>(std::realloc(iqloc_, m_ * sizeof(index_type)));
+    iqloc_ = new index_type[m_];
     std::fill_n(iqloc_, m_, 0);
 
-    ipinv_ = static_cast<index_type*>(std::realloc(ipinv_, m_ * sizeof(index_type)));
+    ipinv_ = new index_type[m_];
     std::fill_n(ipinv_, m_, 0);
 
-    iqinv_ = static_cast<index_type*>(std::realloc(iqinv_, n_ * sizeof(index_type)));
+    iqinv_ = new index_type[n_];
     std::fill_n(iqinv_, n_, 0);
 
-    w_ = static_cast<real_type*>(std::realloc(w_, n_ * sizeof(real_type)));
+    w_ = new real_type[n_];
     std::fill_n(w_, n_, 0);
 
     return 0;
@@ -209,8 +211,7 @@ namespace ReSolve
 
     // TODO: consider handling inform = 7
 
-    // NOTE: this is probably enough to be handled correctly by most callees
-    return -inform;
+    return inform;
   }
 
   int LinSolverDirectLUSOL::refactorize()
@@ -248,8 +249,7 @@ namespace ReSolve
            locr_,
            &inform);
 
-    // NOTE: ditto above
-    return -inform;
+    return inform;
   }
 
   int LinSolverDirectLUSOL::solve(vector_type* x)

diff --git a/resolve/LinSolverDirectLUSOL.hpp b/resolve/LinSolverDirectLUSOL.hpp
@@ -56,18 +56,20 @@ namespace ReSolve
 
       /// @brief Storage used for the matrices
       ///
-      /// Primary storage area used by LUSOL. Used to hold the nonzeros of matrices
+      /// Primary workspace used by LUSOL. Used to hold the nonzeros of matrices
       /// passed along API boundaries and as a general scratch region
       real_type* a_ = nullptr;
 
-      /// @brief Row indices of matrices passed along API boundaries
+      /// @brief Row data of matrices passed along API boundaries, in addition to
+      ///        functioning as additional workspace storage for LUSOL
       index_type* indc_ = nullptr;
 
-      /// @brief Column indices of matrices passed along API boundaries
+      /// @brief Column data of matrices passed along API boundaries, in addition to
+      ///        functioning as additional workspace storage for LUSOL
       index_type* indr_ = nullptr;
 
       /// @brief The number of nonzero elements within the input matrix, A
-      index_type nelem_ = -1;
+      index_type nelem_ = 0;
 
       /// @brief The length of the dynamically-allocated arrays held within `a_`,
       ///        `indc_`, and `indr_`
@@ -103,13 +105,13 @@ namespace ReSolve
       /// utilizing it will return with inform set to 7, and the intended behavior of
       /// the callee is that they should resize `a_`, `indc_`, and `indr_` to at least
       /// the value specified in `luparm_[12]`
-      index_type lena_ = -1;
+      index_type lena_ = 0;
 
       /// @brief The number of rows in the input matrix, A
-      index_type m_ = -1;
+      index_type m_ = 0;
 
       /// @brief The number of columns in the input matrix, A
-      index_type n_ = -1;
+      index_type n_ = 0;
 
       /// @brief Index-typed parameters passed along the API boundary
       index_type luparm_[30] = {0};

diff --git a/resolve/matrix/Coo.hpp b/resolve/matrix/Coo.hpp
@@ -33,9 +33,7 @@ namespace ReSolve { namespace matrix {
       virtual int expand();
 
       // NOTE: see comment in resolve/matrix/Utilities.cpp
-      virtual std::function<
-          std::tuple<std::tuple<index_type, index_type, real_type>, bool>()>
-          elements(memory::MemorySpace);
+      virtual std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> elements(memory::MemorySpace);
   };
 
 }} // namespace ReSolve::matrix
diff --git a/resolve/matrix/Csc.cpp b/resolve/matrix/Csc.cpp
@@ -263,14 +263,11 @@ namespace ReSolve
     return -1;
   }
 
-  std::function<
-      std::tuple<std::tuple<index_type, index_type, real_type>, bool>()>
-  matrix::Csc::elements(memory::MemorySpace _)
+  std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> matrix::Csc::elements(memory::MemorySpace _)
   {
     out::error() << "Sparse::elements(memory::MemorySpace) is not implemented"
                     "for Csc!\n";
-    return []() -> std::tuple<std::tuple<index_type, index_type, real_type>,
-                              bool> {
+    return []() -> std::tuple<std::tuple<index_type, index_type, real_type>, bool> {
       return {{0, 0, 0}, false};
     };
   }

diff --git a/resolve/matrix/Csc.hpp b/resolve/matrix/Csc.hpp
@@ -33,9 +33,7 @@ namespace ReSolve { namespace matrix {
       virtual int expand();
 
       // NOTE: see comment in resolve/matrix/Utilities.cpp
-      virtual std::function<
-          std::tuple<std::tuple<index_type, index_type, real_type>, bool>()>
-          elements(memory::MemorySpace);
+      virtual std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> elements(memory::MemorySpace);
   };
 
 }} // namespace ReSolve::matrix
diff --git a/resolve/matrix/Csr.cpp b/resolve/matrix/Csr.cpp
@@ -617,9 +617,7 @@ namespace ReSolve
     return 0;
   }
 
-  std::function<
-      std::tuple<std::tuple<index_type, index_type, real_type>, bool>()>
-  matrix::Csr::elements(memory::MemorySpace memory_space)
+  std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> matrix::Csr::elements(memory::MemorySpace memory_space)
   {
     std::shared_ptr<index_type> i(new index_type(0)), j(new index_type(0));
     std::shared_ptr<index_type> n(new index_type(getNnz()));
@@ -628,15 +626,12 @@ namespace ReSolve
     real_type* values = getValues(memory_space);
 
     if (rows == nullptr || columns == nullptr || values == nullptr) {
-      return []() -> std::tuple<std::tuple<index_type, index_type, real_type>,
-                                bool> {
+      return []() -> std::tuple<std::tuple<index_type, index_type, real_type>, bool> {
         return {{0, 0, 0}, false};
       };
     }
 
-    return
-        [=]()
-            -> std::tuple<std::tuple<index_type, index_type, real_type>, bool> {
+    return [=]() -> std::tuple<std::tuple<index_type, index_type, real_type>, bool> {
           if (*j == *n) {
             return {{0, 0, 0}, false};
           }

diff --git a/resolve/matrix/Csr.hpp b/resolve/matrix/Csr.hpp
@@ -54,9 +54,7 @@ namespace ReSolve { namespace matrix {
         virtual int expand();
 
         // NOTE: see comment in resolve/matrix/Utilities.cpp
-        virtual std::function<
-            std::tuple<std::tuple<index_type, index_type, real_type>, bool>()>
-            elements(memory::MemorySpace);
+        virtual std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> elements(memory::MemorySpace);
 
       int coo2csr(matrix::Coo* mat, memory::MemorySpace memspace);
   };

diff --git a/resolve/matrix/Sparse.hpp b/resolve/matrix/Sparse.hpp
@@ -108,13 +108,9 @@ namespace ReSolve { namespace matrix {
       friend int matrix::expand(matrix::Sparse*);
 
       // NOTE: temporary. see comment in resolve/matrix/Utilities.cpp
-      virtual std::function<
-          std::tuple<std::tuple<index_type, index_type, real_type>, bool>()>
-          elements(memory::MemorySpace) = 0;
+      virtual std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> elements(memory::MemorySpace) = 0;
 
       // NOTE: ditto (same note on `virtual void elements()`)
-      friend std::function<
-          std::tuple<std::tuple<index_type, index_type, real_type>, bool>()>
-      matrix::elements(matrix::Sparse*);
+      friend std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> matrix::elements(matrix::Sparse*);
   };
 }} // namespace ReSolve::matrix
diff --git a/resolve/matrix/Utilities.cpp b/resolve/matrix/Utilities.cpp
@@ -25,9 +25,7 @@ namespace ReSolve
 
   // NOTE: i don't like these too much. they're just "good enough"
 
-  std::function<
-      std::tuple<std::tuple<index_type, index_type, real_type>, bool>()>
-  matrix::elements(matrix::Sparse* A)
+  std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> matrix::elements(matrix::Sparse* A)
   {
     return A->elements(memory::HOST);
   }

diff --git a/resolve/matrix/Utilities.hpp b/resolve/matrix/Utilities.hpp
@@ -29,8 +29,6 @@ namespace ReSolve
     ///
     /// This interfaced is to be used as a bandage. It will be removed in a
     /// later version
-    std::function<
-        std::tuple<std::tuple<index_type, index_type, real_type>, bool>()>
-    elements(Sparse*);
+    std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> elements(Sparse*);
   } // namespace matrix
 } // namespace ReSolve
diff --git a/tests/unit/matrix/LUSOLTests.hpp b/tests/unit/matrix/LUSOLTests.hpp
@@ -92,15 +92,33 @@ namespace ReSolve
         //     [     0      0      2      0      2      0      3      3      0]
         //     [     2      0      0      0      0      0      0      5      1]
         //     [     0      0      7      0      8      0      0      0      4]
-        std::vector<index_type> rowsA_ = {0, 0, 0, 1, 1, 1, 2, 2, 2,
-                                          3, 3, 3, 4, 5, 5, 5, 5, 6,
-                                          6, 6, 6, 7, 7, 7, 8, 8, 8};
-        std::vector<index_type> colsA_ = {0, 4, 6, 1, 3, 5, 0, 4, 7,
-                                          3, 5, 8, 0, 1, 3, 5, 6, 2,
-                                          4, 6, 7, 0, 7, 8, 2, 4, 8};
-        std::vector<real_type> valsA_ = {2., 1., 3., 7., 5., 4., 1., 3., 2.,
-                                         3., 2., 8., 1., 4., 5., 1., 6., 2.,
-                                         2., 3., 3., 2., 5., 1., 7., 8., 4.};
+        std::vector<index_type> rowsA_ = {0, 0, 0,
+                                          1, 1, 1,
+                                          2, 2, 2,
+                                          3, 3, 3,
+                                          4,
+                                          5, 5, 5, 5,
+                                          6, 6, 6, 6,
+                                          7, 7, 7,
+                                          8, 8, 8};
+        std::vector<index_type> colsA_ = {0, 4, 6,
+                                          1, 3, 5,
+                                          0, 4, 7,
+                                          3, 5, 8,
+                                          0,
+                                          1, 3, 5, 6,
+                                          2, 4, 6, 7,
+                                          0, 7, 8,
+                                          2, 4, 8};
+        std::vector<real_type> valsA_ = {2., 1., 3.,
+                                         7., 5., 4.,
+                                         1., 3., 2.,
+                                         3., 2., 8.,
+                                         1.,
+                                         4., 5., 1., 6.,
+                                         2., 2., 3., 3.,
+                                         2., 5., 1.,
+                                         7., 8., 4.};
 
         /// @brief Creates a test matrix
         matrix::Coo* createMatrix()

diff --git a/tests/unit/matrix/MatrixExpansionTests.hpp b/tests/unit/matrix/MatrixExpansionTests.hpp
@@ -65,25 +65,20 @@ namespace ReSolve
                                    {0, 4, 3.0},
                                    {4, 0, 3.0}};
 
-        bool validateAnswer(
-            matrix::Coo* A,
-            std::vector<std::tuple<index_type, index_type, real_type>> target)
+        bool validateAnswer(matrix::Coo* A,
+                            std::vector<std::tuple<index_type, index_type, real_type>> target)
         {
           return validateAnswer(matrix::elements(A), target);
         }
 
-        bool validateAnswer(
-            matrix::Csr* A,
-            std::vector<std::tuple<index_type, index_type, real_type>> target)
+        bool validateAnswer(matrix::Csr* A,
+                            std::vector<std::tuple<index_type, index_type, real_type>> target)
         {
           return validateAnswer(matrix::elements(A), target);
         }
 
-        bool validateAnswer(
-            std::function<
-                std::tuple<std::tuple<index_type, index_type, real_type>,
-                           bool>()> f,
-            std::vector<std::tuple<index_type, index_type, real_type>> target)
+        bool validateAnswer(std::function<std::tuple<std::tuple<index_type, index_type, real_type>, bool>()> f,
+                            std::vector<std::tuple<index_type, index_type, real_type>> target)
         {
           bool ok;
           std::tuple<index_type, index_type, real_type> t;