#include <cstdlib>
#include <cstdio>
#include "gematrix.h"

template<typename T, typename Index>
T init_value(Index i, Index j, Index m, Index n) {
   return T(j)*T(n) + T(i) + T(1);
}

template<typename T, typename Index>
std::complex<T> init_complex_value(Index i, Index j, Index m, Index n) {
   return std::complex<T>(T(i), T(j));
}

void print_value(double value) {
   std::printf(" %10.1lf", value);
}

void print_value(std::complex<double> value) {
   std::printf(" (%4.1lf, %4.1lf)", value.real(), value.imag());
}

template<typename Matrix>
void print_matrix(const Matrix& A) {
   for (std::size_t i = 0; i < A.m; ++i) {
      std::printf("  ");
      for (std::size_t j = 0; j < A.n; ++j) {
	 print_value(A(i, j));
      }
      std::printf("\n");
   }
}

template<typename T, typename Index>
struct IncreasingRandomValues {
   std::mt19937 mt;
   std::uniform_real_distribution<T> uniform;
   T last_value;
   
   IncreasingRandomValues() :
      mt(std::random_device()()), uniform(0, 1<<16), last_value(0) {
   }
   IncreasingRandomValues(unsigned int seed) :
      mt(seed), uniform(1, 1<<16) {
   }
   T operator()(Index i, Index j) {
      last_value += uniform(mt);
      return last_value;
   }
};

int main() {
   using namespace matvec;
   GeMatrix<double> A(3, 7, StorageOrder::ColMajor);
   {
      std::random_device random;
      std::mt19937 mt(random());
      std::uniform_real_distribution<double> uniform(1, 1<<16);
      double last_value = 0;
      initGeMatrix(A, [=,&last_value]
	    (std::size_t i, std::size_t j) mutable -> double {
	 last_value += uniform(mt);
	 return last_value;
      });
   }
   printf("A:\n");
   print_matrix(A);

   GeMatrix<double> B(3, 7, StorageOrder::ColMajor);
   {
      std::random_device random;
      std::mt19937 mt(random());
      std::uniform_real_distribution<double> uniform(1, 1<<16);
      double last_value = 0;
      applyGeMatrix(B, [=,&last_value]
	    (double& Bij, std::size_t i, std::size_t j) mutable {
	 last_value += uniform(mt);
	 Bij = last_value;
      });
   }
   printf("B:\n");
   print_matrix(B);
   {
      double sum = 0;
      applyGeMatrix(B, [&sum]
	    (double Bij, std::size_t i, std::size_t j) mutable {
	 sum += Bij;
      });
      printf("sum of all elements of B: %10.2lf\n", sum);
   }
   {
      double amax = 0;
      applyGeMatrix(B, [&amax]
	    (double& Bij, std::size_t i, std::size_t j) mutable {
	 if (abs(Bij) > amax) amax = abs(Bij);
      });
      printf("absolute max of B: %10.2lf\n", amax);
   }
   GeMatrixCombinedConstView<GeMatrix<double>, GeMatrix<double>>
      C(A, B, [](double a, double b) -> double {
	 return abs(a - b);
      });
   printf("C:\n");
   print_matrix(C);
   {
      double sum = 0;
      applyGeMatrix(C, [&sum]
	    (double Cij, std::size_t i, std::size_t j) mutable {
	 sum += Cij;
      });
      printf("sum of all elements of C: %10.2lf\n", sum);
   }
}