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
main()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
320