ceph/src/boost/libs/numeric/ublas/doc/operations_overview.html

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