#include <thread>
#include <random>
#include <vector>
#include <cstdlib>
#include <hpc/matvec/gematrix.h>
#include <hpc/matvec/apply.h>
#include <hpc/matvec/print.h>
template<typename T>
struct Slices {
Slices(T nof_threads, T problem_size) :
nof_threads((assert(nof_threads > 0), nof_threads)),
problem_size(problem_size),
remainder(problem_size % nof_threads),
slice_size(problem_size / nof_threads) {
}
T offset(T index) {
assert(index < nof_threads);
if (index < remainder) {
return index * (slice_size + 1);
} else {
return remainder * (slice_size + 1) +
(index - remainder) * slice_size;
}
}
T size(T index) {
assert(index < nof_threads);
if (index < remainder) {
return slice_size + 1;
} else {
return slice_size;
}
}
T nof_threads; T problem_size;
T remainder; T slice_size;
};
template<typename MA>
typename std::enable_if<hpc::matvec::IsRealGeMatrix<MA>::value, void>::type
randomInit(MA& A) {
using ElementType = typename MA::ElementType;
using Index = typename MA::Index;
std::random_device random;
std::mt19937 mt(random());
std::uniform_real_distribution<ElementType> uniform(-100,100);
hpc::matvec::apply(A, [&](ElementType& val, Index i, Index j) -> void {
val = uniform(mt);
});
}
int main() {
using namespace hpc::matvec;
GeMatrix<double> A(51, 7);
unsigned int nof_threads = std::thread::hardware_concurrency();
std::vector<std::thread> threads(nof_threads);
Slices<GeMatrix<double>::Index> slices(nof_threads, A.numRows);
for (int index = 0; index < nof_threads; ++index) {
auto firstRow = slices.offset(index);
auto numRows = slices.size(index);
auto A_ = A(firstRow, 0, numRows, A.numCols);
threads[index] = std::thread([=]() mutable { randomInit(A_); });
}
for (int index = 0; index < nof_threads; ++index) {
threads[index].join();
}
/* print a small block of each of the initialized matrices */
print(A, "A");
}