Object-Oriented Programming for Scientific Computing Template - PowerPoint PPT Presentation

Object-Oriented Programming for Scientific Computing Template Metaprogramming Ole Klein Interdisciplinary Center for Scientific Computing Heidelberg University ole.klein@iwr.uni-heidelberg.de 30. Juni 2015 Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 1 / 51

Calculating the Square Root We can calculate the square root as follows using nested intervals: #include <iostream > template <std :: size_t N, std :: size_t L=1, std :: size_t H=N> struct Sqrt { public: enum{ mid = (L+H+1) /2 }; enum{ value = (N<mid*mid)? (std :: size_t)Sqrt <N,L,mid -1 >:: value : (std :: size_t)Sqrt <N,mid ,H >:: value }; }; template <std :: size_t N, std :: size_t M> struct Sqrt <N,M,M> { enum{ value = M }; }; int main () { std ::cout <<Sqrt <9 >:: value <<"�"<<Sqrt <42 >:: value <<std :: endl; } Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 2 / 51

Template Instantiations • Calculating Sqrt<9> first leads to the execution of: Sqrt <9,1,9>:: value = (9 <25) ? Sqrt <9,1,4>:: value : Sqrt <9,5,9>:: value = Sqrt <9,1,4>:: value As a result Sqrt<9,1,4> is calculated next: Sqrt <9,1,4>:: value = (9 <9) ? Sqrt <9,1,2>:: value : Sqrt <9,3,4>:: value = Sqrt <9,3,4>:: value The next recursion step is then: Sqrt <9,3,4>:: value = (9 <16) ? Sqrt <9,3,3>:: value : Sqrt <9,4,3>:: value = Sqrt <9,3,3>:: value = 3 • However, there is a problem with the ternary operator <condition>?<true-path>:<false-path> . The compiler generates not just the relevant part, but also the other that is not used. This means it has to expand the next recursion level on that side as well (although the result will be discarded), which results in the assembly of the complete binary tree. • For the square root of N this leads to 2N template instantiations. This puts a large strain on the resources of the compiler (both runtime and memory), and limits the scope of the technique. Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 3 / 51

Type Selection at Compile Time We can get rid of the unnecessary template instantiations by simply selecting the correct type and evaluating it directly. This can be carried out with a small metaprogram that corresponds to an if statement (also called a “compile time type selection”). // D e f i n i t i o n i n c l u d i n g s p e c i a l i z a t i o n f o r the t r u e case template < bool B, typename T1 , typename T2 > s t r u c t I f T h e n E l s e { t y p e d e f T1 ResultT ; } ; // P a r t i a l s p e c i a l i z a t i o n f o r the f a l s e case template < typename T1 , typename T2 > s t r u c t IfThenElse < f a l s e , T1 , T2 > { t y p e d e f T2 ResultT ; } ; #e n d i f Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 4 / 51

Improved Calculation of the Square Root Using our meta- if statement we can implement the square root as follows: template < std : : s i z e t N, std : : s i z e t L=1, std : : s i z e t H =N > s t r u c t Sqrt { p u b l i c : enum { mid = (L+H+1)/2 } ; t y p e d e f typename IfThenElse < (N < mid ∗ mid ) , Sqrt < N, L , mid − 1 > , Sqrt < N, mid ,H > > :: ResultT ResultType ; enum { v a l u e = ResultType : : v a l u e } ; } ; template < std : : s i z e t N, std : : s i z e t M > s t r u c t Sqrt < N,M,M > { enum { v a l u e = M } ; } ; This only requires about log 2 ( N ) template instantiations! Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 5 / 51

Turing Completeness of Template Metaprogramming Template meta programs may include: • State variables: the template parameters. • Loops: using recursion. • Conditional execution: using the ternary operator or template specialization (e.g. the meta- if ). • Integer calculations. This is sufficient to perform any calculation, as long as there isn’t any restriction on the number of recursive instantiations and on the number of state variables (which does not imply that it is useful to calculate everything with template metaprogramming) . Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 6 / 51

Loop Unrolling In the calculation of a scalar product, as in: template <typename T> inline T dot_product (int dim , T* a, T* b) { T result = T(); for (int i=0;i<dim ;++i) { result += a[i]*b[i]; } return result; } the compiler often optimizes the computation for large arrays. However, if small scalar products of the type dp = dot_product (3,a,b); are used most of the time, then this may be written more efficiently. Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 7 / 51

Loop Unrolling // Primary template template <int DIM , typename T> class DotProduct { public: static T result (T* a, T* b) { return *a * *b + DotProduct <DIM -1,T >:: result(a+1,b+1); } }; // Partial specialization as stopping criterion template <typename T> class DotProduct <1,T> { public: static T result (T* a, T* b) { return *a * *b; } }; Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 8 / 51

Loop Unrolling // for simplification template <int DIM , typename T> inline T dot_product (T* a, T* b) { return DotProduct <DIM ,T >:: result(a,b); } Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 9 / 51

Application: Loop Unrolling #include <iostream > #include " loop_unrolling .h" int main () { int a[3] = { 1, 2, 3}; int b[3] = { 5, 6, 7}; std :: cout << "dot_product <3>(a,b)�=�" << dot_product <3>(a,b) << ’\n’; std :: cout << "dot_product <3>(a,a)�=�" << dot_product <3>(a,a) << ’\n’; } Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 10 / 51

Loop Unrolling for Random Access Container // Primary template template <int DIM , typename T, template <typename U,typename=std :: allocator <U> > class vect > struct DotProduct { static T result(const vect <T> &a, const vect <T> &b) { return a[DIM -1]*b[DIM -1] + DotProduct <DIM -1,T,vect >:: result(a,b); } }; // Partial specialization as stopping criterion template <typename T, template <typename U,typename=std :: allocator <U> > class vect > struct DotProduct <1,T,vect > { static T result(const vect <T> &a, const vect <T> &b) { return a[0] * b[0]; } }; Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 11 / 51

Loop Unrolling for Random Access Container // For simplification template <int DIM , typename T, template <typename U,typename=std :: allocator <U> > class vect > inline T dot_product (const vect <T> &a, const vect <T> &b) { return DotProduct <DIM ,T,vect >:: result(a,b); } Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 12 / 51

Application: Loop Unrolling for Random Access Container #include <iostream > #include <vector > #include " loop_unrolling2 .h" int main () { std :: vector <double > a(3 ,3.0); std :: vector <double > b(3 ,5.0); std :: cout << "dot_product <3>(a,b)�=�" << dot_product <3>(a,b) << ’\n’; std :: cout << "dot_product <3>(a,a)�=�" << dot_product <3>(a,a) << ’\n’; } Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 13 / 51

constexpr #include <iostream > int x1 = 7; constexpr int x2 = 7; constexpr int x3 = x1; // Error , x1 is not a constexpr constexpr int x4 = x2; constexpr int Fac(int n) { return n<2 ? 1 : n * Fac(n -1); } int main () { std :: cout << Fac (10) << std :: endl; } • C++11 introduces a simple alternative to template metaprogramming: expressions that are already evaluated at compile time. • In a constexpr only variables or functions which are constexpr themselves may be used. Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 14 / 51

constexpr It must be possible to evaluate a constexpr at compile time: void f(int n) { constexpr int x = Fac(n); // Error , n isn ’t known at // time of translation int f10 = Fac (10); // Correct } int main () { const int ten = 10; int f10 = Fac (10); // Also correct } Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 15 / 51

constexpr This will work even for objects of classes whose constructor is simple enough to be defined as constexpr : struct Point { int x,y; constexpr Point(int xx , int yy) : x(xx), y(yy) {} }; constexpr Point origo (0 ,0); constexpr int z = origo.x; constexpr Point a[] = {Point (0 ,0), Point (1 ,1), Point (2 ,2)}; constexpr x = a[1].x; // x becomes 1 • constexpr functions may not have the return type void and neither variables nor functions may be defined within them (this also applies to constexpr constructors). • The function body can only contain declarations and a single return statement. Ole Klein (IWR) Object-Oriented Programming 30. Juni 2015 16 / 51