#ifndef MATRIX_HPP
#define MATRIX_HPP

#include <cassert>
#include <cstdlib>
#include <memory>
#include <iomanip>

template<typename T>
class Matrix {
   public:
      using ElementType = T;

      Matrix() : m(nullptr), offset_start(0), offset_nextcol(0),
	    offset_nextrow(0), rows(0), columns(0) {
      }
      /* keep default copy constructor */
      Matrix(std::size_t rows, std::size_t columns) :
	    m(new T[rows * columns], [](T *p) { delete[] p; }),
	    offset_start(0), offset_nextcol(1), offset_nextrow(columns),
	    rows(rows), columns(columns) {
      }
      template<typename Initializer>
      Matrix(std::size_t rows, std::size_t columns,
	       Initializer initializer) :
	    m(new T[rows * columns], [](T *p) { delete[] p; }),
	    offset_start(0), offset_nextcol(1), offset_nextrow(columns),
	    rows(rows), columns(columns) {
	 for (std::size_t i = 0; i < rows; ++i) {
	    for (std::size_t j = 0; j < columns; ++j) {
	       m.get()[i * columns + j] = initializer(i, j);
	    }
	 }
      }
      Matrix(std::shared_ptr<T> base, std::size_t rows,
	    std::size_t columns) :
	 m(base), offset_start(0), offset_nextcol(1),
	 offset_nextrow(columns), rows(rows), columns(columns) {
      }
      Matrix(std::shared_ptr<T> base, std::size_t rows,
	    std::size_t columns, std::ptrdiff_t offset_start) :
	 m(base), offset_start(offset_start), offset_nextcol(1),
	 offset_nextrow(columns),
	 rows(rows), columns(columns) {
      }
      Matrix(std::shared_ptr<T> base, std::size_t rows,
	    std::size_t columns, std::ptrdiff_t offset_start,
	    std::ptrdiff_t offset_nextrow) :
	 m(base), offset_start(offset_start), offset_nextcol(1),
	 offset_nextrow(offset_nextrow),
	 rows(rows), columns(columns) {
      }
      Matrix(std::shared_ptr<T> base, std::size_t rows,
	    std::size_t columns, std::ptrdiff_t offset_start,
	    std::ptrdiff_t offset_nextrow, std::ptrdiff_t offset_nextcol) :
	 m(base), offset_start(offset_start), offset_nextcol(offset_nextcol),
	 offset_nextrow(offset_nextrow), rows(rows), columns(columns) {
      }

      Matrix get_block(std::size_t start_row, std::size_t start_column,
	       std::size_t block_rows, unsigned block_columns) {
	 assert(start_row + block_rows <= rows);
	 assert(start_column + block_columns <= columns);
	 return Matrix(m, block_rows, block_columns,
	    offset_start + start_row * offset_nextrow +
	       start_column * offset_nextcol,
	    offset_nextrow, offset_nextcol);
      }

      std::size_t get_rows() const {
	 return rows;
      }
      std::size_t get_columns() const {
	 return columns;
      }
      double& operator()(std::size_t i, std::size_t j) {
	 assert(i < rows && j < columns);
	 return *(m.get()
	    + offset_start + i * offset_nextrow + j * offset_nextcol);
      }
      const double& operator()(std::size_t i, std::size_t j) const {
	 assert(i < rows && j < columns);
	 return *(m.get()
	    + offset_start + i * offset_nextrow + j * offset_nextcol);
      }
   private:
      std::shared_ptr<T> m;
      std::ptrdiff_t offset_start;
      std::ptrdiff_t offset_nextcol;
      std::ptrdiff_t offset_nextrow;
      std::size_t rows;
      std::size_t columns;
};

template<typename T>
std::ostream& operator<<(std::ostream& out, const Matrix<T>& matrix) {
   for (std::size_t i = 0; i < matrix.get_rows(); ++i) {
      for (std::size_t j = 0; j < matrix.get_columns(); ++j) {
	 out << " " << matrix(i, j);
      }
      out << std::endl;
   }
   return out;
}

#endif