#include <cstdio>
#include <hpc/matvec/gematrix.h>
#include <hpc/matvec/apply.h>
#include <hpc/matvec/print.h>
#include <hpc/matvec/asum.h>
#include <hpc/aux/slices.h>
#include <hpc/aux/repool.h>

template<typename MA, typename Pool>
typename std::enable_if<hpc::matvec::IsRealGeMatrix<MA>::value, void>::type
randomInit(MA& A, Pool& pool) {
   using ElementType = typename MA::ElementType;
   using Index = typename MA::Index;
   using EngineType = typename Pool::EngineType;

   std::uniform_real_distribution<double> uniform(-100, 100);
   hpc::aux::RandomEngineGuard<EngineType> guard(pool);
   auto& engine(guard.get());

   hpc::matvec::apply(A, [&](ElementType& val, Index i, Index j) -> void {
      val = uniform(engine);
   });
}

template<typename MA, typename RandomEnginePool>
typename std::enable_if<hpc::matvec::IsRealGeMatrix<MA>::value, void>::type
mt_randomInit(MA& A, RandomEnginePool& repool) {
   using ElementType = typename MA::ElementType;
   using Index = typename MA::Index;
   using namespace hpc::aux;

   #pragma omp parallel for
   for (Index row = 0; row < A.numRows; ++row) {
      auto A_ = A(row, 0, 1, A.numCols);
      randomInit(A_, repool);
   }
}

template<typename MA>
auto mt_asum(MA& A) -> decltype(asum(A)) {
   using ElementType = typename MA::ElementType;
   using Index = typename MA::Index;
   using namespace hpc::aux;
   using Result = decltype(asum(A));

   Result sum = 0;
   #pragma omp parallel for reduction(+:sum)
   for (Index row = 0; row < A.numRows; ++row) {
      auto A_ = A(row, 0, 1, A.numCols);
      sum += asum(A_);
   }
   return sum;
}

int main() {
   using namespace hpc::matvec;
   using namespace hpc::aux;
   RandomEnginePool<std::mt19937> repool(4);

   GeMatrix<double> A(51, 7);
   mt_randomInit(A, repool);
   std::printf("mt_asum: %lg\n", mt_asum(A));
   std::printf("asum: %lg\n", asum(A));
}