]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // |
2 | // Copyright (c) 2010 Athanasios Iliopoulos | |
3 | // | |
4 | // Distributed under the Boost Software License, Version 1.0. (See | |
5 | // accompanying file LICENSE_1_0.txt or copy at | |
6 | // http://www.boost.org/LICENSE_1_0.txt) | |
7 | // | |
8 | ||
9 | #include <boost/numeric/ublas/assignment.hpp> | |
10 | #include <boost/numeric/ublas/vector.hpp> | |
11 | #include <boost/numeric/ublas/vector_proxy.hpp> | |
12 | #include <boost/numeric/ublas/matrix_proxy.hpp> | |
13 | #include <boost/numeric/ublas/vector_sparse.hpp> | |
14 | #include <boost/numeric/ublas/matrix_sparse.hpp> | |
15 | #include <boost/numeric/ublas/io.hpp> | |
16 | #include <boost/numeric/ublas/matrix.hpp> | |
17 | ||
18 | using namespace boost::numeric::ublas; | |
19 | ||
20 | int main() { | |
21 | // Simple vector fill | |
22 | vector<double> a(3); | |
23 | a <<= 0, 1, 2; | |
24 | std::cout << a << std::endl; | |
25 | // [ 0 1 2] | |
26 | ||
27 | // Vector from vector | |
28 | vector<double> b(7); | |
29 | b <<= a, 10, a; | |
30 | std::cout << b << std::endl; | |
31 | // [ 0 1 2 10 0 1 2] | |
32 | ||
33 | // Simple matrix fill | |
34 | matrix<double> A(3,3); | |
35 | A <<= 0, 1, 2, | |
36 | 3, 4, 5, | |
37 | 6, 7, 8; | |
38 | std::cout << A << std::endl; | |
39 | // [ 0 1 2 ] | |
40 | // [ 3 4 5 ] | |
41 | // [ 6 7 8 ] | |
42 | ||
43 | // Matrix from vector | |
44 | A <<= 0, 1, 2, | |
45 | 3, 4, 5, | |
46 | a; | |
47 | std::cout << A << std::endl; | |
48 | // [ 0 1 2 ] | |
49 | // [ 3 4 5 ] | |
50 | // [ 0 1 2 ] | |
51 | ||
52 | // Matrix from vector - column assignment | |
53 | A <<= move(0,2), traverse_policy::by_column(), | |
54 | a; | |
55 | std::cout << A << std::endl; | |
56 | // [ 0 1 0 ] | |
57 | // [ 3 4 1 ] | |
58 | // [ 0 1 2 ] | |
59 | ||
60 | // Another matrix from vector example (watch the wraping); | |
61 | vector<double> c(9); c <<= 1, 2, 3, 4, 5, 6, 7, 8, 9; | |
62 | A <<= c; | |
63 | std::cout << A << std::endl; | |
64 | // [ 1 2 3 ] | |
65 | // [ 4 5 6 ] | |
66 | // [ 7 8 9 ] | |
67 | ||
68 | // If for performance(Benchmarks are not definite about that) or consistency reasons you need to disable wraping: | |
69 | static next_row_manip endr; //This can be defined globally | |
70 | A <<= traverse_policy::by_row_no_wrap(), | |
71 | 1, 2, 3, endr, | |
72 | 4, 5, 6, endr, | |
73 | 7, 8, 9, endr; | |
74 | // [ 1 2 3 ] | |
75 | // [ 4 5 6 ] | |
76 | // [ 7 8 9 ] | |
77 | // If by default you need to disable wraping define | |
78 | // BOOST_UBLAS_DEFAULT_NO_WRAP_POLICY, in the compilation options, | |
79 | // so that you avoid typing the "traverse_policy::by_row_no_wrap()". | |
80 | ||
81 | // Plus and minus assign: | |
82 | A <<= fill_policy::index_plus_assign(), | |
83 | 3,2,1; | |
84 | std::cout << A << std::endl; | |
85 | // [ 4 4 4 ] | |
86 | // [ 4 5 6 ] | |
87 | // [ 7 8 9 ] | |
88 | ||
89 | // Matrix from proxy | |
90 | A <<= 0, 1, 2, | |
91 | project(b, range(3,6)), | |
92 | a; | |
93 | std::cout << A << std::endl; | |
94 | // [ 0 1 2 ] | |
95 | // [10 0 1 ] | |
96 | // [ 6 7 8 ] | |
97 | ||
98 | // Matrix from matrix | |
99 | matrix<double> B(6,6); | |
100 | B <<= A, A, | |
101 | A, A; | |
102 | std::cout << B << std::endl; | |
103 | // [ A A ] | |
104 | // [ A A ] | |
105 | ||
106 | // Matrix range (vector is similar) | |
107 | B = zero_matrix<double>(6,6); | |
108 | matrix_range<matrix<double> > mrB (B, range (1, 4), range (1, 4)); | |
109 | mrB <<= 1,2,3,4,5,6,7,8,9; | |
110 | std::cout << B << std::endl; | |
111 | // [ 0 0 0 0 0 0] | |
112 | // [ 0 1 2 3 0 0] | |
113 | // [ 0 4 5 6 0 0] | |
114 | // [ 0 0 0 0 0 0] | |
115 | // [ 0 0 0 0 0 0] | |
116 | // [ 0 0 0 0 0 0] | |
117 | ||
118 | // Horizontal concatenation can be achieved using this trick: | |
119 | matrix<double> BH(3,9); | |
120 | BH <<= A, A, A; | |
121 | std::cout << BH << std::endl; | |
122 | // [ A A A] | |
123 | ||
124 | // Vertical concatenation can be achieved using this trick: | |
125 | matrix<double> BV(9,3); | |
126 | BV <<= A, | |
127 | A, | |
128 | A; | |
129 | std::cout << BV << std::endl; | |
130 | // [ A ] | |
131 | // [ A ] | |
132 | // [ A ] | |
133 | ||
134 | // Watch the difference when assigning matrices for different traverse policies: | |
135 | matrix<double> BR(9,9, 0); | |
136 | BR <<= traverse_policy::by_row(), // This is the default, so this might as well be omitted. | |
137 | A, A, A; | |
138 | std::cout << BR << std::endl; | |
139 | // [ A A A] | |
140 | // [ 0 0 0] | |
141 | // [ 0 0 0] | |
142 | ||
143 | matrix<double> BC(9,9, 0); | |
144 | BC <<= traverse_policy::by_column(), | |
145 | A, A, A; | |
146 | std::cout << BC << std::endl; | |
147 | // [ A 0 0] | |
148 | // [ A 0 0] | |
149 | // [ A 0 0] | |
150 | ||
151 | // The following will throw a run-time exception in debug mode (matrix mid-assignment wrap is not allowed) : | |
152 | // matrix<double> C(7,7); | |
153 | // C <<= A, A, A; | |
154 | ||
155 | // Matrix from matrix with index manipulators | |
156 | matrix<double> C(6,6,0); | |
157 | C <<= A, move(3,0), A; | |
158 | // [ A 0 ] | |
159 | // [ 0 A ] | |
160 | ||
161 | // A faster way for to construct this dense matrix. | |
162 | matrix<double> D(6,6); | |
163 | D <<= A, zero_matrix<double>(3,3), | |
164 | zero_matrix<double>(3,3), A; | |
165 | // [ A 0 ] | |
166 | // [ 0 A ] | |
167 | ||
168 | // The next_row and next_column index manipulators: | |
169 | // note: next_row and next_column functions return | |
170 | // a next_row_manip and and next_column_manip object. | |
171 | // This is the manipulator we used earlier when we disabled | |
172 | // wrapping. | |
173 | matrix<double> E(2,4,0); | |
174 | E <<= 1, 2, next_row(), | |
175 | 3, 4, next_column(),5; | |
176 | std::cout << E << std::endl; | |
177 | // [ 1 2 0 5 ] | |
178 | // [ 3 4 0 0 ] | |
179 | ||
180 | // The begin1 (moves to the begining of the column) index manipulator, begin2 does the same for the row: | |
181 | matrix<double> F(2,4,0); | |
182 | F <<= 1, 2, next_row(), | |
183 | 3, 4, begin1(),5; | |
184 | std::cout << F << std::endl; | |
185 | // [ 1 2 5 0 ] | |
186 | // [ 3 4 0 0 ] | |
187 | ||
188 | // The move (relative) and move_to(absolute) index manipulators (probably the most useful manipulators): | |
189 | matrix<double> G(2,4,0); | |
190 | G <<= 1, 2, move(0,1), 3, | |
191 | move_to(1,3), 4; | |
192 | std::cout << G << std::endl; | |
193 | // [ 1 2 0 3 ] | |
194 | // [ 0 0 0 4 ] | |
195 | ||
196 | // Static equivallents (faster) when sizes are known at compile time: | |
197 | matrix<double> Gs(2,4,0); | |
198 | Gs <<= 1, 2, move<0,1>(), 3, | |
199 | move_to<1,3>(), 4; | |
200 | std::cout << Gs << std::endl; | |
201 | // [ 1 2 0 3 ] | |
202 | // [ 0 0 0 4 ] | |
203 | ||
204 | // Choice of traverse policy (default is "row by row" traverse): | |
205 | ||
206 | matrix<double> H(2,4,0); | |
207 | H <<= 1, 2, 3, 4, | |
208 | 5, 6, 7, 8; | |
209 | std::cout << H << std::endl; | |
210 | // [ 1 2 3 4 ] | |
211 | // [ 5 6 7 8 ] | |
212 | ||
213 | H <<= traverse_policy::by_column(), | |
214 | 1, 2, 3, 4, | |
215 | 5, 6, 7, 8; | |
216 | std::cout << H << std::endl; | |
217 | // [ 1 3 5 7 ] | |
218 | // [ 2 4 6 8 ] | |
219 | ||
220 | // traverse policy can be changed mid assignment if desired. | |
221 | matrix<double> H1(4,4,0); | |
222 | H1 <<= 1, 2, 3, traverse_policy::by_column(), 1, 2, 3; | |
223 | ||
224 | std::cout << H << std::endl; | |
225 | // [1 2 3 1] | |
226 | // [0 0 0 2] | |
227 | // [0 0 0 3] | |
228 | // [0 0 0 0] | |
229 | ||
230 | // note: fill_policy and traverse_policy are namespaces, so you can use them | |
231 | // by a using statement. | |
232 | ||
233 | // For compressed and coordinate matrix types a push_back or insert fill policy can be chosen for faster assginment: | |
234 | compressed_matrix<double> I(2, 2); | |
235 | I <<= fill_policy::sparse_push_back(), | |
236 | 0, 1, 2, 3; | |
237 | std::cout << I << std::endl; | |
238 | // [ 0 1 ] | |
239 | // [ 2 3 ] | |
240 | ||
241 | coordinate_matrix<double> J(2,2); | |
242 | J<<=fill_policy::sparse_insert(), | |
243 | 1, 2, 3, 4; | |
244 | std::cout << J << std::endl; | |
245 | // [ 1 2 ] | |
246 | // [ 3 4 ] | |
247 | ||
248 | // A sparse matrix from another matrix works as with other types. | |
249 | coordinate_matrix<double> K(3,3); | |
250 | K<<=fill_policy::sparse_insert(), | |
251 | J; | |
252 | std::cout << K << std::endl; | |
253 | // [ 1 2 0 ] | |
254 | // [ 3 4 0 ] | |
255 | // [ 0 0 0 ] | |
256 | ||
257 | // Be careful this will not work: | |
258 | //compressed_matrix<double> J2(4,4); | |
259 | //J2<<=fill_policy::sparse_push_back(), | |
260 | // J,J; | |
261 | // That's because the second J2's elements | |
262 | // are attempted to be assigned at positions | |
263 | // that come before the elements already pushed. | |
264 | // Unfortunatelly that's the only thing you can do in this case | |
265 | // (or of course make a custom agorithm): | |
266 | compressed_matrix<double> J2(4,4); | |
267 | J2<<=fill_policy::sparse_push_back(), | |
268 | J, fill_policy::sparse_insert(), | |
269 | J; | |
270 | ||
271 | std::cout << J2 << std::endl; | |
272 | // [ J J ] | |
273 | // [ 0 0 0 0 ] | |
274 | // [ 0 0 0 0 ] | |
275 | ||
276 | // A different traverse policy doesn't change the result, only they order it is been assigned. | |
277 | coordinate_matrix<double> L(3,3); | |
278 | L<<=fill_policy::sparse_insert(), traverse_policy::by_column(), | |
279 | J; | |
280 | std::cout << L << std::endl; | |
281 | // (same as previous) | |
282 | // [ 1 2 0 ] | |
283 | // [ 3 4 0 ] | |
284 | // [ 0 0 0 ] | |
285 | ||
286 | typedef coordinate_matrix<double>::size_type cmst; | |
287 | const cmst size = 30; | |
288 | //typedef fill_policy::sparse_push_back spb; | |
289 | // Although the above could have been used the following is may be faster if | |
290 | // you use the policy often and for relatively small containers. | |
291 | static fill_policy::sparse_push_back spb; | |
292 | ||
293 | // A block diagonal sparse using a loop: | |
294 | compressed_matrix<double> M(size, size, 4*15); | |
295 | for (cmst i=0; i!=size; i+=J.size1()) | |
296 | M <<= spb, move_to(i,i), J; | |
297 | ||
298 | ||
299 | // If typedef was used above the last expression should start | |
300 | // with M <<= spb()... | |
301 | ||
302 | // Displaying so that blocks can be easily seen: | |
303 | for (unsigned int i=0; i!=M.size1(); i++) { | |
304 | std::cout << M(i,0); | |
305 | for (unsigned int j=1; j!=M.size2(); j++) std::cout << ", " << M(i,j); | |
306 | std::cout << "\n"; | |
307 | } | |
308 | // [ J 0 0 0 ... 0] | |
309 | // [ 0 J 0 0 ... 0] | |
310 | // [ 0 . . . ... 0] | |
311 | // [ 0 0 ... 0 0 J] | |
312 | ||
313 | ||
314 | // A "repeat" trasverser may by provided so that this becomes faster and an on-liner like: | |
315 | // M <<= spb, repeat(0, size, J.size1(), 0, size, J.size1()), J; | |
316 | // An alternate would be to create a :repeater" matrix and vector expression that can be used in other places as well. The latter is probably better, | |
317 | return 0; | |
318 | } | |
319 |