C++
Plantillas de expresiones
Buscar..
Plantillas de expresiones básicas en expresiones algebraicas de elementos
Introducción y motivación.
Las plantillas de expresión (indicadas como ETs a continuación) son una poderosa técnica de metaprogramación de plantilla, utilizada para acelerar los cálculos de expresiones a veces bastante caras. Se usa ampliamente en diferentes dominios, por ejemplo, en la implementación de bibliotecas de álgebra lineal.
Para este ejemplo, considere el contexto de los cálculos algebraicos lineales. Más específicamente, los cálculos que involucran solo operaciones de elementos . Este tipo de cálculos son las aplicaciones más básicas de los ET y sirven como una buena introducción a cómo funcionan los ET internamente.
Veamos un ejemplo motivador. Considere el cálculo de la expresión:
Vector vec_1, vec_2, vec_3;
// Initializing vec_1, vec_2 and vec_3.
Vector result = vec_1 + vec_2*vec_3;
En aras de la simplicidad, asumiré que la clase Vector y operation + (vector plus: element-wise plus operation) y operation * (aquí significa vector producto interno: también element-operation operation) se implementan correctamente, como cómo deberían ser, matemáticamente.
En una implementación convencional sin utilizar ET (u otras técnicas similares), se realizan al menos cinco construcciones de instancias de Vector para obtener el result final:
- Tres instancias correspondientes a
vec_1,vec_2yvec_3. - Una instancia de
Vectortemporal_tmp, que representa el resultado de_tmp = vec_2*vec_3;. - Finalmente, con el uso adecuado de la optimización del valor de retorno , la construcción del
resultfinal en elresult = vec_1 + _tmp;.
La implementación con ET puede eliminar la creación de Vector _tmp temporal en 2, dejando solo cuatro construcciones de instancias de Vector . Más interesante, considere la siguiente expresión que es más compleja:
Vector result = vec_1 + (vec_2*vec3 + vec_1)*(vec_2 + vec_3*vec_1);
También habrá cuatro construcciones de instancias de Vector en total: vec_1, vec_2, vec_3 y result . En otras palabras, en este ejemplo, donde solo están involucradas operaciones de elementos , se garantiza que no se crearán objetos temporales a partir de cálculos intermedios .
¿Cómo funcionan los TE?
Básicamente, las ET para cualquier cálculo algebraico consisten en dos bloques de construcción:
- Expresiones algebraicas puras ( PAE ): son proxies / abstracciones de expresiones algebraicas. Un algebraico puro no realiza cálculos reales, son simplemente abstracciones / modelos del flujo de trabajo de cálculo. Un PAE puede ser un modelo de la entrada o la salida de cualquier expresión algebraica . Las instancias de PAE a menudo se consideran baratas para copiar.
- Evaluaciones perezosas : que son implementaciones de cálculos reales. En el siguiente ejemplo, veremos que para expresiones que solo involucran operaciones de elementos, las evaluaciones perezosas pueden implementar cálculos reales dentro de la operación de acceso indexado en el resultado final, creando así un esquema de evaluación bajo demanda: no se realiza un cálculo Solo si se accede / se solicita el resultado final.
Entonces, específicamente, ¿cómo implementamos ETs en este ejemplo? Vamos a caminar a través de él ahora.
Considere siempre el siguiente fragmento de código:
Vector vec_1, vec_2, vec_3;
// Initializing vec_1, vec_2 and vec_3.
Vector result = vec_1 + vec_2*vec_3;
La expresión para calcular el resultado se puede descomponer en dos subexpresiones:
- Un vector más expresión (denotado como plus_expr )
- Una expresión de producto interno del vector (denotada como innerprod_expr ).
Lo que hacen los ETs es lo siguiente:
En lugar de calcular de inmediato cada subexpresión , los ET primero modelan la expresión completa utilizando una estructura gráfica. Cada nodo en la gráfica representa un PAE . La conexión de borde de los nodos representa el flujo de cálculo real. Así que para la expresión anterior, obtenemos el siguiente gráfico:
result = plus_expr( vec_1, innerprod_expr(vec_2, vec_3) ) / \ / \ / \ / innerprod_expr( vec_2, vec_3 ) / / \ / / \ / / \ vec_1 vec_2 vec_3
El cálculo final se implementa mirando a través de la jerarquía del gráfico : ya que aquí se tratan solo operaciones de elementos , el cálculo de cada valor indexado en el
resultse puede realizar de forma independiente : la evaluación final delresultse puede posponer perezosamente a un elemento Evaluación sabia de cada elemento delresult. En otras palabras, dado que el cálculo de un elemento deresult,elem_res, se puede expresar utilizando los elementos correspondientes envec_1(elem_1),vec_2(elem_2) yvec_3(elem_3) como:elem_res = elem_1 + elem_2*elem_3;
por lo tanto, no es necesario crear un Vector temporal para almacenar el resultado del producto interno intermedio: todo el cómputo de un elemento se puede realizar en conjunto y se puede codificar dentro de la operación de acceso indexado .
Aquí están los códigos de ejemplo en acción.
Archivo vec.hh: wrapper para std :: vector, utilizado para mostrar el registro cuando se llama a una construcción.
#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
Archivo expr.hh: implementación de plantillas de expresión para operaciones de elementos (vector plus y vector inner product)
Vamos a dividirlo en secciones.
- La sección 1 implementa una clase base para todas las expresiones. Emplea el patrón de plantilla curiosamente recurrente ( CRTP ).
- La Sección 2 implementa el primer PAE : un terminal , que es solo una envoltura (referencia constante) de una estructura de datos de entrada que contiene un valor de entrada real para el cálculo.
- La Sección 3 implementa el segundo PAE : binary_operation , que es una plantilla de clase que se usará más tarde para vector_plus y vector_innerprod. Está parametrizado por el tipo de operación , el PAE del lado izquierdo y el PAE del lado derecho . El cálculo real se codifica en el operador de acceso indexado.
- La Sección 4 define las operaciones vector_plus y vector_innerprod como una operación de elementos . También sobrecarga al operador + y * para PAE s: de modo que estas dos operaciones también devuelven 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
Archivo main.cc: test src file
# 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;
}
Aquí hay una salida posible cuando se compila con -O3 -std=c++14 usando 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> > > >
Las observaciones son:
- El uso de ETs logra un aumento de rendimiento bastante significativo en este caso (> 3x).
- Se elimina la creación de objetos vectoriales temporales. Como en el caso de ETs , ctor es llamado solo una vez.
- Boost :: demangle se usó para visualizar el tipo de retorno de ET antes de la conversión: claramente construyó exactamente el mismo gráfico de expresión demostrado anteriormente.
Draw-backs y advertencias
Una desventaja obvia de las ET es la curva de aprendizaje, la complejidad de la implementación y la dificultad de mantenimiento del código. En el ejemplo anterior, donde solo se consideran las operaciones de elementos sabios, la implementación ya contiene una cantidad enorme de placas de preparación, y mucho menos en el mundo real, donde ocurren expresiones algebraicas más complejas en cada cálculo y la independencia de los elementos ya no es válida (por ejemplo, la multiplicación de matrices). ), la dificultad será exponencial.
Otra advertencia de usar ETs es que juegan bien con la palabra clave
auto. Como se mencionó anteriormente, los PAE son esencialmente proxies: y los proxies básicamente no funcionan bien con elauto. Considere el siguiente ejemplo: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.
Aquí, en cada iteración del bucle for, el resultado se volverá a evaluar , ya que el gráfico de expresión en lugar del valor calculado se pasa al bucle for.
Bibliotecas existentes implementando ETs
- boost :: proto es una biblioteca poderosa que te permite definir tus propias reglas y gramáticas para tus propias expresiones y ejecutar usando ETs .
- Eigen es una biblioteca para álgebra lineal que implementa varios cálculos algebraicos de manera eficiente utilizando ETs .
Un ejemplo básico que ilustra plantillas de expresiones.
Una plantilla de expresión es una técnica de optimización en tiempo de compilación utilizada principalmente en computación científica. Su propósito principal es evitar temporarios innecesarios y optimizar los cálculos de bucle usando una sola pasada (normalmente cuando se realizan operaciones en agregados numéricos). Plantillas de expresiones se diseñaron inicialmente con el fin de eludir las ineficiencias de la sobrecarga de operadores ingenuo al implementar numéricos Array o Matrix tipos. Bjarne Stroustrup introdujo una terminología equivalente para las plantillas de expresión, que las llama "operaciones fusionadas" en la última versión de su libro, "El lenguaje de programación C ++".
Antes de sumergirse realmente en las plantillas de expresión, debe comprender por qué las necesita en primer lugar. Para ilustrar esto, considere la muy simple clase de Matrix que se da a continuación:
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;
}
Dada la definición de clase anterior, ahora puede escribir expresiones de matriz como:
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
Como se ilustra arriba, ser capaz de sobrecargar el operator+() proporciona una notación que imita la notación matemática natural de las matrices.
Desafortunadamente, la implementación anterior también es altamente ineficiente en comparación con una versión equivalente "hecha a mano".
Para entender por qué, debe considerar lo que sucede cuando escribe una expresión como Matrix d = a + b + c . De hecho, esto se expande a ((a + b) + c) u operator+(operator+(a, b), c) . En otras palabras, el bucle dentro del operator+() se ejecuta dos veces, mientras que podría haberse realizado fácilmente en una sola pasada. Esto también resulta en la creación de 2 temporarios, lo que degrada aún más el rendimiento. En esencia, al agregar la flexibilidad para usar una notación cercana a su contraparte matemática, también ha hecho que la clase Matrix altamente ineficiente.
Por ejemplo, sin la sobrecarga del operador, podría implementar una suma de Matrix mucho más eficiente utilizando un solo paso:
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;
}
Sin embargo, el ejemplo anterior tiene sus propias desventajas porque crea una interfaz mucho más compleja para la clase Matrix (tendrías que considerar métodos como Matrix::add2() , Matrix::AddMultiply() etc.).
En su lugar, demos un paso atrás y veamos cómo podemos adaptar la sobrecarga del operador para que funcione de una manera más eficiente.
El problema se deriva del hecho de que la expresión Matrix d = a + b + c se evalúa demasiado "con entusiasmo" antes de que haya tenido la oportunidad de construir todo el árbol de expresiones. En otras palabras, lo que realmente desea lograr es evaluar a + b + c en una sola pasada y solo una vez que realmente necesite asignar la expresión resultante a d .
Esta es la idea central detrás de las plantillas de expresión: en lugar de que el operator+() evalúe inmediatamente el resultado de agregar dos instancias de Matrix, devolverá una "plantilla de expresión" para una futura evaluación una vez que se haya construido todo el árbol de expresiones.
Por ejemplo, aquí hay una posible implementación para una plantilla de expresión correspondiente a la suma de 2 tipos:
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;
};
Y aquí está la versión actualizada del operator+()
template <typename LHS, typename RHS>
MatrixSum<LHS, RHS> operator+(const LHS& lhs, const LHS& rhs) {
return MatrixSum<LHS, RHS>(lhs, rhs);
}
Como puede ver, el operator+() ya no devuelve una "evaluación impaciente" del resultado de agregar 2 instancias de Matrix (que sería otra instancia de Matrix), sino una plantilla de expresión que representa la operación de adición. El punto más importante a tener en cuenta es que la expresión aún no se ha evaluado. Simplemente contiene referencias a sus operandos.
De hecho, nada le impide crear una instancia de la plantilla de expresión MatrixSum<> siguiente manera:
MatrixSum<Matrix<double>, Matrix<double> > SumAB(a, b);
Sin embargo, en una etapa posterior, cuando realmente necesite el resultado de la suma, evalúe la expresión d = a + b siguiente manera:
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);
}
}
Como se puede ver, otro de los beneficios de utilizar una plantilla de expresión, es que básicamente ha logrado evaluar la suma de a y b y asignarlo a d en una sola pasada.
Además, nada le impide combinar varias plantillas de expresión. Por ejemplo, a + b + c daría como resultado la siguiente plantilla de expresión:
MatrixSum<MatrixSum<Matrix<double>, Matrix<double> >, Matrix<double> > SumABC(SumAB, c);
Y aquí nuevamente puedes evaluar el resultado final usando una sola pasada:
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);
}
}
Finalmente, la última pieza del rompecabezas es en realidad conectar su plantilla de expresión a la clase Matrix . Esto se logra esencialmente al proporcionar una implementación para Matrix::operator=() , que toma la plantilla de expresión como un argumento y la evalúa en una sola pasada, como lo hizo "manualmente" antes:
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;
};
