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
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
#include <condition_variable>
#include <cstdlib>
#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;
      RandomEnginePool(std::size_t size) :
	    size(size), nof_free_engines(size),
	    inuse(size), engines(size) {
	 std::random_device r;
	 for (std::size_t i = 0; i < size; ++i) {
	    engines[i].seed(r()); inuse[i] = false;
	 }
      }
      T& get() {
	 std::unique_lock<std::mutex> lock(mutex);
	 while (nof_free_engines == 0) {
	    cv.wait(lock);
	 }
	 for (std::size_t i = 0; i < size; ++i) {
	    if (!inuse[i]) {
	       inuse[i] = true; --nof_free_engines;
	       return engines[i];
	    }
	 }
	 assert(false);
	 #pragma GCC diagnostic ignored "-Wreturn-type"
      }
      void free(T& engine) {
	 {
	    std::unique_lock<std::mutex> lock(mutex);
	    bool found = false;
	    for (std::size_t i = 0; i < size; ++i) {
	       if (&engine == &engines[i]) {
		  inuse[i] = false; ++nof_free_engines;
		  found = true; break;
	       }
	    }
	    assert(found);
	 }
	 cv.notify_one();
      }
   private:
      std::mutex mutex;
      std::condition_variable cv;
      std::size_t size;
      std::size_t nof_free_engines;
      std::vector<bool> inuse;
      std::vector<T> engines;
};

template<typename T>
struct RandomEngineGuard {
   using EngineType = T;
   RandomEngineGuard(RandomEnginePool<T>& pool) :
      pool(pool), engine(pool.get()) {
   }
   ~RandomEngineGuard() {
      pool.free(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(2);
   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");
}