#include <iostream>
#include <iomanip>

#define METHOD_OPERATOR

// Matrix is a template of its dimensions and the element type,
// which if omitted defaults to double. If Columns is omitted,
// it defaults to Rows, i.e. a square matrix.
template <int Rows, int Columns = Rows, class Ring = double>
class Matrix
{
  Ring Elements[Rows][Columns];

public:

  // default constructor, gives a zero matrix
  Matrix (void)
  {
    for (int i = 0; i < Rows; i++)
      for (int j = 0; j < Columns; j++)
        Elements[i][j] = 0;
  };

  // another constructor that copies the elements from a plain array
  Matrix (Ring _Elements[Rows][Columns])
  {
    for (int i = 0; i < Rows; i++)
      for (int j = 0; j < Columns; j++)
        Elements[i][j] = _Elements[i][j];
  };

  // operator to allow easy element access in two forms (const and non-const)
  Ring *operator [] (int i) { return Elements[i]; }
  const Ring *operator [] (int i) const { return Elements[i]; }

  // simple print routine
  void Print (void)
  {
    for (int i = 0; i < Rows; i++)
      {
        for (int j = 0; j < Columns; j++)
          std::cout << std::setw (20) << Elements[i][j];
        std::cout << std::endl;
      }
  }

  #ifdef METHOD_OPERATOR
  // multiplication as a method operator
  template <int Columns2>
  Matrix <Rows, Columns2, Ring> operator * (const Matrix <Columns, Columns2, Ring> &m2)
  {
    Matrix <Rows, Columns2, Ring> m;
    for (int i = 0; i < Rows; i++)
      for (int j = 0; j < Columns; j++)
        for (int k = 0; k < Columns2; k++)
          m[i][k] = m[i][k] + Elements[i][j] * m2[j][k];
    return m;
  }
  #endif
};

#ifndef METHOD_OPERATOR
// multiplication as a stand-alone operator
template <int Rows1, int Columns1, int Columns2, class Ring>
Matrix <Rows1, Columns2, Ring> operator * (const Matrix <Rows1,    Columns1, Ring> &m1,
                                           const Matrix <Columns1, Columns2, Ring> &m2)
{
  Matrix <Rows1, Columns2, Ring> m;
  for (int i = 0; i < Rows1; i++)
    for (int j = 0; j < Columns1; j++)
      for (int k = 0; k < Columns2; k++)
        m[i][k] = m[i][k] + m1[i][j] * m2[j][k];
  return m;
}
#endif

// Demo
int main ()
{
  // declare 3 matrices
  Matrix <3, 2> m1 ((double [3][2]) { { 1, 0 }, { 0, 1 }, { 0.5, 0.5 } });
  Matrix <2>    m2 ((double [2][2]) { { 0, 1 }, { 1, 0 } });
  Matrix <2, 3> m3 ((double [2][3]) { { 10, 20, 30 }, { 0.4, 0.5, 0.6 } });

  // multiply them and print the result
  // Note: The result type is a 3x3 matrix, which is never explicitly
  //       declared, but inferred by the compiler.
  (m1 * m2 * m3).Print ();

  // this would give a compile time error because the dimensions don't fit
  // m2 * m1;

  // how to declare an integer matrix (not used in the demo)
  Matrix <3, 2, int> foo;

  return 0;
}