Tutorial for /expression templates.

Overview

The /expression folder contains Expression templates. For the most part, the containers in MyraMath (like Matrix and Vector) are designed to behave like regular types with value semantics. A compound expression that composes Matrix's together, like Z = X + make_complex(A,B) - 0.5*Y;, will result in the construction of several temporaries for subexpressions like 0.5*Y and make_complex(A,B). All of these temporaries are RAII-managed, and do properly release their resources when they "expire at the semicolon" after Z has been populated. However these temporaries might be frustrating for the performance-minded developer, because it's trivial to write a double-for loop with the exact same semantics that would use only O(1) additional stack space (and no heap allocation). Expression templates are a mechanism to combine the best features of both approaches: high readability via compact mathematical notation (no explicit loops), but only requiring O(1) additional space (no Matrix temporaries).

Expression templates are syntactic sugar for performing entrywise mathematical operations upon numeric containers. Expression's can be composed using the usual arithmetic operations (summation, scaling, entrywise multiplication, entrywise quotient). Many of the common transcendental functions (trignometric functions, exponential functions, etc) are overloaded for Expression's. Numeric-valued Expression's can compared to one another using relational operators (greater, less) to yield boolean-valued Expression's, which can be further composed using logical operators (and, or, not, ternary). It's also straightforward to adapt end-user code into an Expression.

Eventually, Expression's get captured into a real containers (like Matrix and Vector) by evaluate()'ing them. Additionally, most containers offer constructors and assignment operators that accept Expression's as right-hand-side values.

Below are links to the sections of this page:

• Basic Expression's
• Expression arithmetic
• Transcendental functions
• Comparisons and boolean Expression's
• User-defined Expression's
• Appendix: Regarding /sparse Expression's
• Appendix: Design philosophy, "always-on" or "opt-in"?

Basic Expression's

Conceptually, an Expression represents a mapping function from integer indices into Number values. An Expression is defined by:

• .. its "Arity", the rank of its index vector. A vector-like expression like e(i) has Arity of 1, while a matrix-like expression like E(i,j) has an Arity of 2
• .. its "size", the range of valid indices. If e is an Arity=1 expression, that reports a size() of 3, then e(0), e(1) and e(2) are well-defined but trying to evaluate any other e(i) is undefined.
• .. its "Number" type, which is baked into a typedef.
• .. its actual numeric values, which are accessed via the evaluate() member function.

Any Expression can be printed via std::cout << e, but note that you need to include expression/print.h to do so. The << operator will print the Arity/size/Number attributes, and also print the results of evaluating the Expression over all possible indices. The example program below (tutorial/expression/basic.cpp) explores some basic Expression's:

#include <myramath/dense/Vector.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/expr.h>

#include <myramath/expression/Expression.h>
#include <myramath/expression/generators.h>
#include <myramath/expression/print.h>
#include <myramath/expression/reshape.h>
#include <myramath/expression/expr.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("expression_basic","[expression]") // int main()
  {
  // Expressions that return constant values.
  auto E1 = make_ConstantExpression(3,2,1.0);
  std::cout << "make_ConstantExpression(3,2,1.0) = " << E1 << std::endl;
  // An Expression that generates random complex values.
  auto E2 = make_RandomExpression<std::complex<double> >(2,4);
  std::cout << "make_RandomExpression<std::complex<double> >(2,4) = " << E2 << std::endl;
  // An Expression that interpolates linearly between two endpoints.
  auto E3 = make_LinspaceExpression(1.0f,2.0f,11);
  std::cout << "make_LinspaceExpression(1.0f,2.0f,11) = " << E3 << std::endl;
  // Can evaluate() an Expression into a container.
  auto x = Vector<float>::evaluate(E3);
  std::cout << "Vector<float>::evaluate() = " << x << std::endl;
  // Adapting a container into an Expression.
  auto A = Matrix<double>::fill_cmajor(3,3,{1,4,2, 8,5,7, 3,9,6});
  std::cout << "A = " << A << std::endl;
  std::cout << "expr(A) = " << expr(A) << std::endl;
  // Expressions can be further manipulated without change to underlying A.
  std::cout << "reshape(expr(A),1,9) = " << reshape(expr(A),1,9) << std::endl;
  std::cout << "A = " << A << std::endl;
  // Expressions can be built from std::vector's and std::initializer_list's
  std::cout << "expr({3.0f,1.0f,2.0f}) = " << expr({3.0f,1.0f,2.0f}) << std::endl;
  }

(Click here to view the output of this program)

make_ConstantExpression(3,2,1.0) = size 3 by 2 arity-2 Expression of double:
[ 1 1 ]
[ 1 1 ]
[ 1 1 ]

make_RandomExpression<std::complex<double> >(2,4) = size 2 by 4 arity-2 Expression of std::complex<double>:
[ (0.380002,0.0108368) (0.182981,0.10957) (-0.243142,-0.484536) (-0.585236,0.252524) ]
[ (-0.319746,0.687704) (-0.862444,-0.180186) (0.759988,-0.361039) (0.961136,-0.82999) ]

make_LinspaceExpression(1.0f,2.0f,11) = size 11 arity-1 Expression of float:
[ 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 ]

Vector<float>::evaluate() = size 11 Vector of float:
[ 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 ]

A = size 3 by 3 Matrix of double:
[ 1 8 3 ]
[ 4 5 9 ]
[ 2 7 6 ]

expr(A) = size 3 by 3 arity-2 Expression of double:
[ 1 8 3 ]
[ 4 5 9 ]
[ 2 7 6 ]

reshape(expr(A),1,9) = size 1 by 9 arity-2 Expression of double:
[ 1 4 2 8 5 7 3 9 6 ]

A = size 3 by 3 Matrix of double:
[ 1 8 3 ]
[ 4 5 9 ]
[ 2 7 6 ]

expr({3.0f,1.0f,2.0f}) = size 3 arity-1 Expression of float:
[ 3 1 2 ]

Some explanatory remarks follow:

• Expression's that generate values from nothing (like ConstantExpression and RandomExpression) serve as basic building blocks for higher-level compositions. [lines 18-22]
• The semantics of the make_LinspaceExpression() [lines 28-33] and reshape() [lines 40-42] should be familiar to Matlab or Python/numpy users.
• All the containers in /dense (that is, Matrix, LowerMatrix, DiagonalMatrix and Vector) provide an evaluate() method that captures an Expression into actual numeric values. [lines 31-33]
• The reverse process, adapting a container A back into an Expression, can be accomplished by calling expr(A). [line 35-38]
• MyraMath also defines expr(std::vector<T>) and expr(std::initializer_list<T>), each returns an Arity-1 Expression of T of matching size() [line 44-45]

Expression arithmetic

Expression's can be composed with the usual arithmetic operators (+,-,*,/), but bear in mind that they always operate entrywise upon their underlying operands. For example, the product of two Expression's E1*E2 behaves like the Matlab operation E1.*E2, not like gemm(). The program below (tutorial/expression/arithmetic.cpp) shows how Expression's can be combined entrywise using the (+,-,*,/) operators:

#include <myramath/expression/expr.h>
#include <myramath/expression/reshape.h>
#include <myramath/expression/arithmetic.h>
#include <myramath/expression/print.h>

#include <tests/myratest.h>

#include <iostream>
#include <typeinfo>

using namespace myra;

ADD_TEST("expression_arithmetic","[expression]") // int main()
  {
  // Make two arbitrary Expression's by reshaping std::initializer_list's.
  auto E1 = reshape(expr({1.0,3.0,4.0,9.0,2.0,6.0}),2,3);
  auto E2 = reshape(expr({7.0,4.0,2.0,8.0,5.0,9.0}),2,3);
  std::cout << "E1 = " << E1 << std::endl;
  std::cout << "E2 = " << E2 << std::endl;
  // Examine various arithmetic compositions of E1 and E2.
  std::cout << "E1+E2 = " << E1+E2 << std::endl;
  std::cout << "E1-E2 = " << E1-E2 << std::endl;
  std::cout << "E1*E2 = " << E1*E2 << std::endl;
  std::cout << "E1/E2 = " << E1/E2 << std::endl;
  std::cout << "-E1 = " << -E1 << std::endl;
  std::cout << "1/E2 = " << 1/E2 << std::endl;
  // Can also compose an Expression with a scalar constant.
  std::cout << "2*E1 = " << 2*E1 << std::endl;
  std::cout << "E2/2 = " << E2/2 << std::endl;
  std::cout << "E1+3 = " << E1+3 << std::endl;
  std::cout << "3+E1 = " << 3+E1 << std::endl;
  std::cout << "E2-4 = " << E2-4 << std::endl;
  std::cout << "4-E2 = " << 4-E2 << std::endl;
  // Can build up arbitrary Expression's, normal operator precedence rules apply.
  std::cout << "(10-2*E1+3*E2)/4 = " << (10-2*E1+3*E2)/4 << std::endl;
  }

(Click here to view the output of this program)

E1 = size 2 by 3 arity-2 Expression of double:
[ 1 4 2 ]
[ 3 9 6 ]

E2 = size 2 by 3 arity-2 Expression of double:
[ 7 2 5 ]
[ 4 8 9 ]

E1+E2 = size 2 by 3 arity-2 Expression of double:
[ 8 6 7 ]
[ 7 17 15 ]

E1-E2 = size 2 by 3 arity-2 Expression of double:
[ -6 2 -3 ]
[ -1 1 -3 ]

E1*E2 = size 2 by 3 arity-2 Expression of double:
[ 7 8 10 ]
[ 12 72 54 ]

E1/E2 = size 2 by 3 arity-2 Expression of double:
[ 0.142857 2 0.4 ]
[ 0.75 1.125 0.666667 ]

-E1 = size 2 by 3 arity-2 Expression of double:
[ -1 -4 -2 ]
[ -3 -9 -6 ]

1/E2 = size 2 by 3 arity-2 Expression of double:
[ 0.142857 0.5 0.2 ]
[ 0.25 0.125 0.111111 ]

2*E1 = size 2 by 3 arity-2 Expression of double:
[ 2 8 4 ]
[ 6 18 12 ]

E2/2 = size 2 by 3 arity-2 Expression of double:
[ 3.5 1 2.5 ]
[ 2 4 4.5 ]

E1+3 = size 2 by 3 arity-2 Expression of double:
[ 4 7 5 ]
[ 6 12 9 ]

3+E1 = size 2 by 3 arity-2 Expression of double:
[ 4 7 5 ]
[ 6 12 9 ]

E2-4 = size 2 by 3 arity-2 Expression of double:
[ 3 -2 1 ]
[ 0 4 5 ]

4-E2 = size 2 by 3 arity-2 Expression of double:
[ -3 2 -1 ]
[ 0 -4 -5 ]

(10-2*E1+3*E2)/4 = size 2 by 3 arity-2 Expression of double:
[ 7.25 2 5.25 ]
[ 4 4 6.25 ]

Note that Expression's can also be combined (added, subtracted, etc) with scalar constants too.

Transcendental functions

Within /expression/functions.h are routines to apply transcendental functions to an Expression in an entrywise fashion. Most of them delegate to standard <cmath> functions of the same name. The most noteworthy functions are:

• Trigonometric functions, like sin(E), cos(E) and atan2(E1,E2).
• Exponential and logarithmic functions, like exp(E), pow(E,1.0/3.0) and log10(E)
• Manipulating complex values, like realpart(E), argument(E) and make_complex(Er,Ei)

The program below (tutorial/expression/functions.cpp) demonstrates how these functions behave. Remember that anything dealing with an Expression always operates in an entrywise fashion:

#include <myramath/dense/Matrix.h>
#include <myramath/dense/Vector.h>
#include <myramath/dense/expr.h>

#include <myramath/expression/Expression.h>
#include <myramath/expression/functions.h>
#include <myramath/expression/arithmetic.h>
#include <myramath/expression/expr.h>
#include <myramath/expression/print.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("expression_functions","[expression]") // int main()
  {
  // Fill Matrix A with random double's over [-pi,pi].
  double pi = std::acos(-1.0);
  auto A = pi*Matrix<double>::random(2,4);
  std::cout << "A = " << A << std::endl;
  // Fill Matrix Z with unit length complex numbers, whose phases are drawn from A.
  auto E = make_complex( cos(expr(A)), sin(expr(A)) );
  auto Z = Matrix< std::complex<double> >::evaluate(E);
  std::cout << "Z = " << Z << std::endl;
  std::cout << "abs(Z) = " << abs(expr(Z)) << std::endl;
  std::cout << "argument(Z) = " << argument(expr(Z)) << std::endl;
  // Fill vector x with perfect squares of integers.
  auto x = Vector<double>::evaluate(pow(expr({1.0,2.0,3.0,4.0,5.0}),2.0));
  std::cout << "x = " << x << std::endl;
  std::cout << "sqrt(x) = " << sqrt(expr(x)) << std::endl;
  // Fill Vector y with the powers of 10
  auto y = Vector<double>::evaluate(pow(10.0,expr({-2.0,-1.0,0.0,+1.0,+2.0})));
  std::cout << "y = " << y << std::endl;
  std::cout << "log10(y) = " << log10(expr(y)) << std::endl;
  }

(Click here to view the output of this program)

A = size 2 by 4 Matrix of double:
[ 1.19381 0.574852 -0.763854 -1.83857 ]
[ 0.0340448 0.344224 -1.52221 0.793327 ]

Z = size 2 by 4 Matrix of std::complex<double>:
[ (0.368119,0.929779) (0.839273,0.543711) (0.722175,-0.69171) (-0.264587,-0.964362) ]
[ (0.999421,0.0340382) (0.941338,0.337466) (0.0485632,-0.99882) (0.701478,0.712691) ]

abs(Z) = size 2 by 4 arity-2 Expression of double:
[ 1 1 1 1 ]
[ 1 1 1 1 ]

argument(Z) = size 2 by 4 arity-2 Expression of double:
[ 1.19381 0.574852 -0.763854 -1.83857 ]
[ 0.0340448 0.344224 -1.52221 0.793327 ]

x = size 5 Vector of double:
[ 1 4 9 16 25 ]

sqrt(x) = size 5 arity-1 Expression of double:
[ 1 2 3 4 5 ]

y = size 5 Vector of double:
[ 0.01 0.1 1 10 100 ]

log10(y) = size 5 arity-1 Expression of double:
[ -2 -1 0 1 2 ]

Note it's relatively easy to add new functions of Expression's, so just ask if you find one missing that you would like added.

Comparisons and boolean Expression's

One Expression can be compared to another (or to a constant value) using relational operators like > or <, yielding an Expression of bool values. Boolean expressions can be further composed using operators like && (logical and), || (logical or) or ! (logical not). Typically the ternary operator is used to turn bool-valued C++ expressions back into numeric values, using a syntax like bool ? trueValue : falseValue. However, the ternary operator cannot be overloaded. Instead a ternary(bool's,trueValue's,falseValue's) function is provided, where all three arguments can be Expression's of the same Arity and size.

The example program below (tutorial/expression/logical.cpp) explores logical expressions (comparison, boolean composition, and ternary operations):

#include <myramath/expression/Expression.h>
#include <myramath/expression/functions.h>
#include <myramath/expression/logical.h>
#include <myramath/expression/expr.h>
#include <myramath/expression/print.h>

#include <tests/myratest.h>

#include <myramath/dense/Vector.h>


#include <iostream>

using namespace myra;

ADD_TEST("expression_logical","[expression]") // int main()
  {
  // Compare expressions E1 and E2.
  auto E1 = expr({8.0,1.0,7.0,4.0,9.0});
  auto E2 = expr({5.0,0.0,2.0,6.0,3.0});
  std::cout << "E1 = " << E1 << std::endl;
  std::cout << "E2 = " << E2 << std::endl;  
  std::cout << "E1 > E2 = " << (E1 > E2) << std::endl;
  // Compare a function of E1 and E2.
  std::cout << "sin(E1) = " << sin(E1) << std::endl;
  std::cout << "sin(E2) = " << sin(E2) << std::endl;
  std::cout << "sin(E1) > sin(E2) = " << (sin(E1) > sin(E2)) << std::endl;
  // Compose logical expressions using && and ||.
  std::cout << "(E1 > E2) && (sin(E1) > sin(E2)) = " 
            << ((E1>E2) && (sin(E1)>sin(E2))) << std::endl;
  std::cout << "(E1 > E2) || (sin(E1) > sin(E2)) = " 
            << ((E1>E2) || (sin(E1)>sin(E2))) << std::endl;
  // Examine the ternary operator.
  std::cout << "E1>E2 ? E1 : E2 = " << ternary(E1>E2, E1, E2) << std::endl;
  // Note, E1>E2 ? E1 : E2 is equivalent to max(E1,E2), so double check.
  std::cout << "max(E1,E2) = " << max(E1,E2) << std::endl;
  }

(Click here to view the output of this program)

E1 = size 5 arity-1 Expression of double:
[ 8 1 7 4 9 ]

E2 = size 5 arity-1 Expression of double:
[ 5 0 2 6 3 ]

E1 > E2 = size 5 arity-1 Expression of bool:
[ 1 1 1 0 1 ]

sin(E1) = size 5 arity-1 Expression of double:
[ 0.989358 0.841471 0.656987 -0.756802 0.412118 ]

sin(E2) = size 5 arity-1 Expression of double:
[ -0.958924 0 0.909297 -0.279415 0.14112 ]

sin(E1) > sin(E2) = size 5 arity-1 Expression of bool:
[ 1 1 0 0 1 ]

(E1 > E2) && (sin(E1) > sin(E2)) = size 5 arity-1 Expression of bool:
[ 1 1 0 0 1 ]

(E1 > E2) || (sin(E1) > sin(E2)) = size 5 arity-1 Expression of bool:
[ 1 1 1 0 1 ]

E1>E2 ? E1 : E2 = size 5 arity-1 Expression of double:
[ 8 1 7 6 9 ]

max(E1,E2) = size 5 arity-1 Expression of double:
[ 8 1 7 6 9 ]

Boolean expressions and ternary() can be used to implement similar functionality to binary masking in Matlab.

User-defined Expression's

Advanced users may discover a need to adapt their own code or data structure as an Expression, see expression/UserExpression.h. For example, the excerpt below (tutorial/expression/user.cpp) adapts a user-defined matrix-like container into an Expression, which can be further composed using any of the routines enumerated above, or captured into a Matrix to feed into downstream algorithms.

#include <myramath/dense/Matrix.h>

#include <myramath/expression/Expression.h>
#include <myramath/expression/UserExpression.h>
#include <myramath/expression/functions.h>
#include <myramath/expression/logical.h>
#include <myramath/expression/print.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

namespace {

// A strawman 3x3 Matrix class, to be adapted into an Expression.
class MyMatrix
  {
  public:

    // Constructor, initializes with zeros.
    MyMatrix()
      {
      for (int m = 0; m < 3; ++m)
        for (int n = 0; n < 3; ++n)
          contents[m][n] = 0.0;
      }

    // Accessing/mutating the (m,n)'th entry.
    const double& operator() (int m, int n) const
      { return contents[m][n]; }
    double& operator()(int m, int n) 
      { return contents[m][n]; }

  private:
    
    // Internal storage.
    double contents[3][3];
    
  };

// For adapting a MyMatrix instance into an Expression.
class MyMatrixExpression
  {
  public:

    // User expressions must provide an Arity constant and Number typedef.
    const static int Arity = 2;
    typedef double Number;

    // User expressions must report their size. This one is 3x3.
    Index<2> size() const
      { return make_Index(3,3); }

    // User expressions must supply an evaluate() function that returns a Number.
    Number evaluate(Index<2> ij) const
      { return (*m)(ij[0],ij[1]); }

    // User expressions should capture any needed environmental references via constructor arguments.    
    MyMatrixExpression(const MyMatrix& in_m)
      : m(&in_m) { }

  private:

    // User routines can carry member data, but they must be copy-constructible.
    const MyMatrix* m;
    
  };

} // namespace

ADD_TEST("expression_user","[expression]") // int main()
  {
  // Initialize a MyMatrix instance.
  MyMatrix M;
  M(0,0) = 1.0; M(0,1) = 2.0; M(0,2) = 9.0;
  M(1,0) = 4.0; M(1,1) = 6.0; M(1,2) = 5.0;
  M(2,0) = 7.0; M(2,1) = 3.0; M(2,2) = 8.0;
  // Adapt it using make_UserExpression()
  auto E = make_UserExpression( MyMatrixExpression(M) );
  std::cout << "E = " << E << std::endl;
  // Expression E can be used within furthers compositions, and evaluate()'d
  std::cout << "E > 4.0 = " << (E > 4.0) << std::endl;
  auto A = Matrix<double>::evaluate( pow(E,2.0) );
  std::cout << "Matrix::evaluate(pow(E,2)) = " << A << std::endl;
  }

(Click here to view the output of this program)

E = size 3 by 3 arity-2 Expression of double:
[ 1 2 9 ]
[ 4 6 5 ]
[ 7 3 8 ]

E > 4.0 = size 3 by 3 arity-2 Expression of bool:
[ 0 0 1 ]
[ 0 1 1 ]
[ 1 0 1 ]

Matrix::evaluate(pow(E,2)) = size 3 by 3 Matrix of double:
[ 1 4 81 ]
[ 16 36 25 ]
[ 49 9 64 ]

User expressions should be encapsulated within a class with the following signatures:

• A typedef for the Number type that it returns (float for example, or std::complex<double>)
• A const static int that holds its Arity
• A size() method, which should return an Index<Arity>
• An evaluate(i) method, which should take an Index<Arity> i and return the corresponding Number

Appendix: Regarding /sparse Expression's

This tutorial has only discussed Expression's for /dense containers. Some of the containers in /sparse can also be adapted into Expression's using expr(). For example, a SparseMatrix can be adapted into an Arity-2 Expression of Number's, and a PatternBuilder can be adapted into an Arity-2 Expression of bool's. Any structural zeros in these containers become explicit zero/false values. Such Expression's composed using the usual Expression arithmetic, and can be captured back into /dense containers using evaluate().

But to be truthful, Expression's are not particularly useful for populating /sparse containers. The naive implementation of A = SparseMatrix::evaluate(E), which would loop over all indices in E.size(), would be asymptotically inefficient. Note this algorithm would require work proportional to the size of E, instead of work proportional to the actual number of nonzeros in E. Recalling back to the /sparse tutorial, explicit looping over all (i,j) indices was discouraged in favor of SparseMatrixIterator, for essentially the same reason. The bottom line is that SparseMatrix does not provide an evaluate(E) method. Instead, it provides an evaluate(P,E) method, where P prescribes the Pattern where E should be sampled. This tends to spoil the "temporary-free" rationale behind Expression's to begin with (because the user is forced to introduce a temporary to track the Pattern of E, whenever E does anything nontrivial like compare or sum two other sparse Expression's). Often, using a SparseMatrixBuilder is the better choice.

That said, the example below (tutorial/expression/sparse.cpp) shows how to use Expression's for /sparse containers.

#include <myramath/dense/intRange.h>

#include <myramath/sparse/SparseMatrix.h>
#include <myramath/sparse/PatternRange.h>
#include <myramath/sparse/expr.h>

#include <myramath/expression/Expression.h>
#include <myramath/expression/functions.h>
#include <myramath/expression/print.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("expression_sparse","[expression]") // int main()
  {
  // Populate a SparseMatrix from coordinate format, (i,j,value) triplets.
  auto A = SparseMatrix<double>::fill(
    { 0 ,  5 ,  1 ,  2 ,  0 ,  4 ,  1 ,  3 },    // i's
    { 0 ,  5 ,  3 ,  4 ,  1 ,  3 ,  2 ,  2 },    // j's
    {1.0, 1.7, 2.4, 3.1, 0.8, 9.2, 1.4, 4.6}  ); // value's
  std::cout << "A = " << A << std::endl;
  auto P = A.pattern();
  // Examine expr(A)
  std::cout << "expr(A) = " << expr(A) << std::endl;
  std::cout << "expr(P) = " << expr(P) << std::endl;
  // Populate B with the squares of A's values.
  auto B = SparseMatrix<double>::evaluate(P, pow(expr(A),2.0) );
  std::cout << "B = SparseMatrix::evaluate(pow(expr(A),2.0)) = " << B << std::endl;
  std::cout << "expr(B) = " << expr(B) << std::endl;
  }

(Click here to view the output of this program)

A = size 6 by 6 SparseMatrix of double (8 nz):
A(:,j0) = 1 nz [ i0,1 ]
A(:,j1) = 1 nz [ i0,0.8 ]
A(:,j2) = 2 nz [ i1,1.4 i3,4.6 ]
A(:,j3) = 2 nz [ i1,2.4 i4,9.2 ]
A(:,j4) = 1 nz [ i2,3.1 ]
A(:,j5) = 1 nz [ i5,1.7 ]

expr(A) = size 6 by 6 arity-2 Expression of double:
[ 1 0.8 0 0 0 0 ]
[ 0 0 1.4 2.4 0 0 ]
[ 0 0 0 0 3.1 0 ]
[ 0 0 4.6 0 0 0 ]
[ 0 0 0 9.2 0 0 ]
[ 0 0 0 0 0 1.7 ]

expr(P) = size 6 by 6 arity-2 Expression of bool:
[ 1 1 0 0 0 0 ]
[ 0 0 1 1 0 0 ]
[ 0 0 0 0 1 0 ]
[ 0 0 1 0 0 0 ]
[ 0 0 0 1 0 0 ]
[ 0 0 0 0 0 1 ]

B = SparseMatrix::evaluate(pow(expr(A),2.0)) = size 6 by 6 SparseMatrix of double (8 nz):
A(:,j0) = 1 nz [ i0,1 ]
A(:,j1) = 1 nz [ i0,0.64 ]
A(:,j2) = 2 nz [ i1,1.96 i3,21.16 ]
A(:,j3) = 2 nz [ i1,5.76 i4,84.64 ]
A(:,j4) = 1 nz [ i2,9.61 ]
A(:,j5) = 1 nz [ i5,2.89 ]

expr(B) = size 6 by 6 arity-2 Expression of double:
[ 1 0.64 0 0 0 0 ]
[ 0 0 1.96 5.76 0 0 ]
[ 0 0 0 0 9.61 0 ]
[ 0 0 21.16 0 0 0 ]
[ 0 0 0 84.64 0 0 ]
[ 0 0 0 0 0 2.89 ]

Appendix: Design philosophy, "always-on" or "opt-in"?

MyraMath's expression templates (ET's) are relatively weak when compared to those of other C++ linear algebra libraries (such as Eigen, MTL4, Blaze, Armadillo, or FLENS, to name a few). Within an ET framework, all kinds of performance optimizations can be embedded (loop unrolling, vectorization, algorithmic blocking). Often, such optimizations lead to C++ code that can credibly compete with vendor-tuned BLAS implementations in common benchmarks (like axpy(), gemv(), gemm()), especially at small problem sizes. MyraMath ET's do none of this, if you want to achieve vendor-tuned performance then you should call BLAS primitives (like axpy(), gemv(), gemm()) explicitly.

Note that pervasive/"always-on" use of ET's does carry drawbacks, in particular they interfere with "value-semantics". Within an ET framework a subexpression like B+C might not yield a "regular type", but rather represent a deferred calculation that is only triggered upon assignment. In particular, a statement like auto A = B+C might not do what one expects! (ET's interact poorly with auto). MyraMath places very high importance on maintaining "value-semantics" and compatibility with auto, so its Expression templates are not "always-on" but rather an "opt-in" feature accessible only by calling expr(). If B and C are explicit containers then auto A = B+C is always an explicit container too, you only get an Expression if you ask for it via auto A = expr(B)+expr(C).

MyraMath's ET's really serve just one purpose: "opt-in" syntactic sugar for populating a /dense container via evaluate() without triggering the construction of extra temporaries.

Go back to API Tutorials