1
      2
      3
      4
      5
      6
      7
      8
      9
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
     20
     21
     22
     23
     24
     25
     26
     27
     28
     29
     30
     31
     32
     33
     34
     35
     36
     37
     38
     39
     40
     41
     42
     43
     44
     45
     46
     47
     48
     49
     50
     51
     52
     53
     54
     55
     56
     57
     58
     59
     60
     61
     62
     63
     64
     65
     66
     67
     68
     69
     70
     71
     72
     73
     74
     75
     76
     77
     78
     79
     80
     81
     82
     83
     84
     85
     86
     87
     88
     89
     90
     91
     92
     93
     94
     95
     96
#include <cstdlib>
#include <deque>
#include <mutex>
#include <random>
#include <thread>
#include <vector>
#include <hpc/aux/slices.hpp>
#include <hpc/matvec/gematrix.hpp>
#include <hpc/matvec/iterators.hpp>
#include <hpc/matvec/print.hpp>

using namespace hpc;
using namespace hpc::aux;
using namespace hpc::matvec;

template<typename T>
struct RandomEnginePool {
      using EngineType = T;
      T get() {
         /* check if we have a free engine in the unused deque */
         {
            std::lock_guard<std::mutex> lock(mutex);
            if (unused.size() > 0) {
               T rg = std::move(unused.front());
               unused.pop_front();
               return rg;
            }
         }
         /* prepare new random generator */
         return T(r());
      }
      void free(T&& engine) {
         std::lock_guard<std::mutex> lock(mutex);
         unused.push_back(engine);
      }
   private:
      std::mutex mutex;
      std::random_device r;
      std::deque<T> unused;
};

template<typename T>
struct RandomEngineGuard {
   using EngineType = T;
   RandomEngineGuard(RandomEnginePool<T>& pool) :
      pool(pool), engine(pool.get()) {
   }
   ~RandomEngineGuard() {
      pool.free(std::move(engine));
   }
   T& get() {
      return engine;
   }
   RandomEnginePool<T>& pool;
   T engine;
};

template <
   template<typename> class MatrixA, typename T,
   typename POOL,
   Require< Ge<MatrixA<T>> > = true
>
void randomInit(MatrixA<T>& A, POOL& pool) {
   using EngineType = typename POOL::EngineType;
   RandomEngineGuard<EngineType> guard(pool);
   std::uniform_real_distribution<double> uniform(-100, 100);
   auto& engine = guard.get();

   for (auto [i, j, Aij]: A) {
      Aij = uniform(engine);
      (void) i; (void) j; // suppress gcc warning
   }
}

int main() {
   RandomEnginePool<std::mt19937> pool;
   GeMatrix<double> A(51, 7);
   std::size_t nof_threads = std::thread::hardware_concurrency();

   std::vector<std::thread> threads(nof_threads);
   UniformSlices<std::size_t> slices(nof_threads, A.numRows());
   for (std::size_t index = 0; index < nof_threads; ++index) {
      auto firstRow = slices.offset(index);
      auto numRows = slices.size(index);
      threads[index] = std::thread([
	       A_ = A.view(firstRow, 0, numRows, A.numCols()),
	       &pool
	    ]() mutable {
	 randomInit(A_, pool);
      });
   }
   for (std::size_t index = 0; index < nof_threads; ++index) {
      threads[index].join();
   }
   print(A, " %7.2f");
}