Shaked/random error fix

CHANGELOG.md

@@ -4,6 +4,10 @@
 
 - Added classes and tests for permutation, Ruiz scaling, Cholesky factorization, Schur complement conjugate gradient and matrix multiplication and addition.
 
+- Changed random number generation int tests to be C++ style and fixed-seed, to avoid random failures.
+
+- Added Schur Complement Conjugate Gradient class.
+
 ## Changes to Re::Solve since release 0.99.2
 
 - Added cmake-format.

resolve/Common.hpp

@@ -14,11 +14,12 @@ namespace ReSolve
 
   namespace constants
   {
-    constexpr real_type ZERO      = 0.0;
-    constexpr real_type ONE       = 1.0;
-    constexpr real_type TWO       = 2.0;
-    constexpr real_type HALF      = 0.5;
-    constexpr real_type MINUS_ONE = -1.0;
+    constexpr real_type  ZERO      = 0.0;
+    constexpr real_type  ONE       = 1.0;
+    constexpr real_type  TWO       = 2.0;
+    constexpr real_type  HALF      = 0.5;
+    constexpr real_type  MINUS_ONE = -1.0;
+    constexpr index_type SEED      = 12345;
 
     constexpr real_type MACHINE_EPSILON = std::numeric_limits<real_type>::epsilon();
   } // namespace constants

resolve/hykkt/cholesky/CholeskySolverHip.cpp

@@ -60,6 +60,7 @@ namespace ReSolve
       }
       A_chol_ = convertToCholmod(A);
       A_      = A;
+      mem_.deviceSynchronize();
     }
 
     void CholeskySolverHip::symbolicAnalysis()
@@ -70,6 +71,7 @@ namespace ReSolve
       {
         out::error() << "Cholesky symbolic analysis failed with status: " << Common_.status << "\n";
       }
+      mem_.deviceSynchronize();
     }
 
     /**
@@ -156,6 +158,7 @@ namespace ReSolve
           out::error() << "Refactorization step failed with status: " << status << "\n";
         }
         L_->setUpdated(memory::DEVICE);
+        mem_.deviceSynchronize();
       }
     }
 

tests/unit/hykkt/HykktCholeskyTests.hpp

@@ -4,6 +4,7 @@
 #include <algorithm>
 #include <cholmod.h>
 #include <iterator>
+#include <random>
 #include <sstream>
 #include <string>
 #include <vector>
@@ -26,8 +27,8 @@ namespace ReSolve
     class HykktCholeskyTests : public TestBase
     {
     public:
-      HykktCholeskyTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler)
-        : memspace_(memspace), matrixHandler_(matrixHandler)
+      HykktCholeskyTests(memory::MemorySpace memspace, MatrixHandler& matrixHandler, std::mt19937& generator)
+        : memspace_(memspace), matrixHandler_(matrixHandler), generator_(generator)
       {
         cholmod_start(&Common);
       }
@@ -127,7 +128,6 @@ namespace ReSolve
         solver.addMatrixInfo(A);
         solver.symbolicAnalysis();
         solver.numericalFactorization();
-
         // Generate a random vector x_expected and compute b = A * x_expected
         vector::Vector* x_expected = randomVector(n);
 
@@ -242,7 +242,8 @@ namespace ReSolve
           // reset values
           for (size_t j = 0; j < L->nzmax; j++)
           {
-            static_cast<double*>(L->x)[j] = 2.0 * static_cast<double>(rand()) / RAND_MAX - 1.0;
+            std::uniform_real_distribution<double> distribution(-1.0, 1.0);
+            static_cast<double*>(L->x)[j] = 2.0 * distribution(generator_) - 1.0;
           }
           cholmod_free_sparse(&L_tr, &Common);
           cholmod_free_sparse(&A_chol, &Common);
@@ -266,6 +267,7 @@ namespace ReSolve
     private:
       ReSolve::memory::MemorySpace memspace_;
       MatrixHandler&               matrixHandler_;
+      std::mt19937&                generator_;
 
       cholmod_common Common;
 
@@ -281,14 +283,16 @@ namespace ReSolve
           L_p[i + 1] = L_p[i];
           for (size_t j = i; j < n; ++j)
           {
-            // with probability 'density', add a non-zero entry
-            if (i == j || static_cast<double>(rand()) / RAND_MAX < density)
+            std::uniform_real_distribution<double> rand_dist(0.0, 1.0);
+            if (i == j || rand_dist(generator_) < density)
             {
               L_i.push_back((int) j);
               // force diagonal entry to be non-zero
               double value = 2.0;
               if (i != j)
               {
+                std::uniform_real_distribution<double> value_dist(-1.0, 1.0);
+                value = 2.0 * value_dist(generator_);
                 value = 2.0 * static_cast<double>(rand()) / RAND_MAX - 1.0;
               }
               L_x.push_back(value);
@@ -311,9 +315,10 @@ namespace ReSolve
       {
         vector::Vector* v = new vector::Vector(n);
         v->allocate(memory::HOST);
+        std::uniform_real_distribution<double> distribution(0.0, 1.0);
         for (index_type i = 0; i < n; ++i)
         {
-          v->getData(memory::HOST)[i] = static_cast<double>(rand()) / RAND_MAX;
+          v->getData(memory::HOST)[i] = distribution(generator_);
         }
         v->setDataUpdated(memory::HOST);
         if (memspace_ == memory::DEVICE)

tests/unit/hykkt/HykktSCCGTests.hpp

@@ -6,6 +6,7 @@
 #pragma once
 
 #include <filesystem>
+#include <random>
 
 #include <resolve/MemoryUtils.hpp>
 #include <resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp>
@@ -35,9 +36,10 @@ namespace ReSolve
        * @param[in] memspace Memory space for the test (HOST or DEVICE).
        * @param[in] matrix_handler Reference to a matrix handler for the selected backend.
        * @param[in] vector_handler Reference to a vector handler for the selected backend.
+       * @param[in] generator Reference to a C++ random number generator.
        */
-      HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrix_handler, VectorHandler& vector_handler)
-        : memspace_(memspace), matrix_handler_(matrix_handler), vector_handler_(vector_handler)
+      HykktSchurComplementConjugateGradientTests(memory::MemorySpace memspace, MatrixHandler& matrix_handler, VectorHandler& vector_handler, std::mt19937& generator)
+        : memspace_(memspace), matrix_handler_(matrix_handler), vector_handler_(vector_handler), generator_(generator)
       {
       }
 
@@ -130,16 +132,18 @@ namespace ReSolve
       memory::MemorySpace memspace_;       ///< Memory space used by the test.
       MatrixHandler&      matrix_handler_; ///< Backend-specific matrix handler.
       VectorHandler&      vector_handler_; ///< Backend-specific vector handler.
+      std::mt19937&       generator_;      ///< C++ random number generator.
 
       /**
-       * @brief Generate a random vector of doubles between 0 and RAND_MAX. Copied from HykktCholeskyTests.hpp.
+       * @brief Generate a random vector of doubles between 0 and 1. Copied from HykktCholeskyTests.hpp.
        * @param[in] vec Target vector to write to.
        */
       void randomVector(vector::Vector* vec)
       {
+        std::uniform_real_distribution<double> distribution(0.0, 1.0);
         for (index_type i = 0; i < vec->getSize(); ++i)
         {
-          vec->getData(memory::HOST)[i] = static_cast<double>(rand()) / RAND_MAX;
+          vec->getData(memory::HOST)[i] = distribution(generator_);
         }
         vec->setDataUpdated(memory::HOST);
         if (memspace_ == memory::DEVICE)

tests/unit/hykkt/runHykktCholeskyTests.cpp

@@ -7,8 +7,10 @@
  */
 #include <fstream>
 #include <iostream>
+#include <random>
 #include <string>
 
+#include "resolve/Common.hpp"
 #include "tests/unit/hykkt/HykktCholeskyTests.hpp"
 
 /**
@@ -24,9 +26,9 @@ void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace,
 
   WorkspaceType workspace;
   workspace.initializeHandles();
-  ReSolve::MatrixHandler handler(&workspace);
-
-  ReSolve::tests::HykktCholeskyTests test(memspace, handler);
+  ReSolve::MatrixHandler             handler(&workspace);
+  std::mt19937                       generator(ReSolve::constants::SEED); // set random seed for reproducibility
+  ReSolve::tests::HykktCholeskyTests test(memspace, handler, generator);
 
   result += test.minimalCorrectness();
   handler.setValuesChanged(true, memspace);

tests/unit/hykkt/runHykktSCCGTests.cpp

@@ -5,6 +5,7 @@
  */
 #include <fstream>
 #include <iostream>
+#include <random>
 #include <string>
 
 #include <resolve/matrix/MatrixHandler.hpp>
@@ -35,7 +36,8 @@ void runTests(const std::string& backend, ReSolve::memory::MemorySpace memspace,
   workspace.initializeHandles();
   ReSolve::MatrixHandler                                     matrix_handler(&workspace);
   ReSolve::VectorHandler                                     vector_handler(&workspace);
-  ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, matrix_handler, vector_handler);
+  std::mt19937                                               generator(ReSolve::constants::SEED);
+  ReSolve::tests::HykktSchurComplementConjugateGradientTests test(memspace, matrix_handler, vector_handler, generator);
 
   result += test.SCCGTest();