1<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" 2 "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd"> 3<html xmlns="http://www.w3.org/1999/xhtml"> 4<head> 5<meta name="generator" content= 6"HTML Tidy for Linux/x86 (vers 1st March 2004), see www.w3.org" /> 7<meta name="GENERATOR" content="Quanta Plus" /> 8<meta http-equiv="Content-Type" content= 9"text/html; charset=us-ascii" /> 10<link rel="stylesheet" href="../../../../boost.css" type="text/css"/> 11<link rel="stylesheet" href="ublas.css" type="text/css" /> 12<script type="text/javascript" src="js/jquery-1.3.2.min.js" async="async" ></script> 13<script type="text/javascript" src="js/jquery.toc-gw.js" async="async" ></script> 14<title>uBLAS operations overview</title> 15</head> 16<body> 17<h1><img src="../../../../boost.png" align="middle" />Overview of Tensor, Matrix and Vector Operations</h1> 18<div class="toc" id="toc"></div> 19 20<dl> 21<dt>Contents:</dt> 22<dd><a href="#blas">Basic Linear Algebra</a></dd> 23<dd><a href="#advanced">Advanced Functions</a></dd> 24<dd><a href="#sub">Submatrices, Subvectors</a></dd> 25<dd><a href="#speed">Speed Improvements</a></dd> 26</dl> 27 28<h2>Definitions</h2> 29 30<table style="" summary="notation"> 31<tr><td><code>X, Y, Z</code></td> 32<td> are tensors</td></tr> 33<tr><td><code>A, B, C</code></td> 34<td> are matrices</td></tr> 35<tr><td><code>u, v, w</code></td> 36<td>are vectors</td></tr> 37<tr><td><code>i, j, k</code></td> 38<td>are integer values</td></tr> 39<tr><td><code>t, t1, t2</code></td> 40<td>are scalar values</td></tr> 41<tr><td><code>r, r1, r2</code></td> 42<td>are <a href="range.html">ranges</a>, e.g. <code>range(0, 3)</code></td></tr> 43<tr><td><code>s, s1, s2</code></td> 44<td>are <a href="range.html#slice">slices</a>, e.g. <code>slice(0, 1, 3)</code></td></tr> 45</table> 46 47<h2><a name="blas">Basic Linear Algebra</a></h2> 48 49<h3>standard operations: addition, subtraction, multiplication by a 50scalar</h3> 51<pre><code> 52X = Y + Z; X = Y - Z; X = -Y; 53C = A + B; C = A - B; C = -A; 54w = u + v; w = u - v; w = -u; 55X = t * Y; Y = X * t; X = Y / t; 56C = t * A; C = A * t; C = A / t; 57w = t * u; w = u * t; w = u / t; 58</code></pre> 59 60<h3>computed assignments</h3> 61<pre><code> 62X += Y; X -= Y; 63C += A; C -= A; 64w += u; w -= u; 65X *= t; X /= t; 66C *= t; C /= t; 67w *= t; w /= t; 68</code></pre> 69 70<h3>inner, outer and other products</h3> 71<pre><code> 72t = inner_prod(u, v); 73C = outer_prod(u, v); 74w = prod(A, u); w = prod(u, A); w = prec_prod(A, u); w = prec_prod(u, A); 75C = prod(A, B); C = prec_prod(A, B); 76w = element_prod(u, v); w = element_div(u, v); 77C = element_prod(A, B); C = element_div(A, B); 78</code></pre> 79 80<h3>tensor products</h3> 81<pre><code> 82Z = prod(X, v, t); 83Z = prod(X, A, t); 84Z = prod(X, Y, p); 85Z = prod(X, Y, pa, pb); 86t = inner_prod(X, Y); 87Z = outer_prod(X, Y); 88</code></pre> 89 90<h3>transformations</h3> 91<pre><code> 92w = conj(u); w = real(u); w = imag(u); 93C = trans(A); C = conj(A); C = herm(A); C = real(A); C = imag(A); 94Z = trans(X); Z = conj(X); Z = real(X); Z = imag(X); 95</code></pre> 96 97 98<h2><a name="advanced">Advanced functions</a></h2> 99 100<h3>norms</h3> 101 102<pre><code> 103t = norm_inf(v); i = index_norm_inf(v); 104t = norm_1(v); t = norm_2(v); 105t = norm_2_square(v); 106t = norm_inf(A); i = index_norm_inf(A); 107t = norm_1(A); t = norm_frobenius(A); 108t = norm(X); 109</code></pre> 110 111<h3>products</h3> 112 113<pre><code> 114axpy_prod(A, u, w, true); // w = A * u 115axpy_prod(A, u, w, false); // w += A * u 116axpy_prod(u, A, w, true); // w = trans(A) * u 117axpy_prod(u, A, w, false); // w += trans(A) * u 118axpy_prod(A, B, C, true); // C = A * B 119axpy_prod(A, B, C, false); // C += A * B 120</code></pre> 121<p><em>Note:</em> The last argument (<code>bool init</code>) of 122<code>axpy_prod</code> is optional. Currently it defaults to 123<code>true</code>, but this may change in the future. Setting the 124<code>init</code> to <code>true</code> is equivalent to calling 125<code>w.clear()</code> before <code>axpy_prod</code>. 126There are some specialisation for products of compressed matrices that give a 127large speed up compared to <code>prod</code>.</p> 128<pre><code> 129w = block_prod<matrix_type, 64> (A, u); // w = A * u 130w = block_prod<matrix_type, 64> (u, A); // w = trans(A) * u 131C = block_prod<matrix_type, 64> (A, B); // C = A * B 132</code></pre> 133<p><em>Note:</em> The blocksize can be any integer. However, the 134actual speed depends very significantly on the combination of blocksize, 135CPU and compiler. The function <code>block_prod</code> is designed 136for large dense matrices.</p> 137<h3>rank-k updates</h3> 138<pre><code> 139opb_prod(A, B, C, true); // C = A * B 140opb_prod(A, B, C, false); // C += A * B 141</code></pre> 142<p><em>Note:</em> The last argument (<code>bool init</code>) of 143<code>opb_prod</code> is optional. Currently it defaults to 144<code>true</code>, but this may change in the future. This function 145may give a speedup if <code>A</code> has less columns than rows, 146because the product is computed as a sum of outer products.</p> 147 148<h2><a name="sub">Submatrices, Subvectors</a></h2> 149<p>Accessing submatrices and subvectors via <b>proxies</b> using <code>project</code> functions:</p> 150<pre><code> 151w = project(u, r); // the subvector of u specifed by the index range r 152w = project(u, s); // the subvector of u specifed by the index slice s 153C = project(A, r1, r2); // the submatrix of A specified by the two index ranges r1 and r2 154C = project(A, s1, s2); // the submatrix of A specified by the two index slices s1 and s2 155w = row(A, i); w = column(A, j); // a row or column of matrix as a vector 156</code></pre> 157<p>Assigning to submatrices and subvectors via <b>proxies</b> using <code>project</code> functions:</p> 158<pre><code> 159project(u, r) = w; // assign the subvector of u specifed by the index range r 160project(u, s) = w; // assign the subvector of u specifed by the index slice s 161project(A, r1, r2) = C; // assign the submatrix of A specified by the two index ranges r1 and r2 162project(A, s1, s2) = C; // assign the submatrix of A specified by the two index slices s1 and s2 163row(A, i) = w; column(A, j) = w; // a row or column of matrix as a vector 164</code></pre> 165<p><em>Note:</em> A range <code>r = range(start, stop)</code> 166contains all indices <code>i</code> with <code>start <= i < 167stop</code>. A slice is something more general. The slice 168<code>s = slice(start, stride, size)</code> contains the indices 169<code>start, start+stride, ..., start+(size-1)*stride</code>. The 170stride can be 0 or negative! If <code>start >= stop</code> for a range 171or <code>size == 0</code> for a slice then it contains no elements.</p> 172<p>Sub-ranges and sub-slices of vectors and matrices can be created directly with the <code>subrange</code> and <code>sublice</code> functions:</p> 173<pre><code> 174w = subrange(u, 0, 2); // the 2 element subvector of u 175w = subslice(u, 0, 1, 2); // the 2 element subvector of u 176C = subrange(A, 0,2, 0,3); // the 2x3 element submatrix of A 177C = subslice(A, 0,1,2, 0,1,3); // the 2x3 element submatrix of A 178subrange(u, 0, 2) = w; // assign the 2 element subvector of u 179subslice(u, 0, 1, 2) = w; // assign the 2 element subvector of u 180subrange(A, 0,2, 0,3) = C; // assign the 2x3 element submatrix of A 181subrange(A, 0,1,2, 0,1,3) = C; // assigne the 2x3 element submatrix of A 182</code></pre> 183<p>There are to more ways to access some matrix elements as a 184vector:</p> 185<pre><code>matrix_vector_range<matrix_type> (A, r1, r2); 186matrix_vector_slice<matrix_type> (A, s1, s2); 187</code></pre> 188<p><em>Note:</em> These matrix proxies take a sequence of elements 189of a matrix and allow you to access these as a vector. In 190particular <code>matrix_vector_slice</code> can do this in a very 191general way. <code>matrix_vector_range</code> is less useful as the 192elements must lie along a diagonal.</p> 193<p><em>Example:</em> To access the first two elements of a sub 194column of a matrix we access the row with a slice with stride 1 and 195the column with a slice with stride 0 thus:<br /> 196<code>matrix_vector_slice<matrix_type> (A, slice(0,1,2), 197slice(0,0,2)); 198</code></p> 199 200<h2><a name="speed">Speed improvements</a></h2> 201<h3><a name='noalias'>Matrix / Vector assignment</a></h3> 202<p>If you know for sure that the left hand expression and the right 203hand expression have no common storage, then assignment has 204no <em>aliasing</em>. A more efficient assignment can be specified 205in this case:</p> 206<pre><code>noalias(C) = prod(A, B); 207</code></pre> 208<p>This avoids the creation of a temporary matrix that is required in a normal assignment. 209'noalias' assignment requires that the left and right hand side be size conformant.</p> 210 211<h3>Sparse element access</h3> 212<p>The matrix element access function <code>A(i1,i2)</code> or the equivalent vector 213element access functions (<code>v(i) or v[i]</code>) usually create 'sparse element proxies' 214when applied to a sparse matrix or vector. These <em>proxies</em> allow access to elements 215without having to worry about nasty C++ issues where references are invalidated.</p> 216<p>These 'sparse element proxies' can be implemented more efficiently when applied to <code>const</code> 217objects. 218Sadly in C++ there is no way to distinguish between an element access on the left and right hand side of 219an assignment. Most often elements on the right hand side will not be changed and therefore it would 220be better to use the <code>const</code> proxies. We can do this by making the matrix or vector 221<code>const</code> before accessing it's elements. For example:</p> 222<pre><code>value = const_cast<const VEC>(v)[i]; // VEC is the type of V 223</code></pre> 224<p>If more then one element needs to be accessed <code>const_iterator</code>'s should be used 225in preference to <code>iterator</code>'s for the same reason. For the more daring 'sparse element proxies' 226can be completely turned off in uBLAS by defining the configuration macro <code>BOOST_UBLAS_NO_ELEMENT_PROXIES</code>. 227</p> 228 229 230<h3>Controlling the complexity of nested products</h3> 231 232<p>What is the complexity (the number of add and multiply operations) required to compute the following? 233</p> 234<pre> 235 R = prod(A, prod(B,C)); 236</pre> 237<p>Firstly the complexity depends on matrix size. Also since prod is transitive (not commutative) 238the bracket order affects the complexity. 239</p> 240<p>uBLAS evaluates expressions without matrix or vector temporaries and honours 241the bracketing structure. However avoiding temporaries for nested product unnecessarly increases the complexity. 242Conversly by explictly using temporary matrices the complexity of a nested product can be reduced. 243</p> 244<p>uBLAS provides 3 alternative syntaxes for this purpose: 245</p> 246<pre> 247 temp_type T = prod(B,C); R = prod(A,T); // Preferable if T is preallocated 248</pre> 249<pre> 250 prod(A, temp_type(prod(B,C)); 251</pre> 252<pre> 253 prod(A, prod<temp_type>(B,C)); 254</pre> 255<p>The 'temp_type' is important. Given A,B,C are all of the same type. Say 256matrix<float>, the choice is easy. However if the value_type is mixed (int with float or double) 257or the matrix type is mixed (sparse with symmetric) the best solution is not so obvious. It is up to you! It 258depends on numerical properties of A and the result of the prod(B,C). 259</p> 260 261<hr /> 262<p>Copyright (©) 2000-2007 Joerg Walter, Mathias Koch, Gunter 263Winkler, Michael Stevens<br /> 264 Use, modification and distribution are subject to the 265 Boost Software License, Version 1.0. 266 (See accompanying file LICENSE_1_0.txt 267 or copy at <a href="http://www.boost.org/LICENSE_1_0.txt"> 268 http://www.boost.org/LICENSE_1_0.txt 269 </a>). 270</p> 271<script type="text/javascript"> 272(function($) { 273 $('#toc').toc(); 274})(jQuery); 275</script> 276</body> 277</html> 278
