#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(" %4.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");
   }
}

int main() {
   using namespace matvec;
   GeMatrix<double> A(3, 7, StorageOrder::ColMajor);
   initGeMatrix(A, init_value<double, std::size_t>);
   printf("A:\n");
   print_matrix(A);
   GeMatrix<std::complex<double>> B(2, 3, StorageOrder::ColMajor);
   initGeMatrix(B, init_complex_value<double, std::size_t>);
   printf("B:\n");
   print_matrix(B);
}