#include <blaze/Math.h>
#include <iostream>
#include <random>
#include "gemm.h"

// fill rectangular matrix with random values
template <typename MATRIX>
void
fill(MATRIX &A)
{
    typedef typename MATRIX::ElementType  T;

    std::random_device                  random;
    std::default_random_engine          mt(random());
    std::uniform_real_distribution<T>   uniform(-100,100);

    for (std::size_t i=0; i<(~A).rows(); ++i) {
        for (std::size_t j=0; j<(~A).columns(); ++j) {
            A(i,j) = uniform(mt);
        }
    }
}

int
main()
{
    typedef double ElementType;

    constexpr auto SOA = blaze::rowMajor;
    constexpr auto SOB = blaze::rowMajor;
    constexpr auto SOC = blaze::rowMajor;

    std::size_t m = 5;
    std::size_t n = 5;
    std::size_t k = 5;

    ElementType alpha = 1;
    ElementType beta  = 0;


    blaze::DynamicMatrix<double, SOA>   A(m, k);
    blaze::DynamicMatrix<double, SOB>   B(k, n);
    blaze::DynamicMatrix<double, SOC>   C(m, n);

    fill(A);
    fill(B);
    fill(C);

    std::cout << "A = " << A << std::endl;
    std::cout << "B = " << B << std::endl;
    std::cout << "C = " << C << std::endl;

    std::cout << "C <- " << beta << "*C + " << alpha << "*A*B"
              << std::endl;
    foo::gemm(alpha, A, B, beta, C);
    std::cout << "C = " << C << std::endl;

    std::cout << "C <- " << beta << "*C + " << alpha << "*(A^T+A)*(2*B)"
              << std::endl;
    foo::gemm(alpha, blaze::trans(A)+A, 2*B, beta, C);
    std::cout << "C = " << C << std::endl;
}