Metaprogramming for math library
Jiří Chmelík, Marek Trtík
PA199

‹#›
uMetaprogram is a code whose compilation generates a new code whose functionality we really want.
u
uFirst metaprogram:
uErwin Unruh, 1994.
uDuring compilation of his program the compiler computed prime numbers.
u
uOur first metaprogram will be much simpler…
Metaprogramming

‹#›
Our first metaprogram


‹#›
uLet’s consider this implementation of the dot product:
u float dot_product(int dim, float* u, float* v) {
u float result = 0.0f;
u for (int i=0; i<dim; ++i)
u result += u[i] * v[i];
u return result;
u }
u float u[3] = {1,2,3}; // A 3D vector
u float v[3] = {4,5,6}; // A 3D vector
u std::cout << dot_product(3, u, v);  // Print their dot product.
uFaster code woud be: std::cout << u[0]*v[0] + u[1]*v[1] + u[2]*v[2];
Unrolling loops: Motivation

‹#›
utemplate<int DIM> struct DotProduct { // Recursion rule
u static inline float result(float* u, float* v) {
u return (*u) * (*v) + DotProduct<DIM-1>::result(u+1, v+1);
u }
u};
u
utemplate<> struct DotProduct<1> { // Base case
u static inline float result(float* u, float* v) {
u return (*u) * (*v);
u }
u};
Unrolling loops: Solution

‹#›
u template<int DIM> inline float dot_product(float* u, float* v) {
u return DotProduct<DIM>::result(u, v);
u }
uThe compiler then converts the call
u dot_product<3>(u,v)
u into code:
uDotProduct<3>::result(u,v);
u(*u) * (*v) + DotProduct<2>::result(u+1,v+1);
u(*u) * (*v) + ((*(u+1)) * (*(v+1)) + DotProduct<1>::result(u+2,v+2));
u(*u) * (*v) + ((*(u+1)) * (*(v+1)) + (*(u+2)) * (*(v+2)));
Unrolling loops: Solution

‹#›
uWhen designing a math library, it is desirable to provide this syntax:
u Vector u(1000), v(1000); // Define vectors
u … // Initialise them
u u = 1.2f * u + u * v; // Perform operations.
u
uIn C++ we can achieve this syntax via operator overloading.
u
uQuestion: Could there be a performance issue with such approach?
u
uLet’s start with a naïve solution…
Expression templates: Motivation

‹#›
uclass SVector {  // “Simple” Vector of any dimension.
u int size_;
u float* data_;
u void copy(SVector const& o) { for (int i=0; i<size(); ++i) data_[i] = o.data_[i]; }
upublic:
u explicit SVector(int size) : size_(size), data_(new float[size_]) {}
u ~SVector() { delete [] data_; }
u SVector(SVector const& o) : size_(o.size()), data_(new float[size_]) { copy(o); }
u SVector& operator=(SVector const& o) { copy(o); return *this; }
u int size() const { return size_; }
u float operator[](int const i) const { return data_[i]; }
u float& operator[](int const i) { return data_[i]; }
u};
Expression templates: Motivation

‹#›
u SVector operator+(SVector const& u, SVector const& v) {
u SVector result(u.size()); // A temporary for result computation.
u for (int i=0; i<u.size(); ++i)
u result[i] = u[i] + v[i];
u return result;
u }
u
u SVector operator*(SVector const& u, SVector const& v) {
u SVector result(u.size()); // A temporary for result computation.
u for (int i=0; i<u.size(); ++i)
u result[i] = u[i] * v[i];
u return result;
u }
Expression templates: Motivation

‹#›
u SVector operator*(float const a, SVector const& u) {
u SVector result(u.size()); // A temporary for result computation.
u for (int i=0; i<u.size(); ++i)
u result[i] = a * u[i];
u return result;
u }
Expression templates: Motivation

‹#›
uLet’s now use the code:
u SVector u(1000), v(1000);
u u = 1.2f * u + u * v;
uHow efficient is this code?
utmp1 = 1.2f * u;     // tmp1 == “result” variable in operator*.
uoperator* creates “result” (allocation of “data_”). Then 1000 iterations in the loop.
utmp2 = u * v;     // tmp2 == “result” variable in operator*.
uoperator* creates “result” (allocation of “data_”). Then 1000 iterations in the loop.
utmp3 = tmp1 + tmp2; // tmp3 == “result” variable in operator+.
uoperator+ creates “result” (allocation of “data_”). Then 1000 iterations in the loop.
uu = tmp3;
uoperator= performs 1000 iterations in the loop.
uThen 3 deallocations of arrays “data_” in tmp1, tmp2, and tmp3.
Expression templates: Motivation

‹#›
uWe can certainly do better (optimal code):
u for (int i=0; i<u.size(); ++i)
u u[i] = 1.2f * u[i] + u[i] * v[i]; // Temporaries are CPU registers!
uHow efficient is this code?
uNO RAM TEMPORARIES! (No allocations, no redundant RAM loads/stores).
uJust 1000 iterations in the loop.
u
uCan we achieve this efficiency while preserving the syntax?
uYES – using EXPRESSION TEMPLATES.
Expression templates: Motivation

‹#›
uWe let the compiler to transform the syntax tree of the expression      “u = 1.2f * u + u * v” by
hierarchy of simple template classes:
Expression templates: Idea
1.2f
u
u
v
*
*
+
SVector
SVector
SVector
Scalar
Mul<Scalar,SVector>
Mul<SVector,SVector>
Add<Mul<Scalar,SVector>,
    Mul<SVector,SVector> >

‹#›
u class Scalar {
u float const& s_;
u public:
u Scalar(float const& s) : s_(s) {}
u int size() const { return 0; }
u float operator[](int const) const { return s_; }
u };
Expression templates: Scalar

‹#›
uInside operator classes Add and Mul we must distinguish how we reference operands:
uScalar operand – by value.
uOther operands – by reference.
u
u template<typename T> struct ArgType {
u typedef T const& result;
u };
u template<> struct ArgType<Scalar> {
u typedef Scalar result;
u };
Expression templates: ArgType

‹#›
u template<typename OP1, typename OP2>
u class Mul {
u typename ArgType<OP1>::result op1_;
u typename ArgType<OP2>::result op2_;
u public:
u Mul(OP1 const& op1, OP2 const& op2) : op1_(op1), op2_(op2) {}
u int size() const { return op1.size()!=0? op1.size() : op2.size(); }
u float operator[](int const i) const { return op1_[i] * op2_[i]; }
u };
Expression templates: Mul

‹#›
u template<typename OP1, typename OP2>
u class Add {
u typename ArgType<OP1>::result op1_;
u typename ArgType<OP2>::result op2_;
u public:
u Add(OP1 const& op1, OP2 const& op2) : op1_(op1), op2_(op2) {}
u int size() const { return op1.size()!=0? op1.size() : op2.size(); }
u float operator[](int const i) const { return op1_[i] + op2_[i]; }
u };
Expression templates: Add

‹#›
uHow do we define operators? How about this:
u template<typename T1, typename T2>
u Mul<T1,T2> operator*(T1 const& u, T2 const& v) { … } // (1)
uNot good, parameters are too general, e.g.:
u struct A {} a;
u struct B : public A {} b;
u template<typename T>
u A operator*(T const& s, A const& a) { … } // (2)
u 3 * a; // Always calls (2)
u 3 * b;  // Calls (2) only when (1) is not defined.
Expression templates: Expression class

‹#›
uWe also do not want define all possible versions (too much work):
u Mul< SVector, SVector > operator*(SVector const& u, SVector const& v);
u template<typename T> Mul<SVector,Mul<T> > operator*(SVector const& u, Mul<SVector,T> const& v);
u template<typename T> Mul<Mul<T>, SVector> operator*(Mul<Svector,T> const& u, SVector const& v);
u template<typename T> Mul< SVector,Add<T> > operator*(SVector const& u, Add<SVector,T> const& v);
u …
uSolution: We introduce a class “Expr” for wrapping our expressions:
u template<typename Rep> class Expr;
uWe can safely declare the operator as:
u template<typename OP1, typename OP2>
u Expr<Mul<OP1,OP2> > operator*(
u Expr<OP1> const& u, Expr<OP2> const& v);
u
Expression templates: Expression class

‹#›
u template<typename Rep>
u class Expr { // Expression
u Rep rep_;
u public:
u explicit Expr(int size) : rep_(size) {}
u Expr(Rep const& rep) : rep_(rep) {}
u Rep const& rep() const { return rep_; } // Unwrap the instance.
u int size() const { return rep_.size(); }
u float operator[](int const i) const { return rep_[i]; }
u // Defined later:
u template<typename Rep2> Expr& operator=(Expr<Rep2> const& o);
u };
Expression templates: Expression

‹#›
u template<typename OP1, typename OP2>
u Expr<Mul<OP1,OP2> > operator*(Expr<OP1> const& u, Expr<OP2> const& v) {
u return Expr<Mul<OP1,OP2> >(Mul<OP1, OP2>(u.rep(), v.rep()));
u }
u
u template<typename OP2>
u Expr<Mul<Scalar,OP2> > operator*(float const& s, Expr<OP2> const& u) {
u return Expr<Mul<Scalar,OP2> >(Mul<Scalar, OP2>(Scalar(s), u.rep()));
u }
u
u template<typename OP1, typename OP2>
u Expr<Add<OP1,OP2> > operator+(Expr<OP1> const& u, Expr<OP2> const& v) {
u return Expr<Add<OP1,OP2> >(Add<OP1, OP2>(u.rep(), v.rep()));
u }
Expression templates: Operators
Wrap the result
Actual operation
Unwrap operands

‹#›
utemplate<typename Rep>
utemplate<typename Rep2>
uExpr<Rep>& Expr<Rep>::operator=(Expr<Rep2> const& o) {
u for (int i=0; i<size(); ++i)
u rep_[i] = o[i];
u return *this;
u}
uo[i] ==> Add[i]   (Expr<Add<…> > forwards operator[] to ‘rep’)
u       ==> Add::op1[i]               + Add::op2[i]
u ==> Mul::op1[i] * Mul::op2[i] + Mul::op1[i] * Mul::op2[i]
u ==> Scalar[i]   * SVector[i]  + SVector[i]  * SVector[i]
u ==> 1.2f        * u[i]        + u[i]        * v[i]
u
Expression templates: Expr::operator=
Add<Mul<Scalar,SVector>,
    Mul<SVector,SVector> >
Compiler converts ‘o[i]’ to
‘1.2f*u[i] + u[i]*v[i]’
=> we get the efficient version!

‹#›
uWe must declare our vector variables as:
u Expr<SVector> u(1000), v(1000); // Not nice!
uTherefore, we usually define:
u typedef Expr<SVector> Vector;
u
uLimitation: Does not work for operations where RAM temporaries are really needed, e.g., ‘x = A*x’,
where x is a vector and A is a matrix.
uBut ‘y = A*x’ would work just fine (assuming x,y are different vectors).
u
uUsage rule of thumb:
uFor small fixed-size vectors use metaprogramming, e.g., unrolling loops.
uFor large dynamic-size vectors use expression templates.
Expression templates: Notes & Limitations

‹#›
u[1] D.Vandevoorde, N.M.Josuttis; C++ Temlates, The Complete Guide; Addison-Wesley,2006.
References