1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Using CUDA (or OpenMP, TBB, ...) via Thrust</title> 5<link rel="stylesheet" href="../../../../../../../doc/src/boostbook.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../../index.html" title="Chapter 1. Boost.Numeric.Odeint"> 8<link rel="up" href="../tutorial.html" title="Tutorial"> 9<link rel="prev" href="self_expanding_lattices.html" title="Self expanding lattices"> 10<link rel="next" href="using_opencl_via_vexcl.html" title="Using OpenCL via VexCL"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../logo.jpg"></td> 15<td align="center"><a href="../../../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="self_expanding_lattices.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="using_opencl_via_vexcl.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h3 class="title"> 27<a name="boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust"></a><a class="link" href="using_cuda__or_openmp__tbb_______via_thrust.html" title="Using CUDA (or OpenMP, TBB, ...) via Thrust">Using 28 CUDA (or OpenMP, TBB, ...) via Thrust</a> 29</h3></div></div></div> 30<div class="toc"><dl class="toc"> 31<dt><span class="section"><a href="using_cuda__or_openmp__tbb_______via_thrust.html#boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust.phase_oscillator_ensemble">Phase 32 oscillator ensemble</a></span></dt> 33<dt><span class="section"><a href="using_cuda__or_openmp__tbb_______via_thrust.html#boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust.large_oscillator_chains">Large 34 oscillator chains</a></span></dt> 35<dt><span class="section"><a href="using_cuda__or_openmp__tbb_______via_thrust.html#boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust.parameter_studies">Parameter 36 studies</a></span></dt> 37</dl></div> 38<p> 39 Modern graphic cards (graphic processing units - GPUs) can be used to speed 40 up the performance of time consuming algorithms by means of massive parallelization. 41 They are designed to execute many operations in parallel. odeint can utilize 42 the power of GPUs by means of CUDA and <a href="http://code.google.com/p/thrust/" target="_top">Thrust</a>, 43 which is a STL-like interface for the native CUDA API. 44 </p> 45<div class="important"><table border="0" summary="Important"> 46<tr> 47<td rowspan="2" align="center" valign="top" width="25"><img alt="[Important]" src="../../../../../../../doc/src/images/important.png"></td> 48<th align="left">Important</th> 49</tr> 50<tr><td align="left" valign="top"><p> 51 Thrust also supports parallelization using OpenMP and Intel Threading Building 52 Blocks (TBB). You can switch between CUDA, OpenMP and TBB parallelizations 53 by a simple compiler switch. Hence, this also provides an easy way to get 54 basic OpenMP parallelization into odeint. The examples discussed below 55 are focused on GPU parallelization, though. 56 </p></td></tr> 57</table></div> 58<p> 59 To use odeint with CUDA a few points have to be taken into account. First 60 of all, the problem has to be well chosen. It makes absolutely no sense to 61 try to parallelize the code for a three dimensional system, it is simply 62 too small and not worth the effort. One single function call (kernel execution) 63 on the GPU is slow but you can do the operation on a huge set of data with 64 only one call. We have experienced that the vector size over which is parallelized 65 should be of the order of <span class="emphasis"><em>10<sup>6</sup></em></span> to make full use of the 66 GPU. Secondly, you have to use <a href="http://code.google.com/p/thrust/" target="_top">Thrust</a>'s 67 algorithms and functors when implementing the rhs the ODE. This might be 68 tricky since it involves some kind of functional programming knowledge. 69 </p> 70<p> 71 Typical applications for CUDA and odeint are large systems, like lattices 72 or discretizations of PDE, and parameter studies. We introduce now three 73 examples which show how the power of GPUs can be used in combination with 74 odeint. 75 </p> 76<div class="important"><table border="0" summary="Important"> 77<tr> 78<td rowspan="2" align="center" valign="top" width="25"><img alt="[Important]" src="../../../../../../../doc/src/images/important.png"></td> 79<th align="left">Important</th> 80</tr> 81<tr><td align="left" valign="top"><p> 82 The full power of CUDA is only available for really large systems where 83 the number of coupled ordinary differential equations is of order <span class="emphasis"><em>N=10<sup>6</sup></em></span> 84 or larger. For smaller systems the CPU is usually much faster. You can 85 also integrate an ensemble of different uncoupled ODEs in parallel as shown 86 in the last example. 87 </p></td></tr> 88</table></div> 89<div class="section"> 90<div class="titlepage"><div><div><h4 class="title"> 91<a name="boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust.phase_oscillator_ensemble"></a><a class="link" href="using_cuda__or_openmp__tbb_______via_thrust.html#boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust.phase_oscillator_ensemble" title="Phase oscillator ensemble">Phase 92 oscillator ensemble</a> 93</h4></div></div></div> 94<p> 95 The first example is the phase oscillator ensemble from the previous section: 96 </p> 97<p> 98 <span class="emphasis"><em>dφ<sub>k</sub> / dt = ω<sub>k</sub> + ε / N Σ<sub>j</sub> sin( φ<sub>j</sub> - φ<sub>k</sub> ).</em></span> 99 </p> 100<p> 101 It has a phase transition at <span class="emphasis"><em>ε = 2</em></span> in the limit of infinite 102 numbers of oscillators <span class="emphasis"><em>N</em></span>. In the case of finite <span class="emphasis"><em>N</em></span> 103 this transition is smeared out but still clearly visible. 104 </p> 105<p> 106 <a href="http://code.google.com/p/thrust/" target="_top">Thrust</a> and CUDA are 107 perfectly suited for such kinds of problems where one needs a large number 108 of particles (oscillators). We start by defining the state type which is 109 a <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">device_vector</span></code>. The content of this vector 110 lives on the GPU. If you are not familiar with this we recommend reading 111 the <span class="emphasis"><em>Getting started</em></span> section on the <a href="http://code.google.com/p/thrust/" target="_top">Thrust</a> 112 website. 113 </p> 114<p> 115</p> 116<pre class="programlisting"><span class="comment">//change this to float if your device does not support double computation</span> 117<span class="keyword">typedef</span> <span class="keyword">double</span> <span class="identifier">value_type</span><span class="special">;</span> 118 119<span class="comment">//change this to host_vector< ... > of you want to run on CPU</span> 120<span class="keyword">typedef</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">device_vector</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">state_type</span><span class="special">;</span> 121<span class="comment">// typedef thrust::host_vector< value_type > state_type;</span> 122</pre> 123<p> 124 </p> 125<p> 126 Thrust follows a functional programming approach. If you want to perform 127 a calculation on the GPU you usually have to call a global function like 128 <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">for_each</span></code>, <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">reduce</span></code>, 129 ... with an appropriate local functor which performs the basic operation. 130 An example is 131</p> 132<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">add_two</span> 133<span class="special">{</span> 134 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">T</span> <span class="special">></span> 135 <span class="identifier">__host__</span> <span class="identifier">__device__</span> 136 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="identifier">T</span> <span class="special">&</span><span class="identifier">t</span> <span class="special">)</span> <span class="keyword">const</span> 137 <span class="special">{</span> 138 <span class="identifier">t</span> <span class="special">+=</span> <span class="identifier">T</span><span class="special">(</span> <span class="number">2</span> <span class="special">);</span> 139 <span class="special">}</span> 140<span class="special">};</span> 141 142<span class="comment">// ...</span> 143 144<span class="identifier">thrust</span><span class="special">::</span><span class="identifier">for_each</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">add_two</span><span class="special">()</span> <span class="special">);</span> 145</pre> 146<p> 147 This code generically adds two to every element in the container <code class="computeroutput"><span class="identifier">x</span></code>. 148 </p> 149<p> 150 For the purpose of integrating the phase oscillator ensemble we need 151 </p> 152<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 153<li class="listitem"> 154 to calculate the system function, hence the r.h.s. of the ODE. 155 </li> 156<li class="listitem"> 157 this involves computing the mean field of the oscillator example, i.e. 158 the values of <span class="emphasis"><em>R</em></span> and <span class="emphasis"><em>θ</em></span> 159 </li> 160</ul></div> 161<p> 162 The mean field is calculated in a class <code class="computeroutput"><span class="identifier">mean_field_calculator</span></code> 163 </p> 164<p> 165</p> 166<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">mean_field_calculator</span> 167<span class="special">{</span> 168 <span class="keyword">struct</span> <span class="identifier">sin_functor</span> <span class="special">:</span> <span class="keyword">public</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">unary_function</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">></span> 169 <span class="special">{</span> 170 <span class="identifier">__host__</span> <span class="identifier">__device__</span> 171 <span class="identifier">value_type</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="identifier">value_type</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span> 172 <span class="special">{</span> 173 <span class="keyword">return</span> <span class="identifier">sin</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">);</span> 174 <span class="special">}</span> 175 <span class="special">};</span> 176 177 <span class="keyword">struct</span> <span class="identifier">cos_functor</span> <span class="special">:</span> <span class="keyword">public</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">unary_function</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">></span> 178 <span class="special">{</span> 179 <span class="identifier">__host__</span> <span class="identifier">__device__</span> 180 <span class="identifier">value_type</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="identifier">value_type</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span> 181 <span class="special">{</span> 182 <span class="keyword">return</span> <span class="identifier">cos</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">);</span> 183 <span class="special">}</span> 184 <span class="special">};</span> 185 186 <span class="keyword">static</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">get_mean</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">)</span> 187 <span class="special">{</span> 188 <span class="identifier">value_type</span> <span class="identifier">sin_sum</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">reduce</span><span class="special">(</span> 189 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_transform_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">sin_functor</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> 190 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_transform_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">sin_functor</span><span class="special">()</span> <span class="special">)</span> <span class="special">);</span> 191 <span class="identifier">value_type</span> <span class="identifier">cos_sum</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">reduce</span><span class="special">(</span> 192 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_transform_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">cos_functor</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> 193 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_transform_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">cos_functor</span><span class="special">()</span> <span class="special">)</span> <span class="special">);</span> 194 195 <span class="identifier">cos_sum</span> <span class="special">/=</span> <span class="identifier">value_type</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">size</span><span class="special">()</span> <span class="special">);</span> 196 <span class="identifier">sin_sum</span> <span class="special">/=</span> <span class="identifier">value_type</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">size</span><span class="special">()</span> <span class="special">);</span> 197 198 <span class="identifier">value_type</span> <span class="identifier">K</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span> <span class="identifier">cos_sum</span> <span class="special">*</span> <span class="identifier">cos_sum</span> <span class="special">+</span> <span class="identifier">sin_sum</span> <span class="special">*</span> <span class="identifier">sin_sum</span> <span class="special">);</span> 199 <span class="identifier">value_type</span> <span class="identifier">Theta</span> <span class="special">=</span> <span class="identifier">atan2</span><span class="special">(</span> <span class="identifier">sin_sum</span> <span class="special">,</span> <span class="identifier">cos_sum</span> <span class="special">);</span> 200 201 <span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">make_pair</span><span class="special">(</span> <span class="identifier">K</span> <span class="special">,</span> <span class="identifier">Theta</span> <span class="special">);</span> 202 <span class="special">}</span> 203<span class="special">};</span> 204</pre> 205<p> 206 </p> 207<p> 208 Inside this class two member structures <code class="computeroutput"><span class="identifier">sin_functor</span></code> 209 and <code class="computeroutput"><span class="identifier">cos_functor</span></code> are defined. 210 They compute the sine and the cosine of a value and they are used within 211 a transform iterator to calculate the sum of <span class="emphasis"><em>sin(φ<sub>k</sub>)</em></span> 212 and <span class="emphasis"><em>cos(φ<sub>k</sub>)</em></span>. The classifiers <code class="computeroutput"><span class="identifier">__host__</span></code> 213 and <code class="computeroutput"><span class="identifier">__device__</span></code> are CUDA 214 specific and define a function or operator which can be executed on the 215 GPU as well as on the CPU. The line 216 </p> 217<p> 218</p> 219<pre class="programlisting"><span class="identifier">value_type</span> <span class="identifier">sin_sum</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">reduce</span><span class="special">(</span> 220 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_transform_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">sin_functor</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> 221 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_transform_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">sin_functor</span><span class="special">()</span> <span class="special">)</span> <span class="special">);</span> 222</pre> 223<p> 224 </p> 225<p> 226 performs the calculation of this sine-sum on the GPU (or on the CPU, depending 227 on your thrust configuration). 228 </p> 229<p> 230 The system function is defined via 231 </p> 232<p> 233</p> 234<pre class="programlisting"><span class="keyword">class</span> <span class="identifier">phase_oscillator_ensemble</span> 235<span class="special">{</span> 236 237<span class="keyword">public</span><span class="special">:</span> 238 239 <span class="keyword">struct</span> <span class="identifier">sys_functor</span> 240 <span class="special">{</span> 241 <span class="identifier">value_type</span> <span class="identifier">m_K</span> <span class="special">,</span> <span class="identifier">m_Theta</span> <span class="special">,</span> <span class="identifier">m_epsilon</span><span class="special">;</span> 242 243 <span class="identifier">sys_functor</span><span class="special">(</span> <span class="identifier">value_type</span> <span class="identifier">K</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="identifier">Theta</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="identifier">epsilon</span> <span class="special">)</span> 244 <span class="special">:</span> <span class="identifier">m_K</span><span class="special">(</span> <span class="identifier">K</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_Theta</span><span class="special">(</span> <span class="identifier">Theta</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_epsilon</span><span class="special">(</span> <span class="identifier">epsilon</span> <span class="special">)</span> <span class="special">{</span> <span class="special">}</span> 245 246 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">Tuple</span> <span class="special">></span> 247 <span class="identifier">__host__</span> <span class="identifier">__device__</span> 248 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="identifier">Tuple</span> <span class="identifier">t</span> <span class="special">)</span> 249 <span class="special">{</span> 250 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span><span class="number">2</span><span class="special">>(</span><span class="identifier">t</span><span class="special">)</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span><span class="number">1</span><span class="special">>(</span><span class="identifier">t</span><span class="special">)</span> <span class="special">+</span> <span class="identifier">m_epsilon</span> <span class="special">*</span> <span class="identifier">m_K</span> <span class="special">*</span> <span class="identifier">sin</span><span class="special">(</span> <span class="identifier">m_Theta</span> <span class="special">-</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span><span class="number">0</span><span class="special">>(</span><span class="identifier">t</span><span class="special">)</span> <span class="special">);</span> 251 <span class="special">}</span> 252 <span class="special">};</span> 253 254 <span class="comment">// ...</span> 255 256 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()</span> <span class="special">(</span> <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">dt</span> <span class="special">)</span> <span class="keyword">const</span> 257 <span class="special">{</span> 258 <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">mean_field</span> <span class="special">=</span> <span class="identifier">mean_field_calculator</span><span class="special">::</span><span class="identifier">get_mean</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">);</span> 259 260 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">for_each</span><span class="special">(</span> 261 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_zip_iterator</span><span class="special">(</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">m_omega</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">dxdt</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">)</span> <span class="special">),</span> 262 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_zip_iterator</span><span class="special">(</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">m_omega</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">dxdt</span><span class="special">.</span><span class="identifier">end</span><span class="special">())</span> <span class="special">)</span> <span class="special">,</span> 263 <span class="identifier">sys_functor</span><span class="special">(</span> <span class="identifier">mean_field</span><span class="special">.</span><span class="identifier">first</span> <span class="special">,</span> <span class="identifier">mean_field</span><span class="special">.</span><span class="identifier">second</span> <span class="special">,</span> <span class="identifier">m_epsilon</span> <span class="special">)</span> 264 <span class="special">);</span> 265 <span class="special">}</span> 266 267 <span class="comment">// ...</span> 268<span class="special">};</span> 269</pre> 270<p> 271 </p> 272<p> 273 This class is used within the <code class="computeroutput"><span class="identifier">do_step</span></code> 274 and <code class="computeroutput"><span class="identifier">integrate</span></code> method. It 275 defines a member structure <code class="computeroutput"><span class="identifier">sys_functor</span></code> 276 for the r.h.s. of each individual oscillator and the <code class="computeroutput"><span class="keyword">operator</span><span class="special">()</span></code> for the use in the steppers and integrators 277 of odeint. The functor computes first the mean field of <span class="emphasis"><em>φ<sub>k</sub></em></span> 278 and secondly calculates the whole r.h.s. of the ODE using this mean field. 279 Note, how nicely <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">tuple</span></code> 280 and <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">zip_iterator</span></code> play together. 281 </p> 282<p> 283 Now we are ready to put everything together. All we have to do for making 284 odeint ready for using the GPU is to parametrize the stepper with the 285 <code class="computeroutput"><span class="identifier">state_type</span></code> and <code class="computeroutput"><span class="identifier">value_type</span></code>: 286 </p> 287<p> 288</p> 289<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">runge_kutta4</span><span class="special"><</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">stepper_type</span><span class="special">;</span> 290</pre> 291<p> 292 </p> 293<div class="note"><table border="0" summary="Note"> 294<tr> 295<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../../../doc/src/images/note.png"></td> 296<th align="left">Note</th> 297</tr> 298<tr><td align="left" valign="top"><p> 299 We have specifically define four template parameters because we have 300 to override the default parameter value <code class="computeroutput"><span class="keyword">double</span></code> 301 with <code class="computeroutput"><span class="identifier">value_type</span></code> to ensure 302 our programs runs properly if we use <code class="computeroutput"><span class="keyword">float</span></code> 303 as fundamental data type. 304 </p></td></tr> 305</table></div> 306<p> 307 You can also use a controlled or dense output stepper, e.g. 308 </p> 309<p> 310</p> 311<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">runge_kutta_dopri5</span><span class="special"><</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">stepper_type</span><span class="special">;</span> 312</pre> 313<p> 314 </p> 315<p> 316 Then, it is straightforward to integrate the phase ensemble by creating 317 an instance of the rhs class and using an integration function: 318 </p> 319<p> 320</p> 321<pre class="programlisting"><span class="identifier">phase_oscillator_ensemble</span> <span class="identifier">ensemble</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">,</span> <span class="number">1.0</span> <span class="special">);</span> 322</pre> 323<p> 324 </p> 325<p> 326</p> 327<pre class="programlisting"><span class="identifier">size_t</span> <span class="identifier">steps1</span> <span class="special">=</span> <span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">make_controlled</span><span class="special">(</span> <span class="number">1.0e-6</span> <span class="special">,</span> <span class="number">1.0e-6</span> <span class="special">,</span> <span class="identifier">stepper_type</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">ref</span><span class="special">(</span> <span class="identifier">ensemble</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="identifier">t_transients</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">);</span> 328</pre> 329<p> 330 </p> 331<p> 332 We have to use <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">ref</span></code> here in order to pass the rhs class 333 as reference and not by value. This ensures that the natural frequencies 334 of each oscillator are not copied when calling <code class="computeroutput"><span class="identifier">integrate_const</span></code>. 335 In the full example the performance and results of the Runge-Kutta-4 and 336 the Dopri5 solver are compared. 337 </p> 338<p> 339 The full example can be found at <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/thrust/phase_oscillator_ensemble.cu" target="_top">phase_oscillator_example.cu</a>. 340 </p> 341</div> 342<div class="section"> 343<div class="titlepage"><div><div><h4 class="title"> 344<a name="boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust.large_oscillator_chains"></a><a class="link" href="using_cuda__or_openmp__tbb_______via_thrust.html#boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust.large_oscillator_chains" title="Large oscillator chains">Large 345 oscillator chains</a> 346</h4></div></div></div> 347<p> 348 The next example is a large, one-dimensional chain of nearest-neighbor 349 coupled phase oscillators with the following equations of motion: 350 </p> 351<p> 352 <span class="emphasis"><em>d φ<sub>k</sub> / dt = ω<sub>k</sub> + sin( φ<sub>k+1</sub> - φ<sub>k</sub> ) + sin( φ<sub>k</sub> - φ<sub>k-1</sub>)</em></span> 353 </p> 354<p> 355 In principle we can use all the techniques from the previous phase oscillator 356 ensemble example, but we have to take special care about the coupling of 357 the oscillators. To efficiently implement the coupling you can use a very 358 elegant way employing Thrust's permutation iterator. A permutation iterator 359 behaves like a normal iterator on a vector but it does not iterate along 360 the usual order of the elements. It rather iterates along some permutation 361 of the elements defined by some index map. To realize the nearest neighbor 362 coupling we create one permutation iterator which travels one step behind 363 a usual iterator and another permutation iterator which travels one step 364 in front. The full system class is: 365 </p> 366<p> 367</p> 368<pre class="programlisting"><span class="comment">//change this to host_vector< ... > if you want to run on CPU</span> 369<span class="keyword">typedef</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">device_vector</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">state_type</span><span class="special">;</span> 370<span class="keyword">typedef</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">device_vector</span><span class="special"><</span> <span class="identifier">size_t</span> <span class="special">></span> <span class="identifier">index_vector_type</span><span class="special">;</span> 371<span class="comment">//typedef thrust::host_vector< value_type > state_type;</span> 372<span class="comment">//typedef thrust::host_vector< size_t > index_vector_type;</span> 373 374<span class="keyword">class</span> <span class="identifier">phase_oscillators</span> 375<span class="special">{</span> 376 377<span class="keyword">public</span><span class="special">:</span> 378 379 <span class="keyword">struct</span> <span class="identifier">sys_functor</span> 380 <span class="special">{</span> 381 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">Tuple</span> <span class="special">></span> 382 <span class="identifier">__host__</span> <span class="identifier">__device__</span> 383 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="identifier">Tuple</span> <span class="identifier">t</span> <span class="special">)</span> <span class="comment">// this functor works on tuples of values</span> 384 <span class="special">{</span> 385 <span class="comment">// first, unpack the tuple into value, neighbors and omega</span> 386 <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">phi</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span><span class="number">0</span><span class="special">>(</span><span class="identifier">t</span><span class="special">);</span> 387 <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">phi_left</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span><span class="number">1</span><span class="special">>(</span><span class="identifier">t</span><span class="special">);</span> <span class="comment">// left neighbor</span> 388 <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">phi_right</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span><span class="number">2</span><span class="special">>(</span><span class="identifier">t</span><span class="special">);</span> <span class="comment">// right neighbor</span> 389 <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">omega</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span><span class="number">3</span><span class="special">>(</span><span class="identifier">t</span><span class="special">);</span> 390 <span class="comment">// the dynamical equation</span> 391 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span><span class="number">4</span><span class="special">>(</span><span class="identifier">t</span><span class="special">)</span> <span class="special">=</span> <span class="identifier">omega</span> <span class="special">+</span> <span class="identifier">sin</span><span class="special">(</span> <span class="identifier">phi_right</span> <span class="special">-</span> <span class="identifier">phi</span> <span class="special">)</span> <span class="special">+</span> <span class="identifier">sin</span><span class="special">(</span> <span class="identifier">phi</span> <span class="special">-</span> <span class="identifier">phi_left</span> <span class="special">);</span> 392 <span class="special">}</span> 393 <span class="special">};</span> 394 395 <span class="identifier">phase_oscillators</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">omega</span> <span class="special">)</span> 396 <span class="special">:</span> <span class="identifier">m_omega</span><span class="special">(</span> <span class="identifier">omega</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_N</span><span class="special">(</span> <span class="identifier">omega</span><span class="special">.</span><span class="identifier">size</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_prev</span><span class="special">(</span> <span class="identifier">omega</span><span class="special">.</span><span class="identifier">size</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_next</span><span class="special">(</span> <span class="identifier">omega</span><span class="special">.</span><span class="identifier">size</span><span class="special">()</span> <span class="special">)</span> 397 <span class="special">{</span> 398 <span class="comment">// build indices pointing to left and right neighbours</span> 399 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">counting_iterator</span><span class="special"><</span><span class="identifier">size_t</span><span class="special">></span> <span class="identifier">c</span><span class="special">(</span> <span class="number">0</span> <span class="special">);</span> 400 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">copy</span><span class="special">(</span> <span class="identifier">c</span> <span class="special">,</span> <span class="identifier">c</span><span class="special">+</span><span class="identifier">m_N</span><span class="special">-</span><span class="number">1</span> <span class="special">,</span> <span class="identifier">m_prev</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()+</span><span class="number">1</span> <span class="special">);</span> 401 <span class="identifier">m_prev</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> <span class="comment">// m_prev = { 0 , 0 , 1 , 2 , 3 , ... , N-1 }</span> 402 403 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">copy</span><span class="special">(</span> <span class="identifier">c</span><span class="special">+</span><span class="number">1</span> <span class="special">,</span> <span class="identifier">c</span><span class="special">+</span><span class="identifier">m_N</span> <span class="special">,</span> <span class="identifier">m_next</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">);</span> 404 <span class="identifier">m_next</span><span class="special">[</span><span class="identifier">m_N</span><span class="special">-</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">m_N</span><span class="special">-</span><span class="number">1</span><span class="special">;</span> <span class="comment">// m_next = { 1 , 2 , 3 , ... , N-1 , N-1 }</span> 405 <span class="special">}</span> 406 407 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()</span> <span class="special">(</span> <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">dt</span> <span class="special">)</span> 408 <span class="special">{</span> 409 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">for_each</span><span class="special">(</span> 410 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_zip_iterator</span><span class="special">(</span> 411 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span> 412 <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> 413 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_permutation_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">m_prev</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> 414 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_permutation_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">m_next</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> 415 <span class="identifier">m_omega</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> 416 <span class="identifier">dxdt</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> 417 <span class="special">)</span> <span class="special">),</span> 418 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_zip_iterator</span><span class="special">(</span> 419 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span> 420 <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> 421 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_permutation_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">m_prev</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> 422 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_permutation_iterator</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">m_next</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> 423 <span class="identifier">m_omega</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> 424 <span class="identifier">dxdt</span><span class="special">.</span><span class="identifier">end</span><span class="special">())</span> <span class="special">)</span> <span class="special">,</span> 425 <span class="identifier">sys_functor</span><span class="special">()</span> 426 <span class="special">);</span> 427 <span class="special">}</span> 428 429<span class="keyword">private</span><span class="special">:</span> 430 431 <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">m_omega</span><span class="special">;</span> 432 <span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">m_N</span><span class="special">;</span> 433 <span class="identifier">index_vector_type</span> <span class="identifier">m_prev</span><span class="special">;</span> 434 <span class="identifier">index_vector_type</span> <span class="identifier">m_next</span><span class="special">;</span> 435<span class="special">};</span> 436</pre> 437<p> 438 </p> 439<p> 440 Note, how easy you can obtain the value for the left and right neighboring 441 oscillator in the system functor using the permutation iterators. But, 442 the call of the <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">for_each</span></code> 443 function looks relatively complicated. Every term of the r.h.s. of the 444 ODE is resembled by one iterator packed in exactly the same way as it is 445 unpacked in the functor above. 446 </p> 447<p> 448 Now we put everything together. We create random initial conditions and 449 decreasing frequencies such that we should get synchronization. We copy 450 the frequencies and the initial conditions onto the device and finally 451 initialize and perform the integration. As result we simply write out the 452 current state, hence the phase of each oscillator. 453 </p> 454<p> 455</p> 456<pre class="programlisting"><span class="comment">// create initial conditions and omegas on host:</span> 457<span class="identifier">vector</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">x_host</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">);</span> 458<span class="identifier">vector</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">omega_host</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">);</span> 459<span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">i</span><span class="special"><</span><span class="identifier">N</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span> 460<span class="special">{</span> 461 <span class="identifier">x_host</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">=</span> <span class="number">2.0</span> <span class="special">*</span> <span class="identifier">pi</span> <span class="special">*</span> <span class="identifier">drand48</span><span class="special">();</span> 462 <span class="identifier">omega_host</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">=</span> <span class="special">(</span> <span class="identifier">N</span> <span class="special">-</span> <span class="identifier">i</span> <span class="special">)</span> <span class="special">*</span> <span class="identifier">epsilon</span><span class="special">;</span> <span class="comment">// decreasing frequencies</span> 463<span class="special">}</span> 464 465<span class="comment">// copy to device</span> 466<span class="identifier">state_type</span> <span class="identifier">x</span> <span class="special">=</span> <span class="identifier">x_host</span><span class="special">;</span> 467<span class="identifier">state_type</span> <span class="identifier">omega</span> <span class="special">=</span> <span class="identifier">omega_host</span><span class="special">;</span> 468 469<span class="comment">// create stepper</span> 470<span class="identifier">runge_kutta4</span><span class="special"><</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">stepper</span><span class="special">;</span> 471 472<span class="comment">// create phase oscillator system function</span> 473<span class="identifier">phase_oscillators</span> <span class="identifier">sys</span><span class="special">(</span> <span class="identifier">omega</span> <span class="special">);</span> 474 475<span class="comment">// integrate</span> 476<span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">stepper</span> <span class="special">,</span> <span class="identifier">sys</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">10.0</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">);</span> 477 478<span class="identifier">thrust</span><span class="special">::</span><span class="identifier">copy</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">ostream_iterator</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">>(</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special">,</span> <span class="string">"\n"</span> <span class="special">)</span> <span class="special">);</span> 479<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 480</pre> 481<p> 482 </p> 483<p> 484 The full example can be found at <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/thrust/phase_oscillator_chain.cu" target="_top">phase_oscillator_chain.cu</a>. 485 </p> 486</div> 487<div class="section"> 488<div class="titlepage"><div><div><h4 class="title"> 489<a name="boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust.parameter_studies"></a><a class="link" href="using_cuda__or_openmp__tbb_______via_thrust.html#boost_numeric_odeint.tutorial.using_cuda__or_openmp__tbb_______via_thrust.parameter_studies" title="Parameter studies">Parameter 490 studies</a> 491</h4></div></div></div> 492<p> 493 Another important use case for <a href="http://code.google.com/p/thrust/" target="_top">Thrust</a> 494 and CUDA are parameter studies of relatively small systems. Consider for 495 example the three-dimensional Lorenz system from the chaotic systems example 496 in the previous section which has three parameters. If you want to study 497 the behavior of this system for different parameters you usually have to 498 integrate the system for many parameter values. Using thrust and odeint 499 you can do this integration in parallel, hence you integrate a whole ensemble 500 of Lorenz systems where each individual realization has a different parameter 501 value. 502 </p> 503<p> 504 In the following we will show how you can use <a href="http://code.google.com/p/thrust/" target="_top">Thrust</a> 505 to integrate the above mentioned ensemble of Lorenz systems. We will vary 506 only the parameter <span class="emphasis"><em>β</em></span> but it is straightforward to vary 507 other parameters or even two or all three parameters. Furthermore, we will 508 use the largest Lyapunov exponent to quantify the behavior of the system 509 (chaoticity). 510 </p> 511<p> 512 We start by defining the range of the parameters we want to study. The 513 state_type is again a <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">device_vector</span><span class="special"><</span> <span class="identifier">value_type</span> 514 <span class="special">></span></code>. 515 </p> 516<p> 517</p> 518<pre class="programlisting"><span class="identifier">vector</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">beta_host</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">);</span> 519<span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">beta_min</span> <span class="special">=</span> <span class="number">0.0</span> <span class="special">,</span> <span class="identifier">beta_max</span> <span class="special">=</span> <span class="number">56.0</span><span class="special">;</span> 520<span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">i</span><span class="special"><</span><span class="identifier">N</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span> 521 <span class="identifier">beta_host</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">beta_min</span> <span class="special">+</span> <span class="identifier">value_type</span><span class="special">(</span> <span class="identifier">i</span> <span class="special">)</span> <span class="special">*</span> <span class="special">(</span> <span class="identifier">beta_max</span> <span class="special">-</span> <span class="identifier">beta_min</span> <span class="special">)</span> <span class="special">/</span> <span class="identifier">value_type</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">-</span> <span class="number">1</span> <span class="special">);</span> 522 523<span class="identifier">state_type</span> <span class="identifier">beta</span> <span class="special">=</span> <span class="identifier">beta_host</span><span class="special">;</span> 524</pre> 525<p> 526 </p> 527<p> 528 The next thing we have to implement is the Lorenz system without perturbations. 529 Later, a system with perturbations is also implemented in order to calculate 530 the Lyapunov exponent. We will use an ansatz where each device function 531 calculates one particular realization of the Lorenz ensemble 532 </p> 533<p> 534</p> 535<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">lorenz_system</span> 536<span class="special">{</span> 537 <span class="keyword">struct</span> <span class="identifier">lorenz_functor</span> 538 <span class="special">{</span> 539 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">T</span> <span class="special">></span> 540 <span class="identifier">__host__</span> <span class="identifier">__device__</span> 541 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="identifier">T</span> <span class="identifier">t</span> <span class="special">)</span> <span class="keyword">const</span> 542 <span class="special">{</span> 543 <span class="comment">// unpack the parameter we want to vary and the Lorenz variables</span> 544 <span class="identifier">value_type</span> <span class="identifier">R</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">3</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">);</span> 545 <span class="identifier">value_type</span> <span class="identifier">x</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">0</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">);</span> 546 <span class="identifier">value_type</span> <span class="identifier">y</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">1</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">);</span> 547 <span class="identifier">value_type</span> <span class="identifier">z</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">2</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">);</span> 548 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">4</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">)</span> <span class="special">=</span> <span class="identifier">sigma</span> <span class="special">*</span> <span class="special">(</span> <span class="identifier">y</span> <span class="special">-</span> <span class="identifier">x</span> <span class="special">);</span> 549 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">5</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">)</span> <span class="special">=</span> <span class="identifier">R</span> <span class="special">*</span> <span class="identifier">x</span> <span class="special">-</span> <span class="identifier">y</span> <span class="special">-</span> <span class="identifier">x</span> <span class="special">*</span> <span class="identifier">z</span><span class="special">;</span> 550 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">6</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">)</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">b</span> <span class="special">*</span> <span class="identifier">z</span> <span class="special">+</span> <span class="identifier">x</span> <span class="special">*</span> <span class="identifier">y</span> <span class="special">;</span> 551 552 <span class="special">}</span> 553 <span class="special">};</span> 554 555 <span class="identifier">lorenz_system</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">N</span> <span class="special">,</span> <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">beta</span> <span class="special">)</span> 556 <span class="special">:</span> <span class="identifier">m_N</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_beta</span><span class="special">(</span> <span class="identifier">beta</span> <span class="special">)</span> <span class="special">{</span> <span class="special">}</span> 557 558 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">State</span> <span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Deriv</span> <span class="special">></span> 559 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">State</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">Deriv</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="identifier">t</span> <span class="special">)</span> <span class="keyword">const</span> 560 <span class="special">{</span> 561 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">for_each</span><span class="special">(</span> 562 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_zip_iterator</span><span class="special">(</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span> 563 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">)</span> <span class="special">,</span> 564 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">)</span> <span class="special">+</span> <span class="identifier">m_N</span> <span class="special">,</span> 565 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">)</span> <span class="special">+</span> <span class="number">2</span> <span class="special">*</span> <span class="identifier">m_N</span> <span class="special">,</span> 566 <span class="identifier">m_beta</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> 567 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">dxdt</span> <span class="special">)</span> <span class="special">,</span> 568 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">dxdt</span> <span class="special">)</span> <span class="special">+</span> <span class="identifier">m_N</span> <span class="special">,</span> 569 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">dxdt</span> <span class="special">)</span> <span class="special">+</span> <span class="number">2</span> <span class="special">*</span> <span class="identifier">m_N</span> <span class="special">)</span> <span class="special">)</span> <span class="special">,</span> 570 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_zip_iterator</span><span class="special">(</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span> 571 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">)</span> <span class="special">+</span> <span class="identifier">m_N</span> <span class="special">,</span> 572 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">)</span> <span class="special">+</span> <span class="number">2</span> <span class="special">*</span> <span class="identifier">m_N</span> <span class="special">,</span> 573 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">)</span> <span class="special">+</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">m_N</span> <span class="special">,</span> 574 <span class="identifier">m_beta</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> 575 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">dxdt</span> <span class="special">)</span> <span class="special">+</span> <span class="identifier">m_N</span> <span class="special">,</span> 576 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">dxdt</span> <span class="special">)</span> <span class="special">+</span> <span class="number">2</span> <span class="special">*</span> <span class="identifier">m_N</span> <span class="special">,</span> 577 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">dxdt</span> <span class="special">)</span> <span class="special">+</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">m_N</span> <span class="special">)</span> <span class="special">)</span> <span class="special">,</span> 578 <span class="identifier">lorenz_functor</span><span class="special">()</span> <span class="special">);</span> 579 <span class="special">}</span> 580 <span class="identifier">size_t</span> <span class="identifier">m_N</span><span class="special">;</span> 581 <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">m_beta</span><span class="special">;</span> 582<span class="special">};</span> 583</pre> 584<p> 585 </p> 586<p> 587 As <code class="computeroutput"><span class="identifier">state_type</span></code> a <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">device_vector</span></code> or a <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a> 588 of a <code class="computeroutput"><span class="identifier">device_vector</span></code> is used. 589 The length of the state is <span class="emphasis"><em>3N</em></span> where <span class="emphasis"><em>N</em></span> 590 is the number of systems. The system is encoded into this vector such that 591 all <span class="emphasis"><em>x</em></span> components come first, then every <span class="emphasis"><em>y</em></span> 592 components and finally every <span class="emphasis"><em>z</em></span> components. Implementing 593 the device function is then a simple task, you only have to decompose the 594 tuple originating from the zip iterators. 595 </p> 596<p> 597 Besides the system without perturbations we furthermore need to calculate 598 the system including linearized equations governing the time evolution 599 of small perturbations. Using the method from above this is straightforward, 600 with a small difficulty that Thrust's tuples have a maximal arity of 10. 601 But this is only a small problem since we can create a zip iterator packed 602 with zip iterators. So the top level zip iterator contains one zip iterator 603 for the state, one normal iterator for the parameter, and one zip iterator 604 for the derivative. Accessing the elements of this tuple in the system 605 function is then straightforward, you unpack the tuple with <code class="computeroutput"><span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><>()</span></code>. 606 We will not show the code here, it is to large. It can be found <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/thrust/lorenz_parameters.cu" target="_top">here</a> and 607 is easy to understand. 608 </p> 609<p> 610 Furthermore, we need an observer which determines the norm of the perturbations, 611 normalizes them and averages the logarithm of the norm. The device functor 612 which is used within this observer is defined 613 </p> 614<p> 615</p> 616<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">lyap_functor</span> 617<span class="special">{</span> 618 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">T</span> <span class="special">></span> 619 <span class="identifier">__host__</span> <span class="identifier">__device__</span> 620 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="identifier">T</span> <span class="identifier">t</span> <span class="special">)</span> <span class="keyword">const</span> 621 <span class="special">{</span> 622 <span class="identifier">value_type</span> <span class="special">&</span><span class="identifier">dx</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">0</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">);</span> 623 <span class="identifier">value_type</span> <span class="special">&</span><span class="identifier">dy</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">1</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">);</span> 624 <span class="identifier">value_type</span> <span class="special">&</span><span class="identifier">dz</span> <span class="special">=</span> <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">2</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">);</span> 625 <span class="identifier">value_type</span> <span class="identifier">norm</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span> <span class="identifier">dx</span> <span class="special">*</span> <span class="identifier">dx</span> <span class="special">+</span> <span class="identifier">dy</span> <span class="special">*</span> <span class="identifier">dy</span> <span class="special">+</span> <span class="identifier">dz</span> <span class="special">*</span> <span class="identifier">dz</span> <span class="special">);</span> 626 <span class="identifier">dx</span> <span class="special">/=</span> <span class="identifier">norm</span><span class="special">;</span> 627 <span class="identifier">dy</span> <span class="special">/=</span> <span class="identifier">norm</span><span class="special">;</span> 628 <span class="identifier">dz</span> <span class="special">/=</span> <span class="identifier">norm</span><span class="special">;</span> 629 <span class="identifier">thrust</span><span class="special">::</span><span class="identifier">get</span><span class="special"><</span> <span class="number">3</span> <span class="special">>(</span> <span class="identifier">t</span> <span class="special">)</span> <span class="special">+=</span> <span class="identifier">log</span><span class="special">(</span> <span class="identifier">norm</span> <span class="special">);</span> 630 <span class="special">}</span> 631<span class="special">};</span> 632</pre> 633<p> 634 </p> 635<p> 636 Note, that this functor manipulates the state, i.e. the perturbations. 637 </p> 638<p> 639 Now we complete the whole code to calculate the Lyapunov exponents. First, 640 we have to define a state vector. This vector contains <span class="emphasis"><em>6N</em></span> 641 entries, the state <span class="emphasis"><em>x,y,z</em></span> and its perturbations <span class="emphasis"><em>dx,dy,dz</em></span>. 642 We initialize them such that <span class="emphasis"><em>x=y=z=10</em></span>, <span class="emphasis"><em>dx=1</em></span>, 643 and <span class="emphasis"><em>dy=dz=0</em></span>. We define a stepper type, a controlled 644 Runge-Kutta Dormand-Prince 5 stepper. We start with some integration to 645 overcome the transient behavior. For this, we do not involve the perturbation 646 and run the algorithm only on the state <span class="emphasis"><em>x,y,z</em></span> without 647 any observer. Note, how <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a> 648 is used for partial integration of the state vector without perturbations 649 (the first half of the whole state). After the transient, the full system 650 with perturbations is integrated and the Lyapunov exponents are calculated 651 and written to <code class="computeroutput"><span class="identifier">stdout</span></code>. 652 </p> 653<p> 654</p> 655<pre class="programlisting"><span class="identifier">state_type</span> <span class="identifier">x</span><span class="special">(</span> <span class="number">6</span> <span class="special">*</span> <span class="identifier">N</span> <span class="special">);</span> 656 657<span class="comment">// initialize x,y,z</span> 658<span class="identifier">thrust</span><span class="special">::</span><span class="identifier">fill</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">N</span> <span class="special">,</span> <span class="number">10.0</span> <span class="special">);</span> 659 660<span class="comment">// initial dx</span> 661<span class="identifier">thrust</span><span class="special">::</span><span class="identifier">fill</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">N</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="number">4</span> <span class="special">*</span> <span class="identifier">N</span> <span class="special">,</span> <span class="number">1.0</span> <span class="special">);</span> 662 663<span class="comment">// initialize dy,dz</span> 664<span class="identifier">thrust</span><span class="special">::</span><span class="identifier">fill</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="number">4</span> <span class="special">*</span> <span class="identifier">N</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">);</span> 665 666 667<span class="comment">// create error stepper, can be used with make_controlled or make_dense_output</span> 668<span class="keyword">typedef</span> <span class="identifier">runge_kutta_dopri5</span><span class="special"><</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">stepper_type</span><span class="special">;</span> 669 670 671<span class="identifier">lorenz_system</span> <span class="identifier">lorenz</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">,</span> <span class="identifier">beta</span> <span class="special">);</span> 672<span class="identifier">lorenz_perturbation_system</span> <span class="identifier">lorenz_perturbation</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">,</span> <span class="identifier">beta</span> <span class="special">);</span> 673<span class="identifier">lyap_observer</span> <span class="identifier">obs</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">,</span> <span class="number">1</span> <span class="special">);</span> 674 675<span class="comment">// calculate transients</span> 676<span class="identifier">integrate_adaptive</span><span class="special">(</span> <span class="identifier">make_controlled</span><span class="special">(</span> <span class="number">1.0e-6</span> <span class="special">,</span> <span class="number">1.0e-6</span> <span class="special">,</span> <span class="identifier">stepper_type</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">lorenz</span> <span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">make_pair</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">N</span> <span class="special">)</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">10.0</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">);</span> 677 678<span class="comment">// calculate the Lyapunov exponents -- the main loop</span> 679<span class="keyword">double</span> <span class="identifier">t</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span> 680<span class="keyword">while</span><span class="special">(</span> <span class="identifier">t</span> <span class="special"><</span> <span class="number">10000.0</span> <span class="special">)</span> 681<span class="special">{</span> 682 <span class="identifier">integrate_adaptive</span><span class="special">(</span> <span class="identifier">make_controlled</span><span class="special">(</span> <span class="number">1.0e-6</span> <span class="special">,</span> <span class="number">1.0e-6</span> <span class="special">,</span> <span class="identifier">stepper_type</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">lorenz_perturbation</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">+</span> <span class="number">1.0</span> <span class="special">,</span> <span class="number">0.1</span> <span class="special">);</span> 683 <span class="identifier">t</span> <span class="special">+=</span> <span class="number">1.0</span><span class="special">;</span> 684 <span class="identifier">obs</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">);</span> 685<span class="special">}</span> 686 687<span class="identifier">vector</span><span class="special"><</span> <span class="identifier">value_type</span> <span class="special">></span> <span class="identifier">lyap</span><span class="special">(</span> <span class="identifier">N</span> <span class="special">);</span> 688<span class="identifier">obs</span><span class="special">.</span><span class="identifier">fill_lyap</span><span class="special">(</span> <span class="identifier">lyap</span> <span class="special">);</span> 689 690<span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">i</span><span class="special"><</span><span class="identifier">N</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span> 691 <span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">beta_host</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special"><<</span> <span class="string">"\t"</span> <span class="special"><<</span> <span class="identifier">lyap</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special"><<</span> <span class="string">"\n"</span><span class="special">;</span> 692</pre> 693<p> 694 </p> 695<p> 696 The full example can be found at <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/thrust/lorenz_parameters.cu" target="_top">lorenz_parameters.cu</a>. 697 </p> 698</div> 699</div> 700<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 701<td align="left"></td> 702<td align="right"><div class="copyright-footer">Copyright © 2009-2015 Karsten Ahnert and Mario Mulansky<p> 703 Distributed under the Boost Software License, Version 1.0. (See accompanying 704 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 705 </p> 706</div></td> 707</tr></table> 708<hr> 709<div class="spirit-nav"> 710<a accesskey="p" href="self_expanding_lattices.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="using_opencl_via_vexcl.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a> 711</div> 712</body> 713</html> 714
