C++
Шаблоны выписки
Поиск…
Базовые шаблоны выражений на элементарных алгебраических выражениях
Введение и мотивация
Шаблоны экспрессии (обозначаемые как ET в следующем) - это мощная технология метапрограммирования шаблонов, используемая для ускорения вычислений иногда довольно дорогостоящих выражений. Он широко используется в разных областях, например, в реализации библиотек линейной алгебры.
В этом примере рассмотрим контекст линейных алгебраических вычислений. Более конкретно, вычисления, включающие только элементарные операции . Такие вычисления являются наиболее распространенными приложениями ET , и они служат хорошим представлением о том, как ET работают внутри страны.
Давайте посмотрим на мотивирующий пример. Рассмотрим вычисление выражения:
Vector vec_1, vec_2, vec_3;
// Initializing vec_1, vec_2 and vec_3.
Vector result = vec_1 + vec_2*vec_3;
Здесь для простоты я предполагаю, что класс Vector и операция + (вектор плюс: элемент-плюс плюс операция) и операция * (здесь означает векторное скалярное произведение: также элементная операция), как правильно реализованы, так и как они должны быть, математически.
В обычной реализации без использования ET (или других аналогичных методов) для получения конечного result имеет место не менее пяти конструкций векторов Vector :
- Три экземпляра, соответствующие
vec_1,vec_2иvec_3. - Временной
Vectorэкземпляр_tmp, представляющий результат_tmp = vec_2*vec_3;, - Наконец, при правильном использовании оптимизации возвращаемого значения построение конечного
resultвresult = vec_1 + _tmp;,
Реализация с использованием ET может исключить создание временного Vector _tmp в 2, оставив при этом только четыре конструкции экземпляров Vector . Более интересно, рассмотрим следующее более сложное выражение:
Vector result = vec_1 + (vec_2*vec3 + vec_1)*(vec_2 + vec_3*vec_1);
Также будет четыре конструкции экземпляров Vector : vec_1, vec_2, vec_3 и result . Другими словами, в этом примере, где задействуются только элементарные операции , гарантируется, что из промежуточных вычислений не создаются временные объекты .
Как работают ET
В принципе, ET для любых алгебраических вычислений состоят из двух строительных блоков:
- Чистые алгебраические выражения ( PAE ): они являются прокси-абстракциями алгебраических выражений. Чистая алгебраика не выполняет фактических вычислений, это просто абстракции / моделирование рабочего потока вычислений. PAE может быть моделью либо ввода, либо вывода любых алгебраических выражений . Случаи PAE часто считаются дешевыми для копирования.
- Ленивые оценки : реализация реальных вычислений. В следующем примере мы увидим, что для выражений, содержащих только элементарные операции, ленивые оценки могут реализовывать фактические вычисления внутри операции индексированного доступа в конечном результате, создавая таким образом схему оценки по требованию: вычисление не выполняется только в том случае, если конечный результат был получен / запрошен.
Итак, конкретно, как мы внедряем ET в этом примере? Пойдем теперь через него.
Рассмотрим всегда следующий фрагмент кода:
Vector vec_1, vec_2, vec_3;
// Initializing vec_1, vec_2 and vec_3.
Vector result = vec_1 + vec_2*vec_3;
Выражение для вычисления результата может быть далее разложено на два подвыражения:
- Векторное выражение (обозначенное как plus_expr )
- Векторное выражение внутреннего произведения (обозначается как innerprod_expr ).
Что делают ГЭТ :
Вместо того, чтобы сразу вычислить каждое подвыражение , ETs сначала моделируют все выражение, используя графическую структуру. Каждый узел в графе представляет собой PAE . Конечное соединение узлов представляет собой фактический поток вычислений. Итак, для приведенного выше выражения мы получаем следующий график:
result = plus_expr( vec_1, innerprod_expr(vec_2, vec_3) ) / \ / \ / \ / innerprod_expr( vec_2, vec_3 ) / / \ / / \ / / \ vec_1 vec_2 vec_3
Окончательное вычисление реализуется путем просмотра иерархии графов : поскольку здесь мы имеем дело только с элементарными операциями, вычисление каждого проиндексированного значения в
resultможет быть выполнено независимо : окончательная оценкаresultможет быть лениво перенесена на элемент- мудрая оценка каждого элементаresult. Другими словами, поскольку вычисление элементаresultelem_resможет быть выражено с использованием соответствующих элементов вvec_1(elem_1),vec_2(elem_2) иvec_3(elem_3) как:elem_res = elem_1 + elem_2*elem_3;
поэтому нет необходимости создавать временный Vector для хранения результата промежуточного внутреннего продукта: весь расчет для одного элемента может быть выполнен в целом и закодирован внутри операции с индексированным доступом .
Вот пример кода в действии.
Файл vec.hh: wrapper для std :: vector, используемый для отображения журнала при вызове конструкции.
#ifndef EXPR_VEC
# define EXPR_VEC
# include <vector>
# include <cassert>
# include <utility>
# include <iostream>
# include <algorithm>
# include <functional>
///
/// This is a wrapper for std::vector. It's only purpose is to print out a log when a
/// vector constructions in called.
/// It wraps the indexed access operator [] and the size() method, which are
/// important for later ETs implementation.
///
// std::vector wrapper.
template<typename ScalarType> class Vector
{
public:
explicit Vector() { std::cout << "ctor called.\n"; };
explicit Vector(int size): _vec(size) { std::cout << "ctor called.\n"; };
explicit Vector(const std::vector<ScalarType> &vec): _vec(vec)
{ std::cout << "ctor called.\n"; };
Vector(const Vector<ScalarType> & vec): _vec{vec()}
{ std::cout << "copy ctor called.\n"; };
Vector(Vector<ScalarType> && vec): _vec(std::move(vec()))
{ std::cout << "move ctor called.\n"; };
Vector<ScalarType> & operator=(const Vector<ScalarType> &) = default;
Vector<ScalarType> & operator=(Vector<ScalarType> &&) = default;
decltype(auto) operator[](int indx) { return _vec[indx]; }
decltype(auto) operator[](int indx) const { return _vec[indx]; }
decltype(auto) operator()() & { return (_vec); };
decltype(auto) operator()() const & { return (_vec); };
Vector<ScalarType> && operator()() && { return std::move(*this); }
int size() const { return _vec.size(); }
private:
std::vector<ScalarType> _vec;
};
///
/// These are conventional overloads of operator + (the vector plus operation)
/// and operator * (the vector inner product operation) without using the expression
/// templates. They are later used for bench-marking purpose.
///
// + (vector plus) operator.
template<typename ScalarType>
auto operator+(const Vector<ScalarType> &lhs, const Vector<ScalarType> &rhs)
{
assert(lhs().size() == rhs().size() &&
"error: ops plus -> lhs and rhs size mismatch.");
std::vector<ScalarType> _vec;
_vec.resize(lhs().size());
std::transform(std::cbegin(lhs()), std::cend(lhs()),
std::cbegin(rhs()), std::begin(_vec),
std::plus<>());
return Vector<ScalarType>(std::move(_vec));
}
// * (vector inner product) operator.
template<typename ScalarType>
auto operator*(const Vector<ScalarType> &lhs, const Vector<ScalarType> &rhs)
{
assert(lhs().size() == rhs().size() &&
"error: ops multiplies -> lhs and rhs size mismatch.");
std::vector<ScalarType> _vec;
_vec.resize(lhs().size());
std::transform(std::cbegin(lhs()), std::cend(lhs()),
std::cbegin(rhs()), std::begin(_vec),
std::multiplies<>());
return Vector<ScalarType>(std::move(_vec));
}
#endif //!EXPR_VEC
Файл expr.hh: реализация шаблонов выражений для элементарных операций (векторный плюс и векторный скалярный продукт)
Давайте разложим его на разделы.
- Раздел 1 реализует базовый класс для всех выражений. В нем используется Curiously Recurring Template Pattern ( CRTP ).
- Раздел 2 реализует первый PAE : терминал , который является просто оболочкой (константной ссылкой) структуры входных данных, содержащей реальное входное значение для вычисления.
- Раздел 3 реализует второй PAE : binary_operation , который является шаблоном класса, который впоследствии используется для vector_plus и vector_innerprod. Он параметризуется типом операции , левым PAE и правым PAE . Фактическое вычисление кодируется в операторе индексированного доступа.
- Раздел 4 определяет операции vector_plus и vector_innerprod как элементную операцию . Он также перегружает оператор + и * для PAE s: так что эти две операции также возвращают PAE .
#ifndef EXPR_EXPR
# define EXPR_EXPR
/// Fwd declaration.
template<typename> class Vector;
namespace expr
{
/// -----------------------------------------
///
/// Section 1.
///
/// The first section is a base class template for all kinds of expression. It
/// employs the Curiously Recurring Template Pattern, which enables its instantiation
/// to any kind of expression structure inheriting from it.
///
/// -----------------------------------------
/// Base class for all expressions.
template<typename Expr> class expr_base
{
public:
const Expr& self() const { return static_cast<const Expr&>(*this); }
Expr& self() { return static_cast<Expr&>(*this); }
protected:
explicit expr_base() {};
int size() const { return self().size_impl(); }
auto operator[](int indx) const { return self().at_impl(indx); }
auto operator()() const { return self()(); };
};
/// -----------------------------------------
///
/// The following section 2 & 3 are abstractions of pure algebraic expressions (PAE).
/// Any PAE can be converted to a real object instance using operator(): it is in
/// this conversion process, where the real computations are done.
///
/// Section 2. Terminal
///
/// A terminal is an abstraction wrapping a const reference to the Vector data
/// structure. It inherits from expr_base, therefore providing a unified interface
/// wrapping a Vector into a PAE.
///
/// It provides the size() method, indexed access through at_impl() and a conversion
/// to referenced object through () operator.
///
/// It might no be necessary for user defined data structures to have a terminal
/// wrapper, since user defined structure can inherit expr_base, therefore eliminates
/// the need to provide such terminal wrapper.
///
/// -----------------------------------------
/// Generic wrapper for underlying data structure.
template<typename DataType> class terminal: expr_base<terminal<DataType>>
{
public:
using base_type = expr_base<terminal<DataType>>;
using base_type::size;
using base_type::operator[];
friend base_type;
explicit terminal(const DataType &val): _val(val) {}
int size_impl() const { return _val.size(); };
auto at_impl(int indx) const { return _val[indx]; };
decltype(auto) operator()() const { return (_val); }
private:
const DataType &_val;
};
/// -----------------------------------------
///
/// Section 3. Binary operation expression.
///
/// This is a PAE abstraction of any binary expression. Similarly it inherits from
/// expr_base.
///
/// It provides the size() method, indexed access through at_impl() and a conversion
/// to referenced object through () operator. Each call to the at_impl() method is
/// a element wise computation.
///
/// -----------------------------------------
/// Generic wrapper for binary operations (that are element-wise).
template<typename Ops, typename lExpr, typename rExpr>
class binary_ops: public expr_base<binary_ops<Ops,lExpr,rExpr>>
{
public:
using base_type = expr_base<binary_ops<Ops,lExpr,rExpr>>;
using base_type::size;
using base_type::operator[];
friend base_type;
explicit binary_ops(const Ops &ops, const lExpr &lxpr, const rExpr &rxpr)
: _ops(ops), _lxpr(lxpr), _rxpr(rxpr) {};
int size_impl() const { return _lxpr.size(); };
/// This does the element-wise computation for index indx.
auto at_impl(int indx) const { return _ops(_lxpr[indx], _rxpr[indx]); };
/// Conversion from arbitrary expr to concrete data type. It evaluates
/// element-wise computations for all indices.
template<typename DataType> operator DataType()
{
DataType _vec(size());
for(int _ind = 0; _ind < _vec.size(); ++_ind)
_vec[_ind] = (*this)[_ind];
return _vec;
}
private: /// Ops and expr are assumed cheap to copy.
Ops _ops;
lExpr _lxpr;
rExpr _rxpr;
};
/// -----------------------------------------
/// Section 4.
///
/// The following two structs defines algebraic operations on PAEs: here only vector
/// plus and vector inner product are implemented.
///
/// First, some element-wise operations are defined : in other words, vec_plus and
/// vec_prod acts on elements in Vectors, but not whole Vectors.
///
/// Then, operator + & * are overloaded on PAEs, such that: + & * operations on PAEs
/// also return PAEs.
///
/// -----------------------------------------
/// Element-wise plus operation.
struct vec_plus_t
{
constexpr explicit vec_plus_t() = default;
template<typename LType, typename RType>
auto operator()(const LType &lhs, const RType &rhs) const
{ return lhs+rhs; }
};
/// Element-wise inner product operation.
struct vec_prod_t
{
constexpr explicit vec_prod_t() = default;
template<typename LType, typename RType>
auto operator()(const LType &lhs, const RType &rhs) const
{ return lhs*rhs; }
};
/// Constant plus and inner product operator objects.
constexpr vec_plus_t vec_plus{};
constexpr vec_prod_t vec_prod{};
/// Plus operator overload on expressions: return binary expression.
template<typename lExpr, typename rExpr>
auto operator+(const lExpr &lhs, const rExpr &rhs)
{ return binary_ops<vec_plus_t,lExpr,rExpr>(vec_plus,lhs,rhs); }
/// Inner prod operator overload on expressions: return binary expression.
template<typename lExpr, typename rExpr>
auto operator*(const lExpr &lhs, const rExpr &rhs)
{ return binary_ops<vec_prod_t,lExpr,rExpr>(vec_prod,lhs,rhs); }
} //!expr
#endif //!EXPR_EXPR
Файл main.cc: файл test src
# include <chrono>
# include <iomanip>
# include <iostream>
# include "vec.hh"
# include "expr.hh"
# include "boost/core/demangle.hpp"
int main()
{
using dtype = float;
constexpr int size = 5e7;
std::vector<dtype> _vec1(size);
std::vector<dtype> _vec2(size);
std::vector<dtype> _vec3(size);
// ... Initialize vectors' contents.
Vector<dtype> vec1(std::move(_vec1));
Vector<dtype> vec2(std::move(_vec2));
Vector<dtype> vec3(std::move(_vec3));
unsigned long start_ms_no_ets =
std::chrono::duration_cast<std::chrono::milliseconds>
(std::chrono::system_clock::now().time_since_epoch()).count();
std::cout << "\nNo-ETs evaluation starts.\n";
Vector<dtype> result_no_ets = vec1 + (vec2*vec3);
unsigned long stop_ms_no_ets =
std::chrono::duration_cast<std::chrono::milliseconds>
(std::chrono::system_clock::now().time_since_epoch()).count();
std::cout << std::setprecision(6) << std::fixed
<< "No-ETs. Time eclapses: " << (stop_ms_no_ets-start_ms_no_ets)/1000.0
<< " s.\n" << std::endl;
unsigned long start_ms_ets =
std::chrono::duration_cast<std::chrono::milliseconds>
(std::chrono::system_clock::now().time_since_epoch()).count();
std::cout << "Evaluation using ETs starts.\n";
expr::terminal<Vector<dtype>> vec4(vec1);
expr::terminal<Vector<dtype>> vec5(vec2);
expr::terminal<Vector<dtype>> vec6(vec3);
Vector<dtype> result_ets = (vec4 + vec5*vec6);
unsigned long stop_ms_ets =
std::chrono::duration_cast<std::chrono::milliseconds>
(std::chrono::system_clock::now().time_since_epoch()).count();
std::cout << std::setprecision(6) << std::fixed
<< "With ETs. Time eclapses: " << (stop_ms_ets-start_ms_ets)/1000.0
<< " s.\n" << std::endl;
auto ets_ret_type = (vec4 + vec5*vec6);
std::cout << "\nETs result's type:\n";
std::cout << boost::core::demangle( typeid(decltype(ets_ret_type)).name() ) << '\n';
return 0;
}
Вот один из возможных результатов при компиляции с -O3 -std=c++14 с использованием GCC 5.3:
ctor called.
ctor called.
ctor called.
No-ETs evaluation starts.
ctor called.
ctor called.
No-ETs. Time eclapses: 0.571000 s.
Evaluation using ETs starts.
ctor called.
With ETs. Time eclapses: 0.164000 s.
ETs result's type:
expr::binary_ops<expr::vec_plus_t, expr::terminal<Vector<float> >, expr::binary_ops<expr::vec_prod_t, expr::terminal<Vector<float> >, expr::terminal<Vector<float> > > >
Наблюдения:
- Использование ETs в этом случае приводит к значительному повышению производительности (> 3x).
- Создание временного векторного объекта устраняется. Как и в случае ГЭ т е р вызывается только один раз.
- Boost :: demangle использовался для визуализации типа возврата ET перед преобразованием: он четко построил точно такой же график выражений, что и выше.
Оттяжки и оговорки
Очевидным недостатком ГЭ является кривая обучения, сложность реализации и сложность обслуживания кода. В приведенном выше примере, где рассматриваются только элементарные операции, реализация содержит уже огромное количество шаблонов, не говоря уже о реальном мире, где более сложные алгебраические выражения встречаются в каждом вычислении, а элементная независимость больше не выполняется (например, матричное умножение ), трудность будет экспоненциальной.
Еще одно предостережение от использования ET - это то, что они хорошо
autoключевым словомauto. Как упоминалось выше, PAE s по существу являются прокси-серверами: и прокси-серверы в основном не хорошо работают сauto. Рассмотрим следующий пример:auto result = ...; // Some expensive expression: // auto returns the expr graph, // NOT the computed value. for(auto i = 0; i < 100; ++i) ScalrType value = result* ... // Some other expensive computations using result.
Здесь в каждой итерации цикла for результат будет переоценен , так как граф выражения вместо вычисленного значения передается в цикл for.
Существующие библиотеки, внедряющие ET
- boost :: proto - это мощная библиотека, позволяющая вам определять свои собственные правила и грамматики для собственных выражений и выполнять их с помощью ET .
- Eigen - библиотека для линейной алгебры, которая эффективно использует различные алгебраические вычисления с использованием ET .
Основной пример, иллюстрирующий шаблоны выражений
Шаблон выражения - это метод оптимизации времени компиляции, используемый в основном в научных вычислениях. Основная цель - избегать ненужных временных рядов и оптимизировать вычисления циклов с использованием одного прохода (как правило, при выполнении операций с числовыми агрегатами). Шаблоны экспрессии были первоначально разработаны для того, чтобы обойти неэффективность наивной перегрузки оператора при реализации численных типов Array или Matrix . Эквивалентная терминология для шаблонов выражений была введена Бьярне Страуступом, который называет их «сплавленными операциями» в последней версии своей книги «Язык программирования C ++».
Прежде чем погрузиться в шаблоны выражений, вы должны понять, почему они вам нужны в первую очередь. Чтобы проиллюстрировать это, рассмотрим очень простой класс Matrix, приведенный ниже:
template <typename T, std::size_t COL, std::size_t ROW>
class Matrix {
public:
using value_type = T;
Matrix() : values(COL * ROW) {}
static size_t cols() { return COL; }
static size_t rows() { return ROW; }
const T& operator()(size_t x, size_t y) const { return values[y * COL + x]; }
T& operator()(size_t x, size_t y) { return values[y * COL + x]; }
private:
std::vector<T> values;
};
template <typename T, std::size_t COL, std::size_t ROW>
Matrix<T, COL, ROW>
operator+(const Matrix<T, COL, ROW>& lhs, const Matrix<T, COL, ROW>& rhs)
{
Matrix<T, COL, ROW> result;
for (size_t y = 0; y != lhs.rows(); ++y) {
for (size_t x = 0; x != lhs.cols(); ++x) {
result(x, y) = lhs(x, y) + rhs(x, y);
}
}
return result;
}
Учитывая предыдущее определение класса, вы можете теперь написать выражения Matrix, такие как:
const std::size_t cols = 2000;
const std::size_t rows = 1000;
Matrix<double, cols, rows> a, b, c;
// initialize a, b & c
for (std::size_t y = 0; y != rows; ++y) {
for (std::size_t x = 0; x != cols; ++x) {
a(x, y) = 1.0;
b(x, y) = 2.0;
c(x, y) = 3.0;
}
}
Matrix<double, cols, rows> d = a + b + c; // d(x, y) = 6
Как показано выше, возможность перегрузить operator+() дает вам обозначение, которое имитирует естественную математическую нотацию для матриц.
К сожалению, предыдущая реализация также очень неэффективна по сравнению с эквивалентной версией «ручной работы».
Чтобы понять, почему, вы должны учитывать, что происходит, когда вы пишете выражение, такое как Matrix d = a + b + c . Это фактически расширяется до ((a + b) + c) или operator+(operator+(a, b), c) . Другими словами, цикл внутри operator+() выполняется дважды, тогда как его можно было бы легко выполнить за один проход. Это также приводит к созданию 2 временных времен, что еще больше ухудшает производительность. По сути, добавив гибкость в использование нотации, близкой к ее математическому аналогу, вы также сделали класс Matrix крайне неэффективным.
Например, без перегрузки оператора вы могли бы реализовать гораздо более эффективное суммирование Matrix с использованием одного прохода:
template<typename T, std::size_t COL, std::size_t ROW>
Matrix<T, COL, ROW> add3(const Matrix<T, COL, ROW>& a,
const Matrix<T, COL, ROW>& b,
const Matrix<T, COL, ROW>& c)
{
Matrix<T, COL, ROW> result;
for (size_t y = 0; y != ROW; ++y) {
for (size_t x = 0; x != COL; ++x) {
result(x, y) = a(x, y) + b(x, y) + c(x, y);
}
}
return result;
}
Однако предыдущий пример имеет свои недостатки, поскольку он создает гораздо более запутанный интерфейс для класса Matrix (вам нужно будет рассмотреть такие методы, как Matrix::add2() , Matrix::AddMultiply() и т. Д.).
Вместо этого сделаем шаг назад и посмотрим, как мы можем адаптировать перегрузку оператора для более эффективного выполнения
Проблема связана с тем, что выражение Matrix d = a + b + c оценивается слишком «нетерпеливо», прежде чем у вас появилась возможность построить все дерево выражений. Другими словами, то, что вы действительно хотите достичь, - это оценить a + b + c за один проход и только один раз, когда вам действительно нужно назначить полученное выражение d .
Это основная идея шаблонов выражений: вместо того, чтобы operator+() сразу оценил результат добавления двух экземпляров Matrix, он вернет «шаблон выражения» для будущей оценки после того, как будет создано все дерево выражений.
Например, вот возможная реализация шаблона выражения, соответствующего суммированию двух типов:
template <typename LHS, typename RHS>
class MatrixSum
{
public:
using value_type = typename LHS::value_type;
MatrixSum(const LHS& lhs, const RHS& rhs) : rhs(rhs), lhs(lhs) {}
value_type operator() (int x, int y) const {
return lhs(x, y) + rhs(x, y);
}
private:
const LHS& lhs;
const RHS& rhs;
};
И вот обновленная версия operator+()
template <typename LHS, typename RHS>
MatrixSum<LHS, RHS> operator+(const LHS& lhs, const LHS& rhs) {
return MatrixSum<LHS, RHS>(lhs, rhs);
}
Как вы можете видеть, operator+() больше не возвращает «нетерпеливую оценку» результата добавления 2 экземпляров Matrix (это будет другой экземпляр Matrix), а вместо этого шаблон выражения, представляющий операцию добавления. Самый важный момент, который следует иметь в виду, заключается в том, что выражение еще не оценено. Он просто содержит ссылки на его операнды.
Фактически, ничто не мешает вам создать экземпляр MatrixSum<> выражения MatrixSum<> следующим образом:
MatrixSum<Matrix<double>, Matrix<double> > SumAB(a, b);
Вы можете, однако, на более позднем этапе, когда вам действительно нужен результат суммирования, оцените выражение d = a + b следующим образом:
for (std::size_t y = 0; y != a.rows(); ++y) {
for (std::size_t x = 0; x != a.cols(); ++x) {
d(x, y) = SumAB(x, y);
}
}
Как вы можете видеть, еще одно преимущество использования шаблона выражения заключается в том, что вам в основном удалось оценить сумму a и b и назначить его d за один проход.
Кроме того, ничто не мешает вам комбинировать несколько шаблонов выражений. Например, a + b + c приведет к следующему шаблону выражения:
MatrixSum<MatrixSum<Matrix<double>, Matrix<double> >, Matrix<double> > SumABC(SumAB, c);
И здесь вы можете оценить окончательный результат с помощью одного прохода:
for (std::size_t y = 0; y != a.rows(); ++y) {
for (std::size_t x = 0; x != a.cols(); ++x) {
d(x, y) = SumABC(x, y);
}
}
Наконец, последний фрагмент головоломки состоит в том, чтобы фактически подключить ваш шаблон выражения к классу Matrix . Это достигается, главным образом, за счет реализации Matrix::operator=() , которая принимает шаблон выражения в качестве аргумента и оценивает его за один проход, как вы делали «вручную» до:
template <typename T, std::size_t COL, std::size_t ROW>
class Matrix {
public:
using value_type = T;
Matrix() : values(COL * ROW) {}
static size_t cols() { return COL; }
static size_t rows() { return ROW; }
const T& operator()(size_t x, size_t y) const { return values[y * COL + x]; }
T& operator()(size_t x, size_t y) { return values[y * COL + x]; }
template <typename E>
Matrix<T, COL, ROW>& operator=(const E& expression) {
for (std::size_t y = 0; y != rows(); ++y) {
for (std::size_t x = 0; x != cols(); ++x) {
values[y * COL + x] = expression(x, y);
}
}
return *this;
}
private:
std::vector<T> values;
};