サーチ…
要素指向の代数式の基本式テンプレート
紹介とモチベーション
式テンプレート (以下では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
とoperation +(vector plus:element-wise plus operation)とoperation *(ここではvector inner product:element-wise operation)が正しく実装されていると仮定します。彼らがどのようにすべきか、数学的に。
ET (または他の同様の技術)を使用しない従来の実装では、 Vector
インスタンスの少なくとも5つの構成が最終result
を得るために行われます。
-
vec_1
、vec_2
、およびvec_3
対応する3つのインスタンス。 - 一時
Vector
インスタンス_tmp
の結果を表す、_tmp = vec_2*vec_3;
。 - 最後に戻り値の最適化を適切に使用すると、
result = vec_1 + _tmp;
最終result
の構築result = vec_1 + _tmp;
。
ETを使用した実装では、2 で一時 Vector _tmp
作成する必要がなくなるため、 Vector
インスタンスの構造が4つだけ残されます。より興味深いことに、より複雑な次の式を考えてみましょう。
Vector result = vec_1 + (vec_2*vec3 + vec_1)*(vec_2 + vec_3*vec_1);
vec_1, vec_2, vec_3
およびresult
合計でVector
インスタンスの構造も4つあります。言い換えると、この例では、要素単位の操作のみが含まれているため、中間計算から一時オブジェクトが作成されないことが保証されています 。
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;
結果を計算する式は、さらに2つの部分式に分解できます。
- ベクトルプラス式( plus_exprと表記)
- ベクトル内積式( innerprod_exprと表記)。
ETのやり方は次のとおりです:
各サブ式をすぐに計算する代わりに、 ETは最初にグラフィカルな構造を使って表現全体をモデル化します。グラフの各ノードは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
の各要素の賢明な評価 。つまり、result
elem_res
要素の計算は、vec_1
(elem_1
)、vec_2
(elem_2
)、vec_3
(elem_3
)の対応する要素を使って次のように表現できます。elem_res = elem_1 + elem_2*elem_3;
したがって、中間の内積の結果を格納するために一時的なVector
を作成する必要はありません。つまり、1つの要素の計算全体をすべて実行し、索引アクセス操作内でエンコードすることができます 。
実際のコード例を示します。
ファイルvec.hh: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を実装しています 。これは、実際の入力値を含む入力データ構造の単なるラッパー(const参照)です。
- セクション3は、2番目のPAE : binary_operationを実装します 。これは、後でvector_plusとvector_innerprodで使用されるクラステンプレートです。 操作のタイプ 、左側のPAE 、右側の PAEによってパラメータ化されています。実際の計算は、インデックス付きアクセス演算子でエンコードされます。
- セクション4では、vector_plusおよびvector_innerprod演算を要素ごとの演算として定義しています。また、 PAEの演算子+と*をオーバーロードします。この2つの演算も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:テスト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;
}
GCC 5.3を使用して-O3 -std=c++14
コンパイルすると、次のような出力が得られ-O3 -std=c++14
。
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> > > >
観察は次のとおりです。
- ETを使用すると、この場合はパフォーマンスが大幅に向上します (> 3倍)。
- 一時的なVectorオブジェクトの作成は不要です。 ETの場合と同様に、ctorは一度だけ呼び出されます。
- Boost :: demangleは、変換前のETの戻り値の型を視覚化するために使用されました。
ドローバックと警告
ETの明らかな欠点は、学習曲線、実装の複雑さ、コードメンテナンスの難しさです。要素単位の操作のみが考慮されている上記の例では、実装には、すべての計算でより複雑な代数式が発生し、要素の独立性がもはや保持されていない現実世界ではなく、すでに膨大な量の定型文が含まれています(たとえば、 )、難易度は指数関数になります。
ETを使用する際のもう一つの注意点は、
auto
キーワードでうまくいくことです。上で述べたように、 PAEは本質的にプロキシです:そして、プロキシは基本的に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
型を実装する際の単純な演算子のオーバーロードの非効率性を回避するために考案されました。 Bjarne Stroustrup氏は、最新のバージョン「The C ++ Programming Language」の「融合操作」と呼ぶ表現テンプレートの同等の用語を導入しました。
実際に表現テンプレートに入る前に、なぜそれらを必要とするのか理解しておく必要があります。これを説明するために、以下に示す非常に単純な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;
}
前のクラス定義を前提とすると、次のような行列式を書くことができるようになりました。
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回実行されますが、1回のパスで簡単に実行できました。この結果、2つの一時的な一時ファイルが作成され、パフォーマンスがさらに低下します。本質的には、数学的な対応に近い表記を使用する柔軟性を追加することで、 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
があまりにも「熱心に」評価されるという事実に起因しています。言い換えれば、実際に達成したいのは、1回のパスでa + b + c
を評価a + b + c
ことです。実際には、結果として得られる表現をd
に割り当てる必要がある場合に限ります。
これは、式テンプレートの背後にあるコアアイデアです。 operator+()
2つのMatrixインスタンスを追加した結果を直ちに評価するのではなく、式ツリー全体が構築された後、将来の評価のために「式テンプレート」を返します。
たとえば、次の2つのタイプの合計に対応する式テンプレートの実装が考えられます。
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+()
更新版があり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<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);
}
}
ご覧のように、式テンプレートを使用するもう1つの利点は、基本的にはa
とb
合計を評価し、1回のパスでd
に割り当てることです。
また、複数の式テンプレートを組み合わせることを止めるものはありません。たとえば、 a + b + c
は次の式テンプレートになります。
MatrixSum<MatrixSum<Matrix<double>, Matrix<double> >, Matrix<double> > SumABC(SumAB, c);
ここでもまた、1回のパスを使用して最終結果を評価することができます。
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=()
実装を提供することによって実現されています。これは、式テンプレートを引数としてとり、前に "手動で"行ったように、1回のパスで評価します。
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;
};