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 href="ublas.css" type="text/css" /> |
---|
11 | <title>uBLAS operations overview</title> |
---|
12 | </head> |
---|
13 | <body> |
---|
14 | <h1><img src="../../../../boost.png" align="middle" /> |
---|
15 | Overview of Matrix and Vector Operations</h1> |
---|
16 | |
---|
17 | <dl> |
---|
18 | <dt>Contents:</dt> |
---|
19 | <dd><a href="#blas">Basic Linear Algebra</a></dd> |
---|
20 | <dd><a href="#advanced">Advanced Functions</a></dd> |
---|
21 | <dd><a href="#sub">Submatrices, Subvectors</a></dd> |
---|
22 | <dd><a href="#speed">Speed Improvements</a></dd> |
---|
23 | </dl> |
---|
24 | |
---|
25 | <h3>Definitions:</h3> |
---|
26 | |
---|
27 | <table style="" summary="notation"> |
---|
28 | <tr><td><code>A, B, C</code></td> |
---|
29 | <td> are matrices</td></tr> |
---|
30 | <tr><td><code>u, v, w</code></td> |
---|
31 | <td>are vectors</td></tr> |
---|
32 | <tr><td><code>i, j, k</code></td> |
---|
33 | <td>are integer values</td></tr> |
---|
34 | <tr><td><code>t, t1, t2</code></td> |
---|
35 | <td>are scalar values</td></tr> |
---|
36 | <tr><td><code>r, r1, r2</code></td> |
---|
37 | <td>are <a href="range.htm">ranges</a>, e.g. <code>range(0, 3)</code></td></tr> |
---|
38 | <tr><td><code>s, s1, s2</code></td> |
---|
39 | <td>are <a href="range.htm#slice">slices</a>, e.g. <code>slice(0, 1, 3)</code></td></tr> |
---|
40 | </table> |
---|
41 | |
---|
42 | <h2><a name="blas">Basic Linear Algebra</a></h2> |
---|
43 | |
---|
44 | <h3>standard operations: addition, subtraction, multiplication by a |
---|
45 | scalar</h3> |
---|
46 | |
---|
47 | <pre><code> |
---|
48 | C = A + B; C = A - B; C = -A; |
---|
49 | w = u + v; w = u - v; w = -u; |
---|
50 | C = t * A; C = A * t; C = A / t; |
---|
51 | w = t * u; w = u * t; w = u / t; |
---|
52 | </code></pre> |
---|
53 | |
---|
54 | <h3>computed assignements</h3> |
---|
55 | |
---|
56 | <pre><code> |
---|
57 | C += A; C -= A; |
---|
58 | w += u; w -= u; |
---|
59 | C *= t; C /= t; |
---|
60 | w *= t; w /= t; |
---|
61 | </code></pre> |
---|
62 | |
---|
63 | <h3>inner, outer and other products</h3> |
---|
64 | |
---|
65 | <pre><code> |
---|
66 | t = inner_prod(u, v); |
---|
67 | C = outer_prod(u, v); |
---|
68 | w = prod(A, u); w = prod(u, A); w = prec_prod(A, u); w = prec_prod(u, A); |
---|
69 | C = prod(A, B); C = prec_prod(A, B); |
---|
70 | w = element_prod(u, v); w = element_div(u, v); |
---|
71 | C = element_prod(A, B); C = element_div(A, B); |
---|
72 | </code></pre> |
---|
73 | |
---|
74 | <h3>transformations</h3> |
---|
75 | |
---|
76 | <pre><code> |
---|
77 | w = conj(u); w = real(u); w = imag(u); |
---|
78 | C = trans(A); C = conj(A); C = herm(A); C = real(A); C = imag(A); |
---|
79 | </code></pre> |
---|
80 | |
---|
81 | <h2><a name="advanced">Advanced functions</a></h2> |
---|
82 | |
---|
83 | <h3>norms</h3> |
---|
84 | |
---|
85 | <pre><code> |
---|
86 | t = norm_inf(v); i = index_norm_inf(v); |
---|
87 | t = norm_1(v); t = norm_2(v); |
---|
88 | t = norm_inf(A); i = index_norm_inf(A); |
---|
89 | t = norm_1(A); t = norm_frobenius(A); |
---|
90 | </code></pre> |
---|
91 | |
---|
92 | <h3>products</h3> |
---|
93 | |
---|
94 | <pre><code> |
---|
95 | axpy_prod(A, u, w, true); // w = A * u |
---|
96 | axpy_prod(A, u, w, false); // w += A * u |
---|
97 | axpy_prod(u, A, w, true); // w = trans(A) * u |
---|
98 | axpy_prod(u, A, w, false); // w += trans(A) * u |
---|
99 | axpy_prod(A, B, C, true); // C = A * B |
---|
100 | axpy_prod(A, B, C, false); // C += A * B |
---|
101 | </code></pre> |
---|
102 | <p><em>Note:</em> The last argument (<code>bool init</code>) of |
---|
103 | <code>axpy_prod</code> is optional. Currently it defaults to |
---|
104 | <code>true</code>, but this may change in the future. Set the |
---|
105 | <code>init</code> to <code>true</code> is equivalent to call |
---|
106 | <code>w.clear()</code> before <code>axpy_prod</code>. Up to now |
---|
107 | there are some specialisation for compressed matrices that give a |
---|
108 | large speed up compared to <code>prod</code>.</p> |
---|
109 | <pre><code> |
---|
110 | w = block_prod<matrix_type, 64> (A, u); // w = A * u |
---|
111 | w = block_prod<matrix_type, 64> (u, A); // w = trans(A) * u |
---|
112 | C = block_prod<matrix_type, 64> (A, B); // w = A * B |
---|
113 | </code></pre> |
---|
114 | <p><em>Note:</em> The blocksize can be any integer. However, the |
---|
115 | total speed depends very strong on the combination of blocksize, |
---|
116 | CPU and compiler. The function <code>block_prod</code> is designed |
---|
117 | for large dense matrices.</p> |
---|
118 | <h3>rank-k updates</h3> |
---|
119 | <pre><code> |
---|
120 | opb_prod(A, B, C, true); // C = A * B |
---|
121 | opb_prod(A, B, C, false); // C += A * B |
---|
122 | </code></pre> |
---|
123 | <p><em>Note:</em> The last argument (<code>bool init</code>) of |
---|
124 | <code>opb_prod</code> is optional. Currently it defaults to |
---|
125 | <code>true</code>, but this may change in the future. This function |
---|
126 | may give a speedup if <code>A</code> has less columns than rows, |
---|
127 | because the product is computed as a sum of outer products.</p> |
---|
128 | |
---|
129 | <h2><a name="sub">Submatrices, Subvectors</a></h2> |
---|
130 | <p>Accessing submatrices and subvectors via <b>proxies</b> using <code>project</code> functions:</p> |
---|
131 | <pre><code> |
---|
132 | w = project(u, r); // the subvector of u specifed by the index range r |
---|
133 | w = project(u, s); // the subvector of u specifed by the index slice s |
---|
134 | C = project(A, r1, r2); // the submatrix of A specified by the two index ranges r1 and r2 |
---|
135 | C = project(A, s1, s2); // the submatrix of A specified by the two index slices s1 and s2 |
---|
136 | w = row(A, i); w = column(A, j); // a row or column of matrix as a vector |
---|
137 | </code></pre> |
---|
138 | <p>Assigning to submatrices and subvectors via <b>proxies</b> using <code>project</code> functions:</p> |
---|
139 | <pre><code> |
---|
140 | project(u, r) = w; // assign the subvector of u specifed by the index range r |
---|
141 | project(u, s) = w; // assign the subvector of u specifed by the index slice s |
---|
142 | project(A, r1, r2) = C; // assign the submatrix of A specified by the two index ranges r1 and r2 |
---|
143 | project(A, s1, s2) = C; // assign the submatrix of A specified by the two index slices s1 and s2 |
---|
144 | row(A, i) = w; column(A, j) = w; // a row or column of matrix as a vector |
---|
145 | </code></pre> |
---|
146 | <p><em>Note:</em> A range <code>r = range(start, stop)</code> |
---|
147 | contains all indices <code>i</code> with <code>start <= i < |
---|
148 | stop</code>. A slice is something more general. The slice |
---|
149 | <code>s = slice(start, stride, size)</code> contains the indices |
---|
150 | <code>start, start+stride, ..., start+(size-1)*stride</code>. The |
---|
151 | stride can be 0 or negative! If <code>start >= stop</code> for a range |
---|
152 | or <code>size == 0</code> for a slice then it contains no elements.</p> |
---|
153 | <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> |
---|
154 | <pre><code> |
---|
155 | w = subrange(u, 0, 2); // the 2 element subvector of u |
---|
156 | w = subslice(u, 0, 1, 2); // the 2 element subvector of u |
---|
157 | C = subrange(A, 0,2, 0,3); // the 2x3 element submatrix of A |
---|
158 | C = subslice(A, 0,1,2, 0,1,3); // the 2x3 element submatrix of A |
---|
159 | subrange(u, 0, 2) = w; // assign the 2 element subvector of u |
---|
160 | subslice(u, 0, 1, 2) = w; // assign the 2 element subvector of u |
---|
161 | subrange(A, 0,2, 0,3) = C; // assign the 2x3 element submatrix of A |
---|
162 | subrange(A, 0,1,2, 0,1,3) = C; // assigne the 2x3 element submatrix of A |
---|
163 | </code></pre> |
---|
164 | <p>There are to more ways to access some matrix elements as a |
---|
165 | vector:</p> |
---|
166 | <pre><code>matrix_vector_range<matrix_type> (A, r1, r2); |
---|
167 | matrix_vector_slice<matrix_type> (A, s1, s2); |
---|
168 | </code></pre> |
---|
169 | <p><em>Note:</em> These matrix proxies take a sequence of elements |
---|
170 | of a matrix and allow you to access these as a vector. In |
---|
171 | particular <code>matrix_vector_slice</code> can do this in a very |
---|
172 | general way. <code>matrix_vector_range</code> is less useful as the |
---|
173 | elements must lie along a diagonal.</p> |
---|
174 | <p><em>Example:</em> To access the first two elements of a sub |
---|
175 | column of a matrix we access the row with a slice with stride 1 and |
---|
176 | the column with a slice with stride 0 thus:<br /> |
---|
177 | <code>matrix_vector_slice<matrix_type> (A, slice(0,1,2), |
---|
178 | slice(0,0,2)); |
---|
179 | </code></p> |
---|
180 | |
---|
181 | <h2><a name="speed">Speed improvements</a></h2> |
---|
182 | <h3><a name='noalias'>Matrix / Vector assignment</a></h3> |
---|
183 | <p>If you know for sure that the left hand expression and the right |
---|
184 | hand expression have no common storage, then assignment has |
---|
185 | no <em>aliasing</em>. A more efficient assignment can be specified |
---|
186 | in this case:</p> |
---|
187 | <pre><code>noalias(C) = prod(A, B); |
---|
188 | </code></pre> |
---|
189 | <p>This avoids the creation of a temporary matrix that is required in a normal assignment. |
---|
190 | 'noalias' assignment requires that the left and right hand side be size conformant.</p> |
---|
191 | |
---|
192 | <h3>Sparse element access</h3> |
---|
193 | <p>The matrix element access function <code>A(i1,i2)</code> or the equivalent vector |
---|
194 | element access functions (<code>v(i) or v[i]</code>) usually create 'sparse element proxies' |
---|
195 | when applied to a sparse matrix or vector. These <em>proxies</em> allow access to elements |
---|
196 | without having to worry about nasty C++ issues where references are invalidated.</p> |
---|
197 | <p>These 'sparse element proxies' can be implemented more efficiently when applied to <code>const</code> |
---|
198 | objects. |
---|
199 | Sadly in C++ there is no way to distinguish between an element access on the left and right hand side of |
---|
200 | an assignment. Most often elements on the right hand side will not be changed and therefore it would |
---|
201 | be better to use the <code>const</code> proxies. We can do this by making the matrix or vector |
---|
202 | <code>const</code> before accessing it's elements. For example:</p> |
---|
203 | <pre><code>value = const_cast<const VEC>(v)[i]; // VEC is the type of V |
---|
204 | </code></pre> |
---|
205 | <p>If more then one element needs to be accessed <code>const_iterator</code>'s should be used |
---|
206 | in preference to <code>iterator</code>'s for the same reason. For the more daring 'sparse element proxies' |
---|
207 | can be completely turned off in uBLAS by defining the configuration macro <code>BOOST_UBLAS_NO_ELEMENT_PROXIES</code>. |
---|
208 | </p> |
---|
209 | |
---|
210 | <hr /> |
---|
211 | <p>Copyright (©) 2000-2004 Joerg Walter, Mathias Koch, Gunter |
---|
212 | Winkler, Michael Stevens<br /> |
---|
213 | Permission to copy, use, modify, sell and distribute this document |
---|
214 | is granted provided this copyright notice appears in all copies. |
---|
215 | This document is provided ``as is'' without express or implied |
---|
216 | warranty, and with no claim as to its suitability for any |
---|
217 | purpose.</p> |
---|
218 | </body> |
---|
219 | </html> |
---|