#include <valarray>
#include <initializer_list>
#include <iterator>
#include <stdexcept>

using std::size_t;
using std::ptrdiff_t;
typedef size_t index_t;
using std::begin;
using std::end;

template<typename T>
using nested_init_list = std::initializer_list<std::initializer_list<T>>;

#include <valarray>
#include <initializer_list>
#include <iterator>
#include <stdexcept>

using std::size_t;
using std::ptrdiff_t;
typedef size_t index_t;

template<typename T>
using nested_init_list = std::initializer_list<std::initializer_list<T>>;

Iterators

Since each matrix class has a different internal structure, I needed to write custom iterator classes. To save some common code between the const and non-const iterators, I made a private nested type in Matrix that acts as a base class for both. The iterators store a pointer to a const Matrix and a pair of indices (i, j). The dereference operator then delegates to Matrix::operator().

To be able to write common code for comparison and arithmetic operators in the base class, I needed to store the matrix pointer in it. I had to decide between making it const or non-const. In the end I decided to keep it const and do a const_cast to Matrix * in the iterator (non-const) class when dereferencing the iterator. If I had done it the other way around, the begin and end member functions couldn't be marked const–or at least they would have to const_cast to Matrix *.

class Matrix::base_iterator
        : public std::iterator<std::random_access_iterator_tag, double> {
public:
    base_iterator &operator++();
    base_iterator &operator--();
    base_iterator &operator+=(ptrdiff_t offset);
    base_iterator &operator-=(ptrdiff_t offset) { return *this += -offset; }
    ptrdiff_t operator-(const base_iterator &rhs) const;

    bool operator==(const base_iterator &rhs) const;
    bool operator!=(const base_iterator &rhs) const { return not(*this == rhs); }
    bool operator<(const base_iterator &rhs) const { return rhs - *this > 0; }
    bool operator>(const base_iterator &rhs) const { return *this - rhs > 0; }
    bool operator<=(const base_iterator &rhs) const { return not(*this > rhs); }
    bool operator>=(const base_iterator &rhs) const { return not(*this < rhs); }

protected:
    base_iterator(const Matrix *mat, index_t i, index_t j) : mat(mat), i(i), j(j) {}

    const Matrix *mat;
    index_t i, j;
};

template<typename MatrixIter>
MatrixIter operator+(const MatrixIter &lhs, ptrdiff_t rhs);

template<typename MatrixIter>
MatrixIter operator-(const MatrixIter &lhs, ptrdiff_t rhs);


class Matrix::const_iterator : public Matrix::base_iterator {
public:
    explicit
    const_iterator(const Matrix *mat, index_t i = 0, index_t j = 0) : base_iterator(mat, i, j) {}

    double operator*() const { return mat->operator()(i, j); };
};


class Matrix::iterator : public Matrix::base_iterator {
public:
    explicit
    iterator(Matrix *mat, index_t i = 0, index_t j = 0) : base_iterator(mat, i, j) {}

    double &operator*() { return const_cast<Matrix *>(mat)->operator()(i, j); };
};

Implementation

Matrix::base_iterator &Matrix::base_iterator::operator++() {
    if (++j == mat->size()) { ++i; j = 0; }
    return *this;
}

Matrix::base_iterator &Matrix::base_iterator::operator+=(ptrdiff_t offset) {
    size_t n = mat->size();
    i += (j += offset)/n;
    j %= n;
    return *this;
}

Matrix::base_iterator &Matrix::base_iterator::operator--() {
    if (j == 0) { --i; j = mat->size() - 1; }
    else --j;
    return *this;
}

ptrdiff_t Matrix::base_iterator::operator-(const Matrix::base_iterator &rhs) const {
    return mat->flattenIndex(i, j) - mat->flattenIndex(rhs.i, rhs.j);
}

bool Matrix::base_iterator::operator==(const Matrix::base_iterator &rhs) const {
    return mat == rhs.mat and i == rhs.i and j == rhs.j;
}

Finally, the Matrix::begin and Matrix::end member functions are just:

Matrix::const_iterator Matrix::begin() const {
    return const_iterator(this, 0, 0);
}

Matrix::const_iterator Matrix::end() const {
    return const_iterator(this, _n, 0);
}

Matrix::iterator Matrix::begin() {
    return iterator(this, 0, 0);
}

Matrix::iterator Matrix::end() {
    return iterator(this, _n, 0);
}