Line data Source code
1 : // Copyright (c) 2018-2025 Made to Order Software Corp. All Rights Reserved
2 : //
3 : // https://snapwebsites.org/
4 : // contact@m2osw.com
5 : //
6 : // This program is free software: you can redistribute it and/or modify
7 : // it under the terms of the GNU General Public License as published by
8 : // the Free Software Foundation, either version 3 of the License, or
9 : // (at your option) any later version.
10 : //
11 : // This program is distributed in the hope that it will be useful,
12 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : // GNU General Public License for more details.
15 : //
16 : // You should have received a copy of the GNU General Public License
17 : // along with this program. If not, see <https://www.gnu.org/licenses/>.
18 :
19 :
20 : /** \file
21 : * \brief Verify the various matrix implementation.
22 : *
23 : * This file implements tests to verify that the matrix templates
24 : * work as expected.
25 : */
26 :
27 : // self
28 : //
29 : #include "catch_main.h"
30 :
31 :
32 :
33 : // ignore the == and != against float warnings
34 : //
35 : #pragma GCC diagnostic ignored "-Wfloat-equal"
36 :
37 :
38 : // snapdev
39 : //
40 : #include <snapdev/matrix.h>
41 :
42 :
43 : // last include
44 : //
45 : #include <snapdev/poison.h>
46 :
47 :
48 :
49 :
50 :
51 : namespace
52 : {
53 700 : auto frand = []()
54 : {
55 700 : return static_cast<double>(rand()) / static_cast<double>(rand());
56 : };
57 : }
58 : // no name namespace
59 :
60 :
61 :
62 4 : CATCH_TEST_CASE("matrix_init", "[matrix][math]")
63 : {
64 : // constructor/copy
65 : // and
66 : // zero/identity
67 : //
68 4 : CATCH_GIVEN("constructor")
69 : {
70 4 : CATCH_START_SECTION("matrix: empty")
71 : {
72 1 : snapdev::matrix<double> empty;
73 :
74 1 : CATCH_REQUIRE(empty.rows() == 0);
75 1 : CATCH_REQUIRE(empty.columns() == 0);
76 :
77 1 : snapdev::matrix<double> copy(empty);
78 :
79 1 : CATCH_REQUIRE(empty.rows() == 0);
80 1 : CATCH_REQUIRE(empty.columns() == 0);
81 1 : }
82 4 : CATCH_END_SECTION()
83 :
84 4 : CATCH_START_SECTION("matrix: 2x2")
85 : {
86 1 : snapdev::matrix<double> m(2, 2);
87 :
88 1 : CATCH_REQUIRE(m.rows() == 2);
89 1 : CATCH_REQUIRE(m.columns() == 2);
90 :
91 : // by default we get an identity
92 : //
93 1 : CATCH_REQUIRE(m[0][0] == 1.0);
94 1 : CATCH_REQUIRE(m[0][1] == 0.0);
95 1 : CATCH_REQUIRE(m[1][0] == 0.0);
96 1 : CATCH_REQUIRE(m[1][1] == 1.0);
97 :
98 1 : double const r00(frand());
99 1 : double const r01(frand());
100 1 : double const r10(frand());
101 1 : double const r11(frand());
102 :
103 1 : m[0][0] = r00;
104 1 : m[0][1] = r01;
105 1 : m[1][0] = r10;
106 1 : m[1][1] = r11;
107 :
108 1 : CATCH_REQUIRE(m[0][0] == r00);
109 1 : CATCH_REQUIRE(m[0][1] == r01);
110 1 : CATCH_REQUIRE(m[1][0] == r10);
111 1 : CATCH_REQUIRE(m[1][1] == r11);
112 :
113 1 : snapdev::matrix<double> copy(m);
114 1 : snapdev::matrix<double> c2;
115 :
116 1 : CATCH_REQUIRE(copy[0][0] == r00);
117 1 : CATCH_REQUIRE(copy[0][1] == r01);
118 1 : CATCH_REQUIRE(copy[1][0] == r10);
119 1 : CATCH_REQUIRE(copy[1][1] == r11);
120 :
121 1 : m.clear();
122 :
123 1 : CATCH_REQUIRE(m[0][0] == 0.0);
124 1 : CATCH_REQUIRE(m[0][1] == 0.0);
125 1 : CATCH_REQUIRE(m[1][0] == 0.0);
126 1 : CATCH_REQUIRE(m[1][1] == 0.0);
127 :
128 : // copy is still intact
129 : //
130 1 : CATCH_REQUIRE(copy[0][0] == r00);
131 1 : CATCH_REQUIRE(copy[0][1] == r01);
132 1 : CATCH_REQUIRE(copy[1][0] == r10);
133 1 : CATCH_REQUIRE(copy[1][1] == r11);
134 :
135 1 : CATCH_REQUIRE(c2.rows() == 0);
136 1 : CATCH_REQUIRE(c2.columns() == 0);
137 :
138 1 : c2 = copy;
139 :
140 1 : CATCH_REQUIRE(c2[0][0] == r00);
141 1 : CATCH_REQUIRE(c2[0][1] == r01);
142 1 : CATCH_REQUIRE(c2[1][0] == r10);
143 1 : CATCH_REQUIRE(c2[1][1] == r11);
144 :
145 1 : c2.swap(m);
146 :
147 1 : CATCH_REQUIRE(m[0][0] == r00);
148 1 : CATCH_REQUIRE(m[0][1] == r01);
149 1 : CATCH_REQUIRE(m[1][0] == r10);
150 1 : CATCH_REQUIRE(m[1][1] == r11);
151 :
152 1 : CATCH_REQUIRE(c2[0][0] == 0.0);
153 1 : CATCH_REQUIRE(c2[0][1] == 0.0);
154 1 : CATCH_REQUIRE(c2[1][0] == 0.0);
155 1 : CATCH_REQUIRE(c2[1][1] == 0.0);
156 1 : }
157 4 : CATCH_END_SECTION()
158 :
159 4 : CATCH_START_SECTION("matrix: 3x3")
160 : {
161 1 : snapdev::matrix<double> m(3, 3);
162 :
163 1 : CATCH_REQUIRE(m.rows() == 3);
164 1 : CATCH_REQUIRE(m.columns() == 3);
165 :
166 : // by default we get an identity
167 : //
168 1 : CATCH_REQUIRE(m[0][0] == 1.0);
169 1 : CATCH_REQUIRE(m[0][1] == 0.0);
170 1 : CATCH_REQUIRE(m[0][2] == 0.0);
171 1 : CATCH_REQUIRE(m[1][0] == 0.0);
172 1 : CATCH_REQUIRE(m[1][1] == 1.0);
173 1 : CATCH_REQUIRE(m[1][2] == 0.0);
174 1 : CATCH_REQUIRE(m[2][0] == 0.0);
175 1 : CATCH_REQUIRE(m[2][1] == 0.0);
176 1 : CATCH_REQUIRE(m[2][2] == 1.0);
177 :
178 1 : double const r00(frand());
179 1 : double const r01(frand());
180 1 : double const r02(frand());
181 1 : double const r10(frand());
182 1 : double const r11(frand());
183 1 : double const r12(frand());
184 1 : double const r20(frand());
185 1 : double const r21(frand());
186 1 : double const r22(frand());
187 :
188 1 : m[0][0] = r00;
189 1 : m[0][1] = r01;
190 1 : m[0][2] = r02;
191 1 : m[1][0] = r10;
192 1 : m[1][1] = r11;
193 1 : m[1][2] = r12;
194 1 : m[2][0] = r20;
195 1 : m[2][1] = r21;
196 1 : m[2][2] = r22;
197 :
198 1 : CATCH_REQUIRE(m[0][0] == r00);
199 1 : CATCH_REQUIRE(m[0][1] == r01);
200 1 : CATCH_REQUIRE(m[0][2] == r02);
201 1 : CATCH_REQUIRE(m[1][0] == r10);
202 1 : CATCH_REQUIRE(m[1][1] == r11);
203 1 : CATCH_REQUIRE(m[1][2] == r12);
204 1 : CATCH_REQUIRE(m[2][0] == r20);
205 1 : CATCH_REQUIRE(m[2][1] == r21);
206 1 : CATCH_REQUIRE(m[2][2] == r22);
207 :
208 1 : snapdev::matrix<double> copy(m);
209 1 : snapdev::matrix<double> c2;
210 :
211 1 : CATCH_REQUIRE(copy[0][0] == r00);
212 1 : CATCH_REQUIRE(copy[0][1] == r01);
213 1 : CATCH_REQUIRE(copy[0][2] == r02);
214 1 : CATCH_REQUIRE(copy[1][0] == r10);
215 1 : CATCH_REQUIRE(copy[1][1] == r11);
216 1 : CATCH_REQUIRE(copy[1][2] == r12);
217 1 : CATCH_REQUIRE(copy[2][0] == r20);
218 1 : CATCH_REQUIRE(copy[2][1] == r21);
219 1 : CATCH_REQUIRE(copy[2][2] == r22);
220 :
221 1 : m.clear();
222 :
223 1 : CATCH_REQUIRE(m[0][0] == 0.0);
224 1 : CATCH_REQUIRE(m[0][1] == 0.0);
225 1 : CATCH_REQUIRE(m[0][2] == 0.0);
226 1 : CATCH_REQUIRE(m[1][0] == 0.0);
227 1 : CATCH_REQUIRE(m[1][1] == 0.0);
228 1 : CATCH_REQUIRE(m[1][2] == 0.0);
229 1 : CATCH_REQUIRE(m[2][0] == 0.0);
230 1 : CATCH_REQUIRE(m[2][1] == 0.0);
231 1 : CATCH_REQUIRE(m[2][2] == 0.0);
232 :
233 : // copy is still intact
234 : //
235 1 : CATCH_REQUIRE(copy[0][0] == r00);
236 1 : CATCH_REQUIRE(copy[0][1] == r01);
237 1 : CATCH_REQUIRE(copy[0][2] == r02);
238 1 : CATCH_REQUIRE(copy[1][0] == r10);
239 1 : CATCH_REQUIRE(copy[1][1] == r11);
240 1 : CATCH_REQUIRE(copy[1][2] == r12);
241 1 : CATCH_REQUIRE(copy[2][0] == r20);
242 1 : CATCH_REQUIRE(copy[2][1] == r21);
243 1 : CATCH_REQUIRE(copy[2][2] == r22);
244 :
245 1 : CATCH_REQUIRE(c2.rows() == 0);
246 1 : CATCH_REQUIRE(c2.columns() == 0);
247 :
248 1 : c2 = copy;
249 :
250 1 : CATCH_REQUIRE(c2[0][0] == r00);
251 1 : CATCH_REQUIRE(c2[0][1] == r01);
252 1 : CATCH_REQUIRE(c2[0][2] == r02);
253 1 : CATCH_REQUIRE(c2[1][0] == r10);
254 1 : CATCH_REQUIRE(c2[1][1] == r11);
255 1 : CATCH_REQUIRE(c2[1][2] == r12);
256 1 : CATCH_REQUIRE(c2[2][0] == r20);
257 1 : CATCH_REQUIRE(c2[2][1] == r21);
258 1 : CATCH_REQUIRE(c2[2][2] == r22);
259 :
260 1 : c2.swap(m);
261 :
262 1 : CATCH_REQUIRE(m[0][0] == r00);
263 1 : CATCH_REQUIRE(m[0][1] == r01);
264 1 : CATCH_REQUIRE(m[0][2] == r02);
265 1 : CATCH_REQUIRE(m[1][0] == r10);
266 1 : CATCH_REQUIRE(m[1][1] == r11);
267 1 : CATCH_REQUIRE(m[1][2] == r12);
268 1 : CATCH_REQUIRE(m[2][0] == r20);
269 1 : CATCH_REQUIRE(m[2][1] == r21);
270 1 : CATCH_REQUIRE(m[2][2] == r22);
271 :
272 1 : CATCH_REQUIRE(c2[0][0] == 0.0);
273 1 : CATCH_REQUIRE(c2[0][1] == 0.0);
274 1 : CATCH_REQUIRE(c2[0][2] == 0.0);
275 1 : CATCH_REQUIRE(c2[1][0] == 0.0);
276 1 : CATCH_REQUIRE(c2[1][1] == 0.0);
277 1 : CATCH_REQUIRE(c2[1][2] == 0.0);
278 1 : CATCH_REQUIRE(c2[2][0] == 0.0);
279 1 : CATCH_REQUIRE(c2[2][1] == 0.0);
280 1 : CATCH_REQUIRE(c2[2][2] == 0.0);
281 1 : }
282 4 : CATCH_END_SECTION()
283 :
284 4 : CATCH_START_SECTION("matrix: 4x4")
285 : {
286 1 : snapdev::matrix<double> m(4, 4);
287 :
288 1 : CATCH_REQUIRE(m.rows() == 4);
289 1 : CATCH_REQUIRE(m.columns() == 4);
290 :
291 : // by default we get an identity
292 : //
293 1 : CATCH_REQUIRE(m[0][0] == 1.0);
294 1 : CATCH_REQUIRE(m[0][1] == 0.0);
295 1 : CATCH_REQUIRE(m[0][2] == 0.0);
296 1 : CATCH_REQUIRE(m[0][3] == 0.0);
297 1 : CATCH_REQUIRE(m[1][0] == 0.0);
298 1 : CATCH_REQUIRE(m[1][1] == 1.0);
299 1 : CATCH_REQUIRE(m[1][2] == 0.0);
300 1 : CATCH_REQUIRE(m[1][3] == 0.0);
301 1 : CATCH_REQUIRE(m[2][0] == 0.0);
302 1 : CATCH_REQUIRE(m[2][1] == 0.0);
303 1 : CATCH_REQUIRE(m[2][2] == 1.0);
304 1 : CATCH_REQUIRE(m[2][3] == 0.0);
305 1 : CATCH_REQUIRE(m[3][0] == 0.0);
306 1 : CATCH_REQUIRE(m[3][1] == 0.0);
307 1 : CATCH_REQUIRE(m[3][2] == 0.0);
308 1 : CATCH_REQUIRE(m[3][3] == 1.0);
309 :
310 1 : double const r00(frand());
311 1 : double const r01(frand());
312 1 : double const r02(frand());
313 1 : double const r03(frand());
314 1 : double const r10(frand());
315 1 : double const r11(frand());
316 1 : double const r12(frand());
317 1 : double const r13(frand());
318 1 : double const r20(frand());
319 1 : double const r21(frand());
320 1 : double const r22(frand());
321 1 : double const r23(frand());
322 1 : double const r30(frand());
323 1 : double const r31(frand());
324 1 : double const r32(frand());
325 1 : double const r33(frand());
326 :
327 1 : m[0][0] = r00;
328 1 : m[0][1] = r01;
329 1 : m[0][2] = r02;
330 1 : m[0][3] = r03;
331 1 : m[1][0] = r10;
332 1 : m[1][1] = r11;
333 1 : m[1][2] = r12;
334 1 : m[1][3] = r13;
335 1 : m[2][0] = r20;
336 1 : m[2][1] = r21;
337 1 : m[2][2] = r22;
338 1 : m[2][3] = r23;
339 1 : m[3][0] = r30;
340 1 : m[3][1] = r31;
341 1 : m[3][2] = r32;
342 1 : m[3][3] = r33;
343 :
344 1 : CATCH_REQUIRE(m[0][0] == r00);
345 1 : CATCH_REQUIRE(m[0][1] == r01);
346 1 : CATCH_REQUIRE(m[0][2] == r02);
347 1 : CATCH_REQUIRE(m[0][3] == r03);
348 1 : CATCH_REQUIRE(m[1][0] == r10);
349 1 : CATCH_REQUIRE(m[1][1] == r11);
350 1 : CATCH_REQUIRE(m[1][2] == r12);
351 1 : CATCH_REQUIRE(m[1][3] == r13);
352 1 : CATCH_REQUIRE(m[2][0] == r20);
353 1 : CATCH_REQUIRE(m[2][1] == r21);
354 1 : CATCH_REQUIRE(m[2][2] == r22);
355 1 : CATCH_REQUIRE(m[2][3] == r23);
356 1 : CATCH_REQUIRE(m[3][0] == r30);
357 1 : CATCH_REQUIRE(m[3][1] == r31);
358 1 : CATCH_REQUIRE(m[3][2] == r32);
359 1 : CATCH_REQUIRE(m[3][3] == r33);
360 :
361 1 : snapdev::matrix<double> copy(m);
362 1 : snapdev::matrix<double> c2;
363 :
364 1 : CATCH_REQUIRE(copy[0][0] == r00);
365 1 : CATCH_REQUIRE(copy[0][1] == r01);
366 1 : CATCH_REQUIRE(copy[0][2] == r02);
367 1 : CATCH_REQUIRE(copy[0][3] == r03);
368 1 : CATCH_REQUIRE(copy[1][0] == r10);
369 1 : CATCH_REQUIRE(copy[1][1] == r11);
370 1 : CATCH_REQUIRE(copy[1][2] == r12);
371 1 : CATCH_REQUIRE(copy[1][3] == r13);
372 1 : CATCH_REQUIRE(copy[2][0] == r20);
373 1 : CATCH_REQUIRE(copy[2][1] == r21);
374 1 : CATCH_REQUIRE(copy[2][2] == r22);
375 1 : CATCH_REQUIRE(copy[2][3] == r23);
376 1 : CATCH_REQUIRE(copy[3][0] == r30);
377 1 : CATCH_REQUIRE(copy[3][1] == r31);
378 1 : CATCH_REQUIRE(copy[3][2] == r32);
379 1 : CATCH_REQUIRE(copy[3][3] == r33);
380 :
381 1 : m.clear();
382 :
383 1 : CATCH_REQUIRE(m[0][0] == 0.0);
384 1 : CATCH_REQUIRE(m[0][1] == 0.0);
385 1 : CATCH_REQUIRE(m[0][2] == 0.0);
386 1 : CATCH_REQUIRE(m[0][3] == 0.0);
387 1 : CATCH_REQUIRE(m[1][0] == 0.0);
388 1 : CATCH_REQUIRE(m[1][1] == 0.0);
389 1 : CATCH_REQUIRE(m[1][2] == 0.0);
390 1 : CATCH_REQUIRE(m[1][3] == 0.0);
391 1 : CATCH_REQUIRE(m[2][0] == 0.0);
392 1 : CATCH_REQUIRE(m[2][1] == 0.0);
393 1 : CATCH_REQUIRE(m[2][2] == 0.0);
394 1 : CATCH_REQUIRE(m[2][3] == 0.0);
395 :
396 : // copy is still intact
397 : //
398 1 : CATCH_REQUIRE(copy[0][0] == r00);
399 1 : CATCH_REQUIRE(copy[0][1] == r01);
400 1 : CATCH_REQUIRE(copy[0][2] == r02);
401 1 : CATCH_REQUIRE(copy[0][3] == r03);
402 1 : CATCH_REQUIRE(copy[1][0] == r10);
403 1 : CATCH_REQUIRE(copy[1][1] == r11);
404 1 : CATCH_REQUIRE(copy[1][2] == r12);
405 1 : CATCH_REQUIRE(copy[1][3] == r13);
406 1 : CATCH_REQUIRE(copy[2][0] == r20);
407 1 : CATCH_REQUIRE(copy[2][1] == r21);
408 1 : CATCH_REQUIRE(copy[2][2] == r22);
409 1 : CATCH_REQUIRE(copy[2][3] == r23);
410 1 : CATCH_REQUIRE(copy[3][0] == r30);
411 1 : CATCH_REQUIRE(copy[3][1] == r31);
412 1 : CATCH_REQUIRE(copy[3][2] == r32);
413 1 : CATCH_REQUIRE(copy[3][3] == r33);
414 :
415 1 : CATCH_REQUIRE(c2.rows() == 0);
416 1 : CATCH_REQUIRE(c2.columns() == 0);
417 :
418 1 : c2 = copy;
419 :
420 1 : CATCH_REQUIRE(c2[0][0] == r00);
421 1 : CATCH_REQUIRE(c2[0][1] == r01);
422 1 : CATCH_REQUIRE(c2[0][2] == r02);
423 1 : CATCH_REQUIRE(c2[0][3] == r03);
424 1 : CATCH_REQUIRE(c2[1][0] == r10);
425 1 : CATCH_REQUIRE(c2[1][1] == r11);
426 1 : CATCH_REQUIRE(c2[1][2] == r12);
427 1 : CATCH_REQUIRE(c2[1][3] == r13);
428 1 : CATCH_REQUIRE(c2[2][0] == r20);
429 1 : CATCH_REQUIRE(c2[2][1] == r21);
430 1 : CATCH_REQUIRE(c2[2][2] == r22);
431 1 : CATCH_REQUIRE(c2[2][3] == r23);
432 1 : CATCH_REQUIRE(c2[3][0] == r30);
433 1 : CATCH_REQUIRE(c2[3][1] == r31);
434 1 : CATCH_REQUIRE(c2[3][2] == r32);
435 1 : CATCH_REQUIRE(c2[3][3] == r33);
436 :
437 : //std::cout << c2 << std::endl;
438 :
439 1 : c2.swap(m);
440 :
441 1 : CATCH_REQUIRE(m[0][0] == r00);
442 1 : CATCH_REQUIRE(m[0][1] == r01);
443 1 : CATCH_REQUIRE(m[0][2] == r02);
444 1 : CATCH_REQUIRE(m[0][3] == r03);
445 1 : CATCH_REQUIRE(m[1][0] == r10);
446 1 : CATCH_REQUIRE(m[1][1] == r11);
447 1 : CATCH_REQUIRE(m[1][2] == r12);
448 1 : CATCH_REQUIRE(m[1][3] == r13);
449 1 : CATCH_REQUIRE(m[2][0] == r20);
450 1 : CATCH_REQUIRE(m[2][1] == r21);
451 1 : CATCH_REQUIRE(m[2][2] == r22);
452 1 : CATCH_REQUIRE(m[2][3] == r23);
453 1 : CATCH_REQUIRE(m[3][0] == r30);
454 1 : CATCH_REQUIRE(m[3][1] == r31);
455 1 : CATCH_REQUIRE(m[3][2] == r32);
456 1 : CATCH_REQUIRE(m[3][3] == r33);
457 :
458 1 : CATCH_REQUIRE(c2[0][0] == 0.0);
459 1 : CATCH_REQUIRE(c2[0][1] == 0.0);
460 1 : CATCH_REQUIRE(c2[0][2] == 0.0);
461 1 : CATCH_REQUIRE(c2[0][3] == 0.0);
462 1 : CATCH_REQUIRE(c2[1][0] == 0.0);
463 1 : CATCH_REQUIRE(c2[1][1] == 0.0);
464 1 : CATCH_REQUIRE(c2[1][2] == 0.0);
465 1 : CATCH_REQUIRE(c2[1][3] == 0.0);
466 1 : CATCH_REQUIRE(c2[2][0] == 0.0);
467 1 : CATCH_REQUIRE(c2[2][1] == 0.0);
468 1 : CATCH_REQUIRE(c2[2][2] == 0.0);
469 1 : CATCH_REQUIRE(c2[2][3] == 0.0);
470 1 : }
471 4 : CATCH_END_SECTION()
472 4 : }
473 4 : }
474 :
475 :
476 8 : CATCH_TEST_CASE("matrix_additive", "[matrix][math]")
477 : {
478 : // create two random 4x4 matrices and make sure the add works
479 : //
480 8 : CATCH_GIVEN("add")
481 : {
482 4 : CATCH_START_SECTION("matrix: a+=<scalar>")
483 : {
484 : // setup A
485 : //
486 1 : snapdev::matrix<double> a(4, 4);
487 :
488 1 : CATCH_REQUIRE(a.rows() == 4);
489 1 : CATCH_REQUIRE(a.columns() == 4);
490 :
491 1 : double const a00(frand());
492 1 : double const a01(frand());
493 1 : double const a02(frand());
494 1 : double const a03(frand());
495 1 : double const a10(frand());
496 1 : double const a11(frand());
497 1 : double const a12(frand());
498 1 : double const a13(frand());
499 1 : double const a20(frand());
500 1 : double const a21(frand());
501 1 : double const a22(frand());
502 1 : double const a23(frand());
503 1 : double const a30(frand());
504 1 : double const a31(frand());
505 1 : double const a32(frand());
506 1 : double const a33(frand());
507 :
508 1 : a[0][0] = a00;
509 1 : a[0][1] = a01;
510 1 : a[0][2] = a02;
511 1 : a[0][3] = a03;
512 1 : a[1][0] = a10;
513 1 : a[1][1] = a11;
514 1 : a[1][2] = a12;
515 1 : a[1][3] = a13;
516 1 : a[2][0] = a20;
517 1 : a[2][1] = a21;
518 1 : a[2][2] = a22;
519 1 : a[2][3] = a23;
520 1 : a[3][0] = a30;
521 1 : a[3][1] = a31;
522 1 : a[3][2] = a32;
523 1 : a[3][3] = a33;
524 :
525 1 : CATCH_REQUIRE(a[0][0] == a00);
526 1 : CATCH_REQUIRE(a[0][1] == a01);
527 1 : CATCH_REQUIRE(a[0][2] == a02);
528 1 : CATCH_REQUIRE(a[0][3] == a03);
529 1 : CATCH_REQUIRE(a[1][0] == a10);
530 1 : CATCH_REQUIRE(a[1][1] == a11);
531 1 : CATCH_REQUIRE(a[1][2] == a12);
532 1 : CATCH_REQUIRE(a[1][3] == a13);
533 1 : CATCH_REQUIRE(a[2][0] == a20);
534 1 : CATCH_REQUIRE(a[2][1] == a21);
535 1 : CATCH_REQUIRE(a[2][2] == a22);
536 1 : CATCH_REQUIRE(a[2][3] == a23);
537 1 : CATCH_REQUIRE(a[3][0] == a30);
538 1 : CATCH_REQUIRE(a[3][1] == a31);
539 1 : CATCH_REQUIRE(a[3][2] == a32);
540 1 : CATCH_REQUIRE(a[3][3] == a33);
541 :
542 : // create a scalar for our test
543 : //
544 1 : double scalar(frand());
545 1 : while(scalar == 0.0)
546 : {
547 0 : scalar = frand();
548 : }
549 :
550 : // run operation A += <scalar>
551 : //
552 1 : a += scalar;
553 :
554 : // this can't fail because we ensure scalar != 0
555 : //
556 1 : CATCH_REQUIRE_FALSE(a[0][0] == a00);
557 1 : CATCH_REQUIRE_FALSE(a[0][1] == a01);
558 1 : CATCH_REQUIRE_FALSE(a[0][2] == a02);
559 1 : CATCH_REQUIRE_FALSE(a[0][3] == a03);
560 1 : CATCH_REQUIRE_FALSE(a[1][0] == a10);
561 1 : CATCH_REQUIRE_FALSE(a[1][1] == a11);
562 1 : CATCH_REQUIRE_FALSE(a[1][2] == a12);
563 1 : CATCH_REQUIRE_FALSE(a[1][3] == a13);
564 1 : CATCH_REQUIRE_FALSE(a[2][0] == a20);
565 1 : CATCH_REQUIRE_FALSE(a[2][1] == a21);
566 1 : CATCH_REQUIRE_FALSE(a[2][2] == a22);
567 1 : CATCH_REQUIRE_FALSE(a[2][3] == a23);
568 1 : CATCH_REQUIRE_FALSE(a[3][0] == a30);
569 1 : CATCH_REQUIRE_FALSE(a[3][1] == a31);
570 1 : CATCH_REQUIRE_FALSE(a[3][2] == a32);
571 1 : CATCH_REQUIRE_FALSE(a[3][3] == a33);
572 :
573 1 : CATCH_REQUIRE(a[0][0] == a00 + scalar);
574 1 : CATCH_REQUIRE(a[0][1] == a01 + scalar);
575 1 : CATCH_REQUIRE(a[0][2] == a02 + scalar);
576 1 : CATCH_REQUIRE(a[0][3] == a03 + scalar);
577 1 : CATCH_REQUIRE(a[1][0] == a10 + scalar);
578 1 : CATCH_REQUIRE(a[1][1] == a11 + scalar);
579 1 : CATCH_REQUIRE(a[1][2] == a12 + scalar);
580 1 : CATCH_REQUIRE(a[1][3] == a13 + scalar);
581 1 : CATCH_REQUIRE(a[2][0] == a20 + scalar);
582 1 : CATCH_REQUIRE(a[2][1] == a21 + scalar);
583 1 : CATCH_REQUIRE(a[2][2] == a22 + scalar);
584 1 : CATCH_REQUIRE(a[2][3] == a23 + scalar);
585 1 : CATCH_REQUIRE(a[3][0] == a30 + scalar);
586 1 : CATCH_REQUIRE(a[3][1] == a31 + scalar);
587 1 : CATCH_REQUIRE(a[3][2] == a32 + scalar);
588 1 : CATCH_REQUIRE(a[3][3] == a33 + scalar);
589 1 : }
590 4 : CATCH_END_SECTION()
591 :
592 4 : CATCH_START_SECTION("matrix: b=a+<scalar>")
593 : {
594 : // setup A
595 : //
596 1 : snapdev::matrix<double> a(4, 4);
597 :
598 1 : CATCH_REQUIRE(a.rows() == 4);
599 1 : CATCH_REQUIRE(a.columns() == 4);
600 :
601 1 : double const a00(frand());
602 1 : double const a01(frand());
603 1 : double const a02(frand());
604 1 : double const a03(frand());
605 1 : double const a10(frand());
606 1 : double const a11(frand());
607 1 : double const a12(frand());
608 1 : double const a13(frand());
609 1 : double const a20(frand());
610 1 : double const a21(frand());
611 1 : double const a22(frand());
612 1 : double const a23(frand());
613 1 : double const a30(frand());
614 1 : double const a31(frand());
615 1 : double const a32(frand());
616 1 : double const a33(frand());
617 :
618 1 : a[0][0] = a00;
619 1 : a[0][1] = a01;
620 1 : a[0][2] = a02;
621 1 : a[0][3] = a03;
622 1 : a[1][0] = a10;
623 1 : a[1][1] = a11;
624 1 : a[1][2] = a12;
625 1 : a[1][3] = a13;
626 1 : a[2][0] = a20;
627 1 : a[2][1] = a21;
628 1 : a[2][2] = a22;
629 1 : a[2][3] = a23;
630 1 : a[3][0] = a30;
631 1 : a[3][1] = a31;
632 1 : a[3][2] = a32;
633 1 : a[3][3] = a33;
634 :
635 1 : CATCH_REQUIRE(a[0][0] == a00);
636 1 : CATCH_REQUIRE(a[0][1] == a01);
637 1 : CATCH_REQUIRE(a[0][2] == a02);
638 1 : CATCH_REQUIRE(a[0][3] == a03);
639 1 : CATCH_REQUIRE(a[1][0] == a10);
640 1 : CATCH_REQUIRE(a[1][1] == a11);
641 1 : CATCH_REQUIRE(a[1][2] == a12);
642 1 : CATCH_REQUIRE(a[1][3] == a13);
643 1 : CATCH_REQUIRE(a[2][0] == a20);
644 1 : CATCH_REQUIRE(a[2][1] == a21);
645 1 : CATCH_REQUIRE(a[2][2] == a22);
646 1 : CATCH_REQUIRE(a[2][3] == a23);
647 1 : CATCH_REQUIRE(a[3][0] == a30);
648 1 : CATCH_REQUIRE(a[3][1] == a31);
649 1 : CATCH_REQUIRE(a[3][2] == a32);
650 1 : CATCH_REQUIRE(a[3][3] == a33);
651 :
652 : // setup B
653 : //
654 1 : snapdev::matrix<double> b(4, 4);
655 :
656 1 : CATCH_REQUIRE(b.rows() == 4);
657 1 : CATCH_REQUIRE(b.columns() == 4);
658 :
659 1 : double const b00(frand());
660 1 : double const b01(frand());
661 1 : double const b02(frand());
662 1 : double const b03(frand());
663 1 : double const b10(frand());
664 1 : double const b11(frand());
665 1 : double const b12(frand());
666 1 : double const b13(frand());
667 1 : double const b20(frand());
668 1 : double const b21(frand());
669 1 : double const b22(frand());
670 1 : double const b23(frand());
671 1 : double const b30(frand());
672 1 : double const b31(frand());
673 1 : double const b32(frand());
674 1 : double const b33(frand());
675 :
676 1 : b[0][0] = b00;
677 1 : b[0][1] = b01;
678 1 : b[0][2] = b02;
679 1 : b[0][3] = b03;
680 1 : b[1][0] = b10;
681 1 : b[1][1] = b11;
682 1 : b[1][2] = b12;
683 1 : b[1][3] = b13;
684 1 : b[2][0] = b20;
685 1 : b[2][1] = b21;
686 1 : b[2][2] = b22;
687 1 : b[2][3] = b23;
688 1 : b[3][0] = b30;
689 1 : b[3][1] = b31;
690 1 : b[3][2] = b32;
691 1 : b[3][3] = b33;
692 :
693 1 : CATCH_REQUIRE(b[0][0] == b00);
694 1 : CATCH_REQUIRE(b[0][1] == b01);
695 1 : CATCH_REQUIRE(b[0][2] == b02);
696 1 : CATCH_REQUIRE(b[0][3] == b03);
697 1 : CATCH_REQUIRE(b[1][0] == b10);
698 1 : CATCH_REQUIRE(b[1][1] == b11);
699 1 : CATCH_REQUIRE(b[1][2] == b12);
700 1 : CATCH_REQUIRE(b[1][3] == b13);
701 1 : CATCH_REQUIRE(b[2][0] == b20);
702 1 : CATCH_REQUIRE(b[2][1] == b21);
703 1 : CATCH_REQUIRE(b[2][2] == b22);
704 1 : CATCH_REQUIRE(b[2][3] == b23);
705 1 : CATCH_REQUIRE(b[3][0] == b30);
706 1 : CATCH_REQUIRE(b[3][1] == b31);
707 1 : CATCH_REQUIRE(b[3][2] == b32);
708 1 : CATCH_REQUIRE(b[3][3] == b33);
709 :
710 : // create a scalar for our test
711 : //
712 1 : double scalar(frand());
713 1 : while(scalar == 0.0)
714 : {
715 0 : scalar = frand();
716 : }
717 :
718 : // run operation B = A + <scalar>
719 : //
720 1 : b = a + scalar;
721 :
722 1 : CATCH_REQUIRE(a[0][0] == a00);
723 1 : CATCH_REQUIRE(a[0][1] == a01);
724 1 : CATCH_REQUIRE(a[0][2] == a02);
725 1 : CATCH_REQUIRE(a[0][3] == a03);
726 1 : CATCH_REQUIRE(a[1][0] == a10);
727 1 : CATCH_REQUIRE(a[1][1] == a11);
728 1 : CATCH_REQUIRE(a[1][2] == a12);
729 1 : CATCH_REQUIRE(a[1][3] == a13);
730 1 : CATCH_REQUIRE(a[2][0] == a20);
731 1 : CATCH_REQUIRE(a[2][1] == a21);
732 1 : CATCH_REQUIRE(a[2][2] == a22);
733 1 : CATCH_REQUIRE(a[2][3] == a23);
734 1 : CATCH_REQUIRE(a[3][0] == a30);
735 1 : CATCH_REQUIRE(a[3][1] == a31);
736 1 : CATCH_REQUIRE(a[3][2] == a32);
737 1 : CATCH_REQUIRE(a[3][3] == a33);
738 :
739 : // this can't fail because we ensure scalar != 0
740 : //
741 1 : CATCH_REQUIRE_FALSE(b[0][0] == b00);
742 1 : CATCH_REQUIRE_FALSE(b[0][1] == b01);
743 1 : CATCH_REQUIRE_FALSE(b[0][2] == b02);
744 1 : CATCH_REQUIRE_FALSE(b[0][3] == b03);
745 1 : CATCH_REQUIRE_FALSE(b[1][0] == b10);
746 1 : CATCH_REQUIRE_FALSE(b[1][1] == b11);
747 1 : CATCH_REQUIRE_FALSE(b[1][2] == b12);
748 1 : CATCH_REQUIRE_FALSE(b[1][3] == b13);
749 1 : CATCH_REQUIRE_FALSE(b[2][0] == b20);
750 1 : CATCH_REQUIRE_FALSE(b[2][1] == b21);
751 1 : CATCH_REQUIRE_FALSE(b[2][2] == b22);
752 1 : CATCH_REQUIRE_FALSE(b[2][3] == b23);
753 1 : CATCH_REQUIRE_FALSE(b[3][0] == b30);
754 1 : CATCH_REQUIRE_FALSE(b[3][1] == b31);
755 1 : CATCH_REQUIRE_FALSE(b[3][2] == b32);
756 1 : CATCH_REQUIRE_FALSE(b[3][3] == b33);
757 :
758 1 : CATCH_REQUIRE(b[0][0] == a00 + scalar);
759 1 : CATCH_REQUIRE(b[0][1] == a01 + scalar);
760 1 : CATCH_REQUIRE(b[0][2] == a02 + scalar);
761 1 : CATCH_REQUIRE(b[0][3] == a03 + scalar);
762 1 : CATCH_REQUIRE(b[1][0] == a10 + scalar);
763 1 : CATCH_REQUIRE(b[1][1] == a11 + scalar);
764 1 : CATCH_REQUIRE(b[1][2] == a12 + scalar);
765 1 : CATCH_REQUIRE(b[1][3] == a13 + scalar);
766 1 : CATCH_REQUIRE(b[2][0] == a20 + scalar);
767 1 : CATCH_REQUIRE(b[2][1] == a21 + scalar);
768 1 : CATCH_REQUIRE(b[2][2] == a22 + scalar);
769 1 : CATCH_REQUIRE(b[2][3] == a23 + scalar);
770 1 : CATCH_REQUIRE(b[3][0] == a30 + scalar);
771 1 : CATCH_REQUIRE(b[3][1] == a31 + scalar);
772 1 : CATCH_REQUIRE(b[3][2] == a32 + scalar);
773 1 : CATCH_REQUIRE(b[3][3] == a33 + scalar);
774 1 : }
775 4 : CATCH_END_SECTION()
776 :
777 4 : CATCH_START_SECTION("matrix: c=a+b")
778 : {
779 : // setup A
780 : //
781 1 : snapdev::matrix<double> a(4, 4);
782 :
783 1 : CATCH_REQUIRE(a.rows() == 4);
784 1 : CATCH_REQUIRE(a.columns() == 4);
785 :
786 1 : double const a00(frand());
787 1 : double const a01(frand());
788 1 : double const a02(frand());
789 1 : double const a03(frand());
790 1 : double const a10(frand());
791 1 : double const a11(frand());
792 1 : double const a12(frand());
793 1 : double const a13(frand());
794 1 : double const a20(frand());
795 1 : double const a21(frand());
796 1 : double const a22(frand());
797 1 : double const a23(frand());
798 1 : double const a30(frand());
799 1 : double const a31(frand());
800 1 : double const a32(frand());
801 1 : double const a33(frand());
802 :
803 1 : a[0][0] = a00;
804 1 : a[0][1] = a01;
805 1 : a[0][2] = a02;
806 1 : a[0][3] = a03;
807 1 : a[1][0] = a10;
808 1 : a[1][1] = a11;
809 1 : a[1][2] = a12;
810 1 : a[1][3] = a13;
811 1 : a[2][0] = a20;
812 1 : a[2][1] = a21;
813 1 : a[2][2] = a22;
814 1 : a[2][3] = a23;
815 1 : a[3][0] = a30;
816 1 : a[3][1] = a31;
817 1 : a[3][2] = a32;
818 1 : a[3][3] = a33;
819 :
820 1 : CATCH_REQUIRE(a[0][0] == a00);
821 1 : CATCH_REQUIRE(a[0][1] == a01);
822 1 : CATCH_REQUIRE(a[0][2] == a02);
823 1 : CATCH_REQUIRE(a[0][3] == a03);
824 1 : CATCH_REQUIRE(a[1][0] == a10);
825 1 : CATCH_REQUIRE(a[1][1] == a11);
826 1 : CATCH_REQUIRE(a[1][2] == a12);
827 1 : CATCH_REQUIRE(a[1][3] == a13);
828 1 : CATCH_REQUIRE(a[2][0] == a20);
829 1 : CATCH_REQUIRE(a[2][1] == a21);
830 1 : CATCH_REQUIRE(a[2][2] == a22);
831 1 : CATCH_REQUIRE(a[2][3] == a23);
832 1 : CATCH_REQUIRE(a[3][0] == a30);
833 1 : CATCH_REQUIRE(a[3][1] == a31);
834 1 : CATCH_REQUIRE(a[3][2] == a32);
835 1 : CATCH_REQUIRE(a[3][3] == a33);
836 :
837 : // setup B
838 : //
839 1 : snapdev::matrix<double> b(4, 4);
840 :
841 1 : CATCH_REQUIRE(b.rows() == 4);
842 1 : CATCH_REQUIRE(b.columns() == 4);
843 :
844 1 : double const b00(frand());
845 1 : double const b01(frand());
846 1 : double const b02(frand());
847 1 : double const b03(frand());
848 1 : double const b10(frand());
849 1 : double const b11(frand());
850 1 : double const b12(frand());
851 1 : double const b13(frand());
852 1 : double const b20(frand());
853 1 : double const b21(frand());
854 1 : double const b22(frand());
855 1 : double const b23(frand());
856 1 : double const b30(frand());
857 1 : double const b31(frand());
858 1 : double const b32(frand());
859 1 : double const b33(frand());
860 :
861 1 : b[0][0] = b00;
862 1 : b[0][1] = b01;
863 1 : b[0][2] = b02;
864 1 : b[0][3] = b03;
865 1 : b[1][0] = b10;
866 1 : b[1][1] = b11;
867 1 : b[1][2] = b12;
868 1 : b[1][3] = b13;
869 1 : b[2][0] = b20;
870 1 : b[2][1] = b21;
871 1 : b[2][2] = b22;
872 1 : b[2][3] = b23;
873 1 : b[3][0] = b30;
874 1 : b[3][1] = b31;
875 1 : b[3][2] = b32;
876 1 : b[3][3] = b33;
877 :
878 1 : CATCH_REQUIRE(b[0][0] == b00);
879 1 : CATCH_REQUIRE(b[0][1] == b01);
880 1 : CATCH_REQUIRE(b[0][2] == b02);
881 1 : CATCH_REQUIRE(b[0][3] == b03);
882 1 : CATCH_REQUIRE(b[1][0] == b10);
883 1 : CATCH_REQUIRE(b[1][1] == b11);
884 1 : CATCH_REQUIRE(b[1][2] == b12);
885 1 : CATCH_REQUIRE(b[1][3] == b13);
886 1 : CATCH_REQUIRE(b[2][0] == b20);
887 1 : CATCH_REQUIRE(b[2][1] == b21);
888 1 : CATCH_REQUIRE(b[2][2] == b22);
889 1 : CATCH_REQUIRE(b[2][3] == b23);
890 1 : CATCH_REQUIRE(b[3][0] == b30);
891 1 : CATCH_REQUIRE(b[3][1] == b31);
892 1 : CATCH_REQUIRE(b[3][2] == b32);
893 1 : CATCH_REQUIRE(b[3][3] == b33);
894 :
895 : // setup C
896 : //
897 1 : snapdev::matrix<double> c(4, 4);
898 :
899 1 : CATCH_REQUIRE(c.rows() == 4);
900 1 : CATCH_REQUIRE(c.columns() == 4);
901 :
902 1 : double const c00(frand());
903 1 : double const c01(frand());
904 1 : double const c02(frand());
905 1 : double const c03(frand());
906 1 : double const c10(frand());
907 1 : double const c11(frand());
908 1 : double const c12(frand());
909 1 : double const c13(frand());
910 1 : double const c20(frand());
911 1 : double const c21(frand());
912 1 : double const c22(frand());
913 1 : double const c23(frand());
914 1 : double const c30(frand());
915 1 : double const c31(frand());
916 1 : double const c32(frand());
917 1 : double const c33(frand());
918 :
919 1 : c[0][0] = c00;
920 1 : c[0][1] = c01;
921 1 : c[0][2] = c02;
922 1 : c[0][3] = c03;
923 1 : c[1][0] = c10;
924 1 : c[1][1] = c11;
925 1 : c[1][2] = c12;
926 1 : c[1][3] = c13;
927 1 : c[2][0] = c20;
928 1 : c[2][1] = c21;
929 1 : c[2][2] = c22;
930 1 : c[2][3] = c23;
931 1 : c[3][0] = c30;
932 1 : c[3][1] = c31;
933 1 : c[3][2] = c32;
934 1 : c[3][3] = c33;
935 :
936 1 : CATCH_REQUIRE(c[0][0] == c00);
937 1 : CATCH_REQUIRE(c[0][1] == c01);
938 1 : CATCH_REQUIRE(c[0][2] == c02);
939 1 : CATCH_REQUIRE(c[0][3] == c03);
940 1 : CATCH_REQUIRE(c[1][0] == c10);
941 1 : CATCH_REQUIRE(c[1][1] == c11);
942 1 : CATCH_REQUIRE(c[1][2] == c12);
943 1 : CATCH_REQUIRE(c[1][3] == c13);
944 1 : CATCH_REQUIRE(c[2][0] == c20);
945 1 : CATCH_REQUIRE(c[2][1] == c21);
946 1 : CATCH_REQUIRE(c[2][2] == c22);
947 1 : CATCH_REQUIRE(c[2][3] == c23);
948 1 : CATCH_REQUIRE(c[3][0] == c30);
949 1 : CATCH_REQUIRE(c[3][1] == c31);
950 1 : CATCH_REQUIRE(c[3][2] == c32);
951 1 : CATCH_REQUIRE(c[3][3] == c33);
952 :
953 : // run operation C = A + B
954 : //
955 1 : c = a + b;
956 :
957 1 : CATCH_REQUIRE(a[0][0] == a00);
958 1 : CATCH_REQUIRE(a[0][1] == a01);
959 1 : CATCH_REQUIRE(a[0][2] == a02);
960 1 : CATCH_REQUIRE(a[0][3] == a03);
961 1 : CATCH_REQUIRE(a[1][0] == a10);
962 1 : CATCH_REQUIRE(a[1][1] == a11);
963 1 : CATCH_REQUIRE(a[1][2] == a12);
964 1 : CATCH_REQUIRE(a[1][3] == a13);
965 1 : CATCH_REQUIRE(a[2][0] == a20);
966 1 : CATCH_REQUIRE(a[2][1] == a21);
967 1 : CATCH_REQUIRE(a[2][2] == a22);
968 1 : CATCH_REQUIRE(a[2][3] == a23);
969 1 : CATCH_REQUIRE(a[3][0] == a30);
970 1 : CATCH_REQUIRE(a[3][1] == a31);
971 1 : CATCH_REQUIRE(a[3][2] == a32);
972 1 : CATCH_REQUIRE(a[3][3] == a33);
973 :
974 1 : CATCH_REQUIRE(b[0][0] == b00);
975 1 : CATCH_REQUIRE(b[0][1] == b01);
976 1 : CATCH_REQUIRE(b[0][2] == b02);
977 1 : CATCH_REQUIRE(b[0][3] == b03);
978 1 : CATCH_REQUIRE(b[1][0] == b10);
979 1 : CATCH_REQUIRE(b[1][1] == b11);
980 1 : CATCH_REQUIRE(b[1][2] == b12);
981 1 : CATCH_REQUIRE(b[1][3] == b13);
982 1 : CATCH_REQUIRE(b[2][0] == b20);
983 1 : CATCH_REQUIRE(b[2][1] == b21);
984 1 : CATCH_REQUIRE(b[2][2] == b22);
985 1 : CATCH_REQUIRE(b[2][3] == b23);
986 1 : CATCH_REQUIRE(b[3][0] == b30);
987 1 : CATCH_REQUIRE(b[3][1] == b31);
988 1 : CATCH_REQUIRE(b[3][2] == b32);
989 1 : CATCH_REQUIRE(b[3][3] == b33);
990 :
991 : // this could fail
992 : //CATCH_REQUIRE_FALSE(c[0][0] == c00);
993 : //CATCH_REQUIRE_FALSE(c[0][1] == c01);
994 : //CATCH_REQUIRE_FALSE(c[0][2] == c02);
995 : //CATCH_REQUIRE_FALSE(c[0][3] == c03);
996 : //CATCH_REQUIRE_FALSE(c[1][0] == c10);
997 : //CATCH_REQUIRE_FALSE(c[1][1] == c11);
998 : //CATCH_REQUIRE_FALSE(c[1][2] == c12);
999 : //CATCH_REQUIRE_FALSE(c[1][3] == c13);
1000 : //CATCH_REQUIRE_FALSE(c[2][0] == c20);
1001 : //CATCH_REQUIRE_FALSE(c[2][1] == c21);
1002 : //CATCH_REQUIRE_FALSE(c[2][2] == c22);
1003 : //CATCH_REQUIRE_FALSE(c[2][3] == c23);
1004 : //CATCH_REQUIRE_FALSE(c[3][0] == c30);
1005 : //CATCH_REQUIRE_FALSE(c[3][1] == c31);
1006 : //CATCH_REQUIRE_FALSE(c[3][2] == c32);
1007 : //CATCH_REQUIRE_FALSE(c[3][3] == c33);
1008 :
1009 1 : CATCH_REQUIRE(c[0][0] == a00 + b00);
1010 1 : CATCH_REQUIRE(c[0][1] == a01 + b01);
1011 1 : CATCH_REQUIRE(c[0][2] == a02 + b02);
1012 1 : CATCH_REQUIRE(c[0][3] == a03 + b03);
1013 1 : CATCH_REQUIRE(c[1][0] == a10 + b10);
1014 1 : CATCH_REQUIRE(c[1][1] == a11 + b11);
1015 1 : CATCH_REQUIRE(c[1][2] == a12 + b12);
1016 1 : CATCH_REQUIRE(c[1][3] == a13 + b13);
1017 1 : CATCH_REQUIRE(c[2][0] == a20 + b20);
1018 1 : CATCH_REQUIRE(c[2][1] == a21 + b21);
1019 1 : CATCH_REQUIRE(c[2][2] == a22 + b22);
1020 1 : CATCH_REQUIRE(c[2][3] == a23 + b23);
1021 1 : CATCH_REQUIRE(c[3][0] == a30 + b30);
1022 1 : CATCH_REQUIRE(c[3][1] == a31 + b31);
1023 1 : CATCH_REQUIRE(c[3][2] == a32 + b32);
1024 1 : CATCH_REQUIRE(c[3][3] == a33 + b33);
1025 1 : }
1026 4 : CATCH_END_SECTION()
1027 :
1028 4 : CATCH_START_SECTION("matrix: a+=b")
1029 : {
1030 : // setup A
1031 : //
1032 1 : snapdev::matrix<double> a(4, 4);
1033 :
1034 1 : CATCH_REQUIRE(a.rows() == 4);
1035 1 : CATCH_REQUIRE(a.columns() == 4);
1036 :
1037 1 : double const a00(frand());
1038 1 : double const a01(frand());
1039 1 : double const a02(frand());
1040 1 : double const a03(frand());
1041 1 : double const a10(frand());
1042 1 : double const a11(frand());
1043 1 : double const a12(frand());
1044 1 : double const a13(frand());
1045 1 : double const a20(frand());
1046 1 : double const a21(frand());
1047 1 : double const a22(frand());
1048 1 : double const a23(frand());
1049 1 : double const a30(frand());
1050 1 : double const a31(frand());
1051 1 : double const a32(frand());
1052 1 : double const a33(frand());
1053 :
1054 1 : a[0][0] = a00;
1055 1 : a[0][1] = a01;
1056 1 : a[0][2] = a02;
1057 1 : a[0][3] = a03;
1058 1 : a[1][0] = a10;
1059 1 : a[1][1] = a11;
1060 1 : a[1][2] = a12;
1061 1 : a[1][3] = a13;
1062 1 : a[2][0] = a20;
1063 1 : a[2][1] = a21;
1064 1 : a[2][2] = a22;
1065 1 : a[2][3] = a23;
1066 1 : a[3][0] = a30;
1067 1 : a[3][1] = a31;
1068 1 : a[3][2] = a32;
1069 1 : a[3][3] = a33;
1070 :
1071 1 : CATCH_REQUIRE(a[0][0] == a00);
1072 1 : CATCH_REQUIRE(a[0][1] == a01);
1073 1 : CATCH_REQUIRE(a[0][2] == a02);
1074 1 : CATCH_REQUIRE(a[0][3] == a03);
1075 1 : CATCH_REQUIRE(a[1][0] == a10);
1076 1 : CATCH_REQUIRE(a[1][1] == a11);
1077 1 : CATCH_REQUIRE(a[1][2] == a12);
1078 1 : CATCH_REQUIRE(a[1][3] == a13);
1079 1 : CATCH_REQUIRE(a[2][0] == a20);
1080 1 : CATCH_REQUIRE(a[2][1] == a21);
1081 1 : CATCH_REQUIRE(a[2][2] == a22);
1082 1 : CATCH_REQUIRE(a[2][3] == a23);
1083 1 : CATCH_REQUIRE(a[3][0] == a30);
1084 1 : CATCH_REQUIRE(a[3][1] == a31);
1085 1 : CATCH_REQUIRE(a[3][2] == a32);
1086 1 : CATCH_REQUIRE(a[3][3] == a33);
1087 :
1088 : // setup B
1089 : //
1090 1 : snapdev::matrix<double> b(4, 4);
1091 :
1092 1 : CATCH_REQUIRE(b.rows() == 4);
1093 1 : CATCH_REQUIRE(b.columns() == 4);
1094 :
1095 1 : double const b00(frand());
1096 1 : double const b01(frand());
1097 1 : double const b02(frand());
1098 1 : double const b03(frand());
1099 1 : double const b10(frand());
1100 1 : double const b11(frand());
1101 1 : double const b12(frand());
1102 1 : double const b13(frand());
1103 1 : double const b20(frand());
1104 1 : double const b21(frand());
1105 1 : double const b22(frand());
1106 1 : double const b23(frand());
1107 1 : double const b30(frand());
1108 1 : double const b31(frand());
1109 1 : double const b32(frand());
1110 1 : double const b33(frand());
1111 :
1112 1 : b[0][0] = b00;
1113 1 : b[0][1] = b01;
1114 1 : b[0][2] = b02;
1115 1 : b[0][3] = b03;
1116 1 : b[1][0] = b10;
1117 1 : b[1][1] = b11;
1118 1 : b[1][2] = b12;
1119 1 : b[1][3] = b13;
1120 1 : b[2][0] = b20;
1121 1 : b[2][1] = b21;
1122 1 : b[2][2] = b22;
1123 1 : b[2][3] = b23;
1124 1 : b[3][0] = b30;
1125 1 : b[3][1] = b31;
1126 1 : b[3][2] = b32;
1127 1 : b[3][3] = b33;
1128 :
1129 1 : CATCH_REQUIRE(b[0][0] == b00);
1130 1 : CATCH_REQUIRE(b[0][1] == b01);
1131 1 : CATCH_REQUIRE(b[0][2] == b02);
1132 1 : CATCH_REQUIRE(b[0][3] == b03);
1133 1 : CATCH_REQUIRE(b[1][0] == b10);
1134 1 : CATCH_REQUIRE(b[1][1] == b11);
1135 1 : CATCH_REQUIRE(b[1][2] == b12);
1136 1 : CATCH_REQUIRE(b[1][3] == b13);
1137 1 : CATCH_REQUIRE(b[2][0] == b20);
1138 1 : CATCH_REQUIRE(b[2][1] == b21);
1139 1 : CATCH_REQUIRE(b[2][2] == b22);
1140 1 : CATCH_REQUIRE(b[2][3] == b23);
1141 1 : CATCH_REQUIRE(b[3][0] == b30);
1142 1 : CATCH_REQUIRE(b[3][1] == b31);
1143 1 : CATCH_REQUIRE(b[3][2] == b32);
1144 1 : CATCH_REQUIRE(b[3][3] == b33);
1145 :
1146 : // run operation A += B
1147 : //
1148 1 : a += b;
1149 :
1150 : // this could fail if any bXX is 0.0
1151 : //
1152 : //CATCH_REQUIRE_FALSE(a[0][0] == a00);
1153 : //CATCH_REQUIRE_FALSE(a[0][1] == a01);
1154 : //CATCH_REQUIRE_FALSE(a[0][2] == a02);
1155 : //CATCH_REQUIRE_FALSE(a[0][3] == a03);
1156 : //CATCH_REQUIRE_FALSE(a[1][0] == a10);
1157 : //CATCH_REQUIRE_FALSE(a[1][1] == a11);
1158 : //CATCH_REQUIRE_FALSE(a[1][2] == a12);
1159 : //CATCH_REQUIRE_FALSE(a[1][3] == a13);
1160 : //CATCH_REQUIRE_FALSE(a[2][0] == a20);
1161 : //CATCH_REQUIRE_FALSE(a[2][1] == a21);
1162 : //CATCH_REQUIRE_FALSE(a[2][2] == a22);
1163 : //CATCH_REQUIRE_FALSE(a[2][3] == a23);
1164 : //CATCH_REQUIRE_FALSE(a[3][0] == a30);
1165 : //CATCH_REQUIRE_FALSE(a[3][1] == a31);
1166 : //CATCH_REQUIRE_FALSE(a[3][2] == a32);
1167 : //CATCH_REQUIRE_FALSE(a[3][3] == a33);
1168 :
1169 1 : CATCH_REQUIRE(b[0][0] == b00);
1170 1 : CATCH_REQUIRE(b[0][1] == b01);
1171 1 : CATCH_REQUIRE(b[0][2] == b02);
1172 1 : CATCH_REQUIRE(b[0][3] == b03);
1173 1 : CATCH_REQUIRE(b[1][0] == b10);
1174 1 : CATCH_REQUIRE(b[1][1] == b11);
1175 1 : CATCH_REQUIRE(b[1][2] == b12);
1176 1 : CATCH_REQUIRE(b[1][3] == b13);
1177 1 : CATCH_REQUIRE(b[2][0] == b20);
1178 1 : CATCH_REQUIRE(b[2][1] == b21);
1179 1 : CATCH_REQUIRE(b[2][2] == b22);
1180 1 : CATCH_REQUIRE(b[2][3] == b23);
1181 1 : CATCH_REQUIRE(b[3][0] == b30);
1182 1 : CATCH_REQUIRE(b[3][1] == b31);
1183 1 : CATCH_REQUIRE(b[3][2] == b32);
1184 1 : CATCH_REQUIRE(b[3][3] == b33);
1185 :
1186 1 : CATCH_REQUIRE(a[0][0] == a00 + b00);
1187 1 : CATCH_REQUIRE(a[0][1] == a01 + b01);
1188 1 : CATCH_REQUIRE(a[0][2] == a02 + b02);
1189 1 : CATCH_REQUIRE(a[0][3] == a03 + b03);
1190 1 : CATCH_REQUIRE(a[1][0] == a10 + b10);
1191 1 : CATCH_REQUIRE(a[1][1] == a11 + b11);
1192 1 : CATCH_REQUIRE(a[1][2] == a12 + b12);
1193 1 : CATCH_REQUIRE(a[1][3] == a13 + b13);
1194 1 : CATCH_REQUIRE(a[2][0] == a20 + b20);
1195 1 : CATCH_REQUIRE(a[2][1] == a21 + b21);
1196 1 : CATCH_REQUIRE(a[2][2] == a22 + b22);
1197 1 : CATCH_REQUIRE(a[2][3] == a23 + b23);
1198 1 : CATCH_REQUIRE(a[3][0] == a30 + b30);
1199 1 : CATCH_REQUIRE(a[3][1] == a31 + b31);
1200 1 : CATCH_REQUIRE(a[3][2] == a32 + b32);
1201 1 : CATCH_REQUIRE(a[3][3] == a33 + b33);
1202 1 : }
1203 4 : CATCH_END_SECTION()
1204 8 : }
1205 :
1206 : // create two random 4x4 matrices and make sure the add works
1207 : //
1208 8 : CATCH_GIVEN("sub")
1209 : {
1210 4 : CATCH_START_SECTION("matrix: b=a-<scalar>")
1211 : {
1212 : // setup A
1213 : //
1214 1 : snapdev::matrix<double> a(4, 4);
1215 :
1216 1 : CATCH_REQUIRE(a.rows() == 4);
1217 1 : CATCH_REQUIRE(a.columns() == 4);
1218 :
1219 1 : double const a00(frand());
1220 1 : double const a01(frand());
1221 1 : double const a02(frand());
1222 1 : double const a03(frand());
1223 1 : double const a10(frand());
1224 1 : double const a11(frand());
1225 1 : double const a12(frand());
1226 1 : double const a13(frand());
1227 1 : double const a20(frand());
1228 1 : double const a21(frand());
1229 1 : double const a22(frand());
1230 1 : double const a23(frand());
1231 1 : double const a30(frand());
1232 1 : double const a31(frand());
1233 1 : double const a32(frand());
1234 1 : double const a33(frand());
1235 :
1236 1 : a[0][0] = a00;
1237 1 : a[0][1] = a01;
1238 1 : a[0][2] = a02;
1239 1 : a[0][3] = a03;
1240 1 : a[1][0] = a10;
1241 1 : a[1][1] = a11;
1242 1 : a[1][2] = a12;
1243 1 : a[1][3] = a13;
1244 1 : a[2][0] = a20;
1245 1 : a[2][1] = a21;
1246 1 : a[2][2] = a22;
1247 1 : a[2][3] = a23;
1248 1 : a[3][0] = a30;
1249 1 : a[3][1] = a31;
1250 1 : a[3][2] = a32;
1251 1 : a[3][3] = a33;
1252 :
1253 1 : CATCH_REQUIRE(a[0][0] == a00);
1254 1 : CATCH_REQUIRE(a[0][1] == a01);
1255 1 : CATCH_REQUIRE(a[0][2] == a02);
1256 1 : CATCH_REQUIRE(a[0][3] == a03);
1257 1 : CATCH_REQUIRE(a[1][0] == a10);
1258 1 : CATCH_REQUIRE(a[1][1] == a11);
1259 1 : CATCH_REQUIRE(a[1][2] == a12);
1260 1 : CATCH_REQUIRE(a[1][3] == a13);
1261 1 : CATCH_REQUIRE(a[2][0] == a20);
1262 1 : CATCH_REQUIRE(a[2][1] == a21);
1263 1 : CATCH_REQUIRE(a[2][2] == a22);
1264 1 : CATCH_REQUIRE(a[2][3] == a23);
1265 1 : CATCH_REQUIRE(a[3][0] == a30);
1266 1 : CATCH_REQUIRE(a[3][1] == a31);
1267 1 : CATCH_REQUIRE(a[3][2] == a32);
1268 1 : CATCH_REQUIRE(a[3][3] == a33);
1269 :
1270 : // setup B
1271 : //
1272 1 : snapdev::matrix<double> b(4, 4);
1273 :
1274 1 : CATCH_REQUIRE(b.rows() == 4);
1275 1 : CATCH_REQUIRE(b.columns() == 4);
1276 :
1277 1 : double const b00(frand());
1278 1 : double const b01(frand());
1279 1 : double const b02(frand());
1280 1 : double const b03(frand());
1281 1 : double const b10(frand());
1282 1 : double const b11(frand());
1283 1 : double const b12(frand());
1284 1 : double const b13(frand());
1285 1 : double const b20(frand());
1286 1 : double const b21(frand());
1287 1 : double const b22(frand());
1288 1 : double const b23(frand());
1289 1 : double const b30(frand());
1290 1 : double const b31(frand());
1291 1 : double const b32(frand());
1292 1 : double const b33(frand());
1293 :
1294 1 : b[0][0] = b00;
1295 1 : b[0][1] = b01;
1296 1 : b[0][2] = b02;
1297 1 : b[0][3] = b03;
1298 1 : b[1][0] = b10;
1299 1 : b[1][1] = b11;
1300 1 : b[1][2] = b12;
1301 1 : b[1][3] = b13;
1302 1 : b[2][0] = b20;
1303 1 : b[2][1] = b21;
1304 1 : b[2][2] = b22;
1305 1 : b[2][3] = b23;
1306 1 : b[3][0] = b30;
1307 1 : b[3][1] = b31;
1308 1 : b[3][2] = b32;
1309 1 : b[3][3] = b33;
1310 :
1311 1 : CATCH_REQUIRE(b[0][0] == b00);
1312 1 : CATCH_REQUIRE(b[0][1] == b01);
1313 1 : CATCH_REQUIRE(b[0][2] == b02);
1314 1 : CATCH_REQUIRE(b[0][3] == b03);
1315 1 : CATCH_REQUIRE(b[1][0] == b10);
1316 1 : CATCH_REQUIRE(b[1][1] == b11);
1317 1 : CATCH_REQUIRE(b[1][2] == b12);
1318 1 : CATCH_REQUIRE(b[1][3] == b13);
1319 1 : CATCH_REQUIRE(b[2][0] == b20);
1320 1 : CATCH_REQUIRE(b[2][1] == b21);
1321 1 : CATCH_REQUIRE(b[2][2] == b22);
1322 1 : CATCH_REQUIRE(b[2][3] == b23);
1323 1 : CATCH_REQUIRE(b[3][0] == b30);
1324 1 : CATCH_REQUIRE(b[3][1] == b31);
1325 1 : CATCH_REQUIRE(b[3][2] == b32);
1326 1 : CATCH_REQUIRE(b[3][3] == b33);
1327 :
1328 : // create a scalar for our test
1329 : //
1330 1 : double scalar(frand());
1331 1 : while(scalar == 0.0)
1332 : {
1333 0 : scalar = frand();
1334 : }
1335 :
1336 : // run operation B = A - <scalar>
1337 : //
1338 1 : b = a - scalar;
1339 :
1340 1 : CATCH_REQUIRE(a[0][0] == a00);
1341 1 : CATCH_REQUIRE(a[0][1] == a01);
1342 1 : CATCH_REQUIRE(a[0][2] == a02);
1343 1 : CATCH_REQUIRE(a[0][3] == a03);
1344 1 : CATCH_REQUIRE(a[1][0] == a10);
1345 1 : CATCH_REQUIRE(a[1][1] == a11);
1346 1 : CATCH_REQUIRE(a[1][2] == a12);
1347 1 : CATCH_REQUIRE(a[1][3] == a13);
1348 1 : CATCH_REQUIRE(a[2][0] == a20);
1349 1 : CATCH_REQUIRE(a[2][1] == a21);
1350 1 : CATCH_REQUIRE(a[2][2] == a22);
1351 1 : CATCH_REQUIRE(a[2][3] == a23);
1352 1 : CATCH_REQUIRE(a[3][0] == a30);
1353 1 : CATCH_REQUIRE(a[3][1] == a31);
1354 1 : CATCH_REQUIRE(a[3][2] == a32);
1355 1 : CATCH_REQUIRE(a[3][3] == a33);
1356 :
1357 : // this can't fail because we ensure scalar != 0
1358 : //
1359 1 : CATCH_REQUIRE_FALSE(b[0][0] == b00);
1360 1 : CATCH_REQUIRE_FALSE(b[0][1] == b01);
1361 1 : CATCH_REQUIRE_FALSE(b[0][2] == b02);
1362 1 : CATCH_REQUIRE_FALSE(b[0][3] == b03);
1363 1 : CATCH_REQUIRE_FALSE(b[1][0] == b10);
1364 1 : CATCH_REQUIRE_FALSE(b[1][1] == b11);
1365 1 : CATCH_REQUIRE_FALSE(b[1][2] == b12);
1366 1 : CATCH_REQUIRE_FALSE(b[1][3] == b13);
1367 1 : CATCH_REQUIRE_FALSE(b[2][0] == b20);
1368 1 : CATCH_REQUIRE_FALSE(b[2][1] == b21);
1369 1 : CATCH_REQUIRE_FALSE(b[2][2] == b22);
1370 1 : CATCH_REQUIRE_FALSE(b[2][3] == b23);
1371 1 : CATCH_REQUIRE_FALSE(b[3][0] == b30);
1372 1 : CATCH_REQUIRE_FALSE(b[3][1] == b31);
1373 1 : CATCH_REQUIRE_FALSE(b[3][2] == b32);
1374 1 : CATCH_REQUIRE_FALSE(b[3][3] == b33);
1375 :
1376 1 : CATCH_REQUIRE(b[0][0] == a00 - scalar);
1377 1 : CATCH_REQUIRE(b[0][1] == a01 - scalar);
1378 1 : CATCH_REQUIRE(b[0][2] == a02 - scalar);
1379 1 : CATCH_REQUIRE(b[0][3] == a03 - scalar);
1380 1 : CATCH_REQUIRE(b[1][0] == a10 - scalar);
1381 1 : CATCH_REQUIRE(b[1][1] == a11 - scalar);
1382 1 : CATCH_REQUIRE(b[1][2] == a12 - scalar);
1383 1 : CATCH_REQUIRE(b[1][3] == a13 - scalar);
1384 1 : CATCH_REQUIRE(b[2][0] == a20 - scalar);
1385 1 : CATCH_REQUIRE(b[2][1] == a21 - scalar);
1386 1 : CATCH_REQUIRE(b[2][2] == a22 - scalar);
1387 1 : CATCH_REQUIRE(b[2][3] == a23 - scalar);
1388 1 : CATCH_REQUIRE(b[3][0] == a30 - scalar);
1389 1 : CATCH_REQUIRE(b[3][1] == a31 - scalar);
1390 1 : CATCH_REQUIRE(b[3][2] == a32 - scalar);
1391 1 : CATCH_REQUIRE(b[3][3] == a33 - scalar);
1392 1 : }
1393 4 : CATCH_END_SECTION()
1394 :
1395 4 : CATCH_START_SECTION("matrix: a-=<scalar>")
1396 : {
1397 : // setup A
1398 : //
1399 1 : snapdev::matrix<double> a(4, 4);
1400 :
1401 1 : CATCH_REQUIRE(a.rows() == 4);
1402 1 : CATCH_REQUIRE(a.columns() == 4);
1403 :
1404 1 : double const a00(frand());
1405 1 : double const a01(frand());
1406 1 : double const a02(frand());
1407 1 : double const a03(frand());
1408 1 : double const a10(frand());
1409 1 : double const a11(frand());
1410 1 : double const a12(frand());
1411 1 : double const a13(frand());
1412 1 : double const a20(frand());
1413 1 : double const a21(frand());
1414 1 : double const a22(frand());
1415 1 : double const a23(frand());
1416 1 : double const a30(frand());
1417 1 : double const a31(frand());
1418 1 : double const a32(frand());
1419 1 : double const a33(frand());
1420 :
1421 1 : a[0][0] = a00;
1422 1 : a[0][1] = a01;
1423 1 : a[0][2] = a02;
1424 1 : a[0][3] = a03;
1425 1 : a[1][0] = a10;
1426 1 : a[1][1] = a11;
1427 1 : a[1][2] = a12;
1428 1 : a[1][3] = a13;
1429 1 : a[2][0] = a20;
1430 1 : a[2][1] = a21;
1431 1 : a[2][2] = a22;
1432 1 : a[2][3] = a23;
1433 1 : a[3][0] = a30;
1434 1 : a[3][1] = a31;
1435 1 : a[3][2] = a32;
1436 1 : a[3][3] = a33;
1437 :
1438 1 : CATCH_REQUIRE(a[0][0] == a00);
1439 1 : CATCH_REQUIRE(a[0][1] == a01);
1440 1 : CATCH_REQUIRE(a[0][2] == a02);
1441 1 : CATCH_REQUIRE(a[0][3] == a03);
1442 1 : CATCH_REQUIRE(a[1][0] == a10);
1443 1 : CATCH_REQUIRE(a[1][1] == a11);
1444 1 : CATCH_REQUIRE(a[1][2] == a12);
1445 1 : CATCH_REQUIRE(a[1][3] == a13);
1446 1 : CATCH_REQUIRE(a[2][0] == a20);
1447 1 : CATCH_REQUIRE(a[2][1] == a21);
1448 1 : CATCH_REQUIRE(a[2][2] == a22);
1449 1 : CATCH_REQUIRE(a[2][3] == a23);
1450 1 : CATCH_REQUIRE(a[3][0] == a30);
1451 1 : CATCH_REQUIRE(a[3][1] == a31);
1452 1 : CATCH_REQUIRE(a[3][2] == a32);
1453 1 : CATCH_REQUIRE(a[3][3] == a33);
1454 :
1455 : // create a scalar for our test
1456 : //
1457 1 : double scalar(frand());
1458 1 : while(scalar == 0.0)
1459 : {
1460 0 : scalar = frand();
1461 : }
1462 :
1463 : // run operation A -= <scalar>
1464 : //
1465 1 : a -= scalar;
1466 :
1467 : // this can't fail because we ensure scalar != 0
1468 : //
1469 1 : CATCH_REQUIRE_FALSE(a[0][0] == a00);
1470 1 : CATCH_REQUIRE_FALSE(a[0][1] == a01);
1471 1 : CATCH_REQUIRE_FALSE(a[0][2] == a02);
1472 1 : CATCH_REQUIRE_FALSE(a[0][3] == a03);
1473 1 : CATCH_REQUIRE_FALSE(a[1][0] == a10);
1474 1 : CATCH_REQUIRE_FALSE(a[1][1] == a11);
1475 1 : CATCH_REQUIRE_FALSE(a[1][2] == a12);
1476 1 : CATCH_REQUIRE_FALSE(a[1][3] == a13);
1477 1 : CATCH_REQUIRE_FALSE(a[2][0] == a20);
1478 1 : CATCH_REQUIRE_FALSE(a[2][1] == a21);
1479 1 : CATCH_REQUIRE_FALSE(a[2][2] == a22);
1480 1 : CATCH_REQUIRE_FALSE(a[2][3] == a23);
1481 1 : CATCH_REQUIRE_FALSE(a[3][0] == a30);
1482 1 : CATCH_REQUIRE_FALSE(a[3][1] == a31);
1483 1 : CATCH_REQUIRE_FALSE(a[3][2] == a32);
1484 1 : CATCH_REQUIRE_FALSE(a[3][3] == a33);
1485 :
1486 1 : CATCH_REQUIRE(a[0][0] == a00 - scalar);
1487 1 : CATCH_REQUIRE(a[0][1] == a01 - scalar);
1488 1 : CATCH_REQUIRE(a[0][2] == a02 - scalar);
1489 1 : CATCH_REQUIRE(a[0][3] == a03 - scalar);
1490 1 : CATCH_REQUIRE(a[1][0] == a10 - scalar);
1491 1 : CATCH_REQUIRE(a[1][1] == a11 - scalar);
1492 1 : CATCH_REQUIRE(a[1][2] == a12 - scalar);
1493 1 : CATCH_REQUIRE(a[1][3] == a13 - scalar);
1494 1 : CATCH_REQUIRE(a[2][0] == a20 - scalar);
1495 1 : CATCH_REQUIRE(a[2][1] == a21 - scalar);
1496 1 : CATCH_REQUIRE(a[2][2] == a22 - scalar);
1497 1 : CATCH_REQUIRE(a[2][3] == a23 - scalar);
1498 1 : CATCH_REQUIRE(a[3][0] == a30 - scalar);
1499 1 : CATCH_REQUIRE(a[3][1] == a31 - scalar);
1500 1 : CATCH_REQUIRE(a[3][2] == a32 - scalar);
1501 1 : CATCH_REQUIRE(a[3][3] == a33 - scalar);
1502 1 : }
1503 4 : CATCH_END_SECTION()
1504 :
1505 4 : CATCH_START_SECTION("matrix: c=a-b")
1506 : {
1507 : // setup A
1508 : //
1509 1 : snapdev::matrix<double> a(4, 4);
1510 :
1511 1 : CATCH_REQUIRE(a.rows() == 4);
1512 1 : CATCH_REQUIRE(a.columns() == 4);
1513 :
1514 1 : double const a00(frand());
1515 1 : double const a01(frand());
1516 1 : double const a02(frand());
1517 1 : double const a03(frand());
1518 1 : double const a10(frand());
1519 1 : double const a11(frand());
1520 1 : double const a12(frand());
1521 1 : double const a13(frand());
1522 1 : double const a20(frand());
1523 1 : double const a21(frand());
1524 1 : double const a22(frand());
1525 1 : double const a23(frand());
1526 1 : double const a30(frand());
1527 1 : double const a31(frand());
1528 1 : double const a32(frand());
1529 1 : double const a33(frand());
1530 :
1531 1 : a[0][0] = a00;
1532 1 : a[0][1] = a01;
1533 1 : a[0][2] = a02;
1534 1 : a[0][3] = a03;
1535 1 : a[1][0] = a10;
1536 1 : a[1][1] = a11;
1537 1 : a[1][2] = a12;
1538 1 : a[1][3] = a13;
1539 1 : a[2][0] = a20;
1540 1 : a[2][1] = a21;
1541 1 : a[2][2] = a22;
1542 1 : a[2][3] = a23;
1543 1 : a[3][0] = a30;
1544 1 : a[3][1] = a31;
1545 1 : a[3][2] = a32;
1546 1 : a[3][3] = a33;
1547 :
1548 1 : CATCH_REQUIRE(a[0][0] == a00);
1549 1 : CATCH_REQUIRE(a[0][1] == a01);
1550 1 : CATCH_REQUIRE(a[0][2] == a02);
1551 1 : CATCH_REQUIRE(a[0][3] == a03);
1552 1 : CATCH_REQUIRE(a[1][0] == a10);
1553 1 : CATCH_REQUIRE(a[1][1] == a11);
1554 1 : CATCH_REQUIRE(a[1][2] == a12);
1555 1 : CATCH_REQUIRE(a[1][3] == a13);
1556 1 : CATCH_REQUIRE(a[2][0] == a20);
1557 1 : CATCH_REQUIRE(a[2][1] == a21);
1558 1 : CATCH_REQUIRE(a[2][2] == a22);
1559 1 : CATCH_REQUIRE(a[2][3] == a23);
1560 1 : CATCH_REQUIRE(a[3][0] == a30);
1561 1 : CATCH_REQUIRE(a[3][1] == a31);
1562 1 : CATCH_REQUIRE(a[3][2] == a32);
1563 1 : CATCH_REQUIRE(a[3][3] == a33);
1564 :
1565 : // setup B
1566 : //
1567 1 : snapdev::matrix<double> b(4, 4);
1568 :
1569 1 : CATCH_REQUIRE(b.rows() == 4);
1570 1 : CATCH_REQUIRE(b.columns() == 4);
1571 :
1572 1 : double const b00(frand());
1573 1 : double const b01(frand());
1574 1 : double const b02(frand());
1575 1 : double const b03(frand());
1576 1 : double const b10(frand());
1577 1 : double const b11(frand());
1578 1 : double const b12(frand());
1579 1 : double const b13(frand());
1580 1 : double const b20(frand());
1581 1 : double const b21(frand());
1582 1 : double const b22(frand());
1583 1 : double const b23(frand());
1584 1 : double const b30(frand());
1585 1 : double const b31(frand());
1586 1 : double const b32(frand());
1587 1 : double const b33(frand());
1588 :
1589 1 : b[0][0] = b00;
1590 1 : b[0][1] = b01;
1591 1 : b[0][2] = b02;
1592 1 : b[0][3] = b03;
1593 1 : b[1][0] = b10;
1594 1 : b[1][1] = b11;
1595 1 : b[1][2] = b12;
1596 1 : b[1][3] = b13;
1597 1 : b[2][0] = b20;
1598 1 : b[2][1] = b21;
1599 1 : b[2][2] = b22;
1600 1 : b[2][3] = b23;
1601 1 : b[3][0] = b30;
1602 1 : b[3][1] = b31;
1603 1 : b[3][2] = b32;
1604 1 : b[3][3] = b33;
1605 :
1606 1 : CATCH_REQUIRE(b[0][0] == b00);
1607 1 : CATCH_REQUIRE(b[0][1] == b01);
1608 1 : CATCH_REQUIRE(b[0][2] == b02);
1609 1 : CATCH_REQUIRE(b[0][3] == b03);
1610 1 : CATCH_REQUIRE(b[1][0] == b10);
1611 1 : CATCH_REQUIRE(b[1][1] == b11);
1612 1 : CATCH_REQUIRE(b[1][2] == b12);
1613 1 : CATCH_REQUIRE(b[1][3] == b13);
1614 1 : CATCH_REQUIRE(b[2][0] == b20);
1615 1 : CATCH_REQUIRE(b[2][1] == b21);
1616 1 : CATCH_REQUIRE(b[2][2] == b22);
1617 1 : CATCH_REQUIRE(b[2][3] == b23);
1618 1 : CATCH_REQUIRE(b[3][0] == b30);
1619 1 : CATCH_REQUIRE(b[3][1] == b31);
1620 1 : CATCH_REQUIRE(b[3][2] == b32);
1621 1 : CATCH_REQUIRE(b[3][3] == b33);
1622 :
1623 : // setup C
1624 : //
1625 1 : snapdev::matrix<double> c(4, 4);
1626 :
1627 1 : CATCH_REQUIRE(c.rows() == 4);
1628 1 : CATCH_REQUIRE(c.columns() == 4);
1629 :
1630 1 : double const c00(frand());
1631 1 : double const c01(frand());
1632 1 : double const c02(frand());
1633 1 : double const c03(frand());
1634 1 : double const c10(frand());
1635 1 : double const c11(frand());
1636 1 : double const c12(frand());
1637 1 : double const c13(frand());
1638 1 : double const c20(frand());
1639 1 : double const c21(frand());
1640 1 : double const c22(frand());
1641 1 : double const c23(frand());
1642 1 : double const c30(frand());
1643 1 : double const c31(frand());
1644 1 : double const c32(frand());
1645 1 : double const c33(frand());
1646 :
1647 1 : c[0][0] = c00;
1648 1 : c[0][1] = c01;
1649 1 : c[0][2] = c02;
1650 1 : c[0][3] = c03;
1651 1 : c[1][0] = c10;
1652 1 : c[1][1] = c11;
1653 1 : c[1][2] = c12;
1654 1 : c[1][3] = c13;
1655 1 : c[2][0] = c20;
1656 1 : c[2][1] = c21;
1657 1 : c[2][2] = c22;
1658 1 : c[2][3] = c23;
1659 1 : c[3][0] = c30;
1660 1 : c[3][1] = c31;
1661 1 : c[3][2] = c32;
1662 1 : c[3][3] = c33;
1663 :
1664 1 : CATCH_REQUIRE(c[0][0] == c00);
1665 1 : CATCH_REQUIRE(c[0][1] == c01);
1666 1 : CATCH_REQUIRE(c[0][2] == c02);
1667 1 : CATCH_REQUIRE(c[0][3] == c03);
1668 1 : CATCH_REQUIRE(c[1][0] == c10);
1669 1 : CATCH_REQUIRE(c[1][1] == c11);
1670 1 : CATCH_REQUIRE(c[1][2] == c12);
1671 1 : CATCH_REQUIRE(c[1][3] == c13);
1672 1 : CATCH_REQUIRE(c[2][0] == c20);
1673 1 : CATCH_REQUIRE(c[2][1] == c21);
1674 1 : CATCH_REQUIRE(c[2][2] == c22);
1675 1 : CATCH_REQUIRE(c[2][3] == c23);
1676 1 : CATCH_REQUIRE(c[3][0] == c30);
1677 1 : CATCH_REQUIRE(c[3][1] == c31);
1678 1 : CATCH_REQUIRE(c[3][2] == c32);
1679 1 : CATCH_REQUIRE(c[3][3] == c33);
1680 :
1681 : // run operation C = A - B
1682 : //
1683 1 : c = a - b;
1684 :
1685 1 : CATCH_REQUIRE(a[0][0] == a00);
1686 1 : CATCH_REQUIRE(a[0][1] == a01);
1687 1 : CATCH_REQUIRE(a[0][2] == a02);
1688 1 : CATCH_REQUIRE(a[0][3] == a03);
1689 1 : CATCH_REQUIRE(a[1][0] == a10);
1690 1 : CATCH_REQUIRE(a[1][1] == a11);
1691 1 : CATCH_REQUIRE(a[1][2] == a12);
1692 1 : CATCH_REQUIRE(a[1][3] == a13);
1693 1 : CATCH_REQUIRE(a[2][0] == a20);
1694 1 : CATCH_REQUIRE(a[2][1] == a21);
1695 1 : CATCH_REQUIRE(a[2][2] == a22);
1696 1 : CATCH_REQUIRE(a[2][3] == a23);
1697 1 : CATCH_REQUIRE(a[3][0] == a30);
1698 1 : CATCH_REQUIRE(a[3][1] == a31);
1699 1 : CATCH_REQUIRE(a[3][2] == a32);
1700 1 : CATCH_REQUIRE(a[3][3] == a33);
1701 :
1702 1 : CATCH_REQUIRE(b[0][0] == b00);
1703 1 : CATCH_REQUIRE(b[0][1] == b01);
1704 1 : CATCH_REQUIRE(b[0][2] == b02);
1705 1 : CATCH_REQUIRE(b[0][3] == b03);
1706 1 : CATCH_REQUIRE(b[1][0] == b10);
1707 1 : CATCH_REQUIRE(b[1][1] == b11);
1708 1 : CATCH_REQUIRE(b[1][2] == b12);
1709 1 : CATCH_REQUIRE(b[1][3] == b13);
1710 1 : CATCH_REQUIRE(b[2][0] == b20);
1711 1 : CATCH_REQUIRE(b[2][1] == b21);
1712 1 : CATCH_REQUIRE(b[2][2] == b22);
1713 1 : CATCH_REQUIRE(b[2][3] == b23);
1714 1 : CATCH_REQUIRE(b[3][0] == b30);
1715 1 : CATCH_REQUIRE(b[3][1] == b31);
1716 1 : CATCH_REQUIRE(b[3][2] == b32);
1717 1 : CATCH_REQUIRE(b[3][3] == b33);
1718 :
1719 : // this could fail
1720 : //CATCH_REQUIRE_FALSE(c[0][0] == c00);
1721 : //CATCH_REQUIRE_FALSE(c[0][1] == c01);
1722 : //CATCH_REQUIRE_FALSE(c[0][2] == c02);
1723 : //CATCH_REQUIRE_FALSE(c[0][3] == c03);
1724 : //CATCH_REQUIRE_FALSE(c[1][0] == c10);
1725 : //CATCH_REQUIRE_FALSE(c[1][1] == c11);
1726 : //CATCH_REQUIRE_FALSE(c[1][2] == c12);
1727 : //CATCH_REQUIRE_FALSE(c[1][3] == c13);
1728 : //CATCH_REQUIRE_FALSE(c[2][0] == c20);
1729 : //CATCH_REQUIRE_FALSE(c[2][1] == c21);
1730 : //CATCH_REQUIRE_FALSE(c[2][2] == c22);
1731 : //CATCH_REQUIRE_FALSE(c[2][3] == c23);
1732 : //CATCH_REQUIRE_FALSE(c[3][0] == c30);
1733 : //CATCH_REQUIRE_FALSE(c[3][1] == c31);
1734 : //CATCH_REQUIRE_FALSE(c[3][2] == c32);
1735 : //CATCH_REQUIRE_FALSE(c[3][3] == c33);
1736 :
1737 1 : CATCH_REQUIRE(c[0][0] == a00 - b00);
1738 1 : CATCH_REQUIRE(c[0][1] == a01 - b01);
1739 1 : CATCH_REQUIRE(c[0][2] == a02 - b02);
1740 1 : CATCH_REQUIRE(c[0][3] == a03 - b03);
1741 1 : CATCH_REQUIRE(c[1][0] == a10 - b10);
1742 1 : CATCH_REQUIRE(c[1][1] == a11 - b11);
1743 1 : CATCH_REQUIRE(c[1][2] == a12 - b12);
1744 1 : CATCH_REQUIRE(c[1][3] == a13 - b13);
1745 1 : CATCH_REQUIRE(c[2][0] == a20 - b20);
1746 1 : CATCH_REQUIRE(c[2][1] == a21 - b21);
1747 1 : CATCH_REQUIRE(c[2][2] == a22 - b22);
1748 1 : CATCH_REQUIRE(c[2][3] == a23 - b23);
1749 1 : CATCH_REQUIRE(c[3][0] == a30 - b30);
1750 1 : CATCH_REQUIRE(c[3][1] == a31 - b31);
1751 1 : CATCH_REQUIRE(c[3][2] == a32 - b32);
1752 1 : CATCH_REQUIRE(c[3][3] == a33 - b33);
1753 1 : }
1754 4 : CATCH_END_SECTION()
1755 :
1756 4 : CATCH_START_SECTION("matrix: a-=b")
1757 : {
1758 : // setup A
1759 : //
1760 1 : snapdev::matrix<double> a(4, 4);
1761 :
1762 1 : CATCH_REQUIRE(a.rows() == 4);
1763 1 : CATCH_REQUIRE(a.columns() == 4);
1764 :
1765 1 : double const a00(frand());
1766 1 : double const a01(frand());
1767 1 : double const a02(frand());
1768 1 : double const a03(frand());
1769 1 : double const a10(frand());
1770 1 : double const a11(frand());
1771 1 : double const a12(frand());
1772 1 : double const a13(frand());
1773 1 : double const a20(frand());
1774 1 : double const a21(frand());
1775 1 : double const a22(frand());
1776 1 : double const a23(frand());
1777 1 : double const a30(frand());
1778 1 : double const a31(frand());
1779 1 : double const a32(frand());
1780 1 : double const a33(frand());
1781 :
1782 1 : a[0][0] = a00;
1783 1 : a[0][1] = a01;
1784 1 : a[0][2] = a02;
1785 1 : a[0][3] = a03;
1786 1 : a[1][0] = a10;
1787 1 : a[1][1] = a11;
1788 1 : a[1][2] = a12;
1789 1 : a[1][3] = a13;
1790 1 : a[2][0] = a20;
1791 1 : a[2][1] = a21;
1792 1 : a[2][2] = a22;
1793 1 : a[2][3] = a23;
1794 1 : a[3][0] = a30;
1795 1 : a[3][1] = a31;
1796 1 : a[3][2] = a32;
1797 1 : a[3][3] = a33;
1798 :
1799 1 : CATCH_REQUIRE(a[0][0] == a00);
1800 1 : CATCH_REQUIRE(a[0][1] == a01);
1801 1 : CATCH_REQUIRE(a[0][2] == a02);
1802 1 : CATCH_REQUIRE(a[0][3] == a03);
1803 1 : CATCH_REQUIRE(a[1][0] == a10);
1804 1 : CATCH_REQUIRE(a[1][1] == a11);
1805 1 : CATCH_REQUIRE(a[1][2] == a12);
1806 1 : CATCH_REQUIRE(a[1][3] == a13);
1807 1 : CATCH_REQUIRE(a[2][0] == a20);
1808 1 : CATCH_REQUIRE(a[2][1] == a21);
1809 1 : CATCH_REQUIRE(a[2][2] == a22);
1810 1 : CATCH_REQUIRE(a[2][3] == a23);
1811 1 : CATCH_REQUIRE(a[3][0] == a30);
1812 1 : CATCH_REQUIRE(a[3][1] == a31);
1813 1 : CATCH_REQUIRE(a[3][2] == a32);
1814 1 : CATCH_REQUIRE(a[3][3] == a33);
1815 :
1816 : // setup B
1817 : //
1818 1 : snapdev::matrix<double> b(4, 4);
1819 :
1820 1 : CATCH_REQUIRE(b.rows() == 4);
1821 1 : CATCH_REQUIRE(b.columns() == 4);
1822 :
1823 1 : double const b00(frand());
1824 1 : double const b01(frand());
1825 1 : double const b02(frand());
1826 1 : double const b03(frand());
1827 1 : double const b10(frand());
1828 1 : double const b11(frand());
1829 1 : double const b12(frand());
1830 1 : double const b13(frand());
1831 1 : double const b20(frand());
1832 1 : double const b21(frand());
1833 1 : double const b22(frand());
1834 1 : double const b23(frand());
1835 1 : double const b30(frand());
1836 1 : double const b31(frand());
1837 1 : double const b32(frand());
1838 1 : double const b33(frand());
1839 :
1840 1 : b[0][0] = b00;
1841 1 : b[0][1] = b01;
1842 1 : b[0][2] = b02;
1843 1 : b[0][3] = b03;
1844 1 : b[1][0] = b10;
1845 1 : b[1][1] = b11;
1846 1 : b[1][2] = b12;
1847 1 : b[1][3] = b13;
1848 1 : b[2][0] = b20;
1849 1 : b[2][1] = b21;
1850 1 : b[2][2] = b22;
1851 1 : b[2][3] = b23;
1852 1 : b[3][0] = b30;
1853 1 : b[3][1] = b31;
1854 1 : b[3][2] = b32;
1855 1 : b[3][3] = b33;
1856 :
1857 1 : CATCH_REQUIRE(b[0][0] == b00);
1858 1 : CATCH_REQUIRE(b[0][1] == b01);
1859 1 : CATCH_REQUIRE(b[0][2] == b02);
1860 1 : CATCH_REQUIRE(b[0][3] == b03);
1861 1 : CATCH_REQUIRE(b[1][0] == b10);
1862 1 : CATCH_REQUIRE(b[1][1] == b11);
1863 1 : CATCH_REQUIRE(b[1][2] == b12);
1864 1 : CATCH_REQUIRE(b[1][3] == b13);
1865 1 : CATCH_REQUIRE(b[2][0] == b20);
1866 1 : CATCH_REQUIRE(b[2][1] == b21);
1867 1 : CATCH_REQUIRE(b[2][2] == b22);
1868 1 : CATCH_REQUIRE(b[2][3] == b23);
1869 1 : CATCH_REQUIRE(b[3][0] == b30);
1870 1 : CATCH_REQUIRE(b[3][1] == b31);
1871 1 : CATCH_REQUIRE(b[3][2] == b32);
1872 1 : CATCH_REQUIRE(b[3][3] == b33);
1873 :
1874 : // run operation A -= B
1875 : //
1876 1 : a -= b;
1877 :
1878 : // this could fail if one of bXX is 0.0
1879 : //
1880 : //CATCH_REQUIRE(a[0][0] == a00);
1881 : //CATCH_REQUIRE(a[0][1] == a01);
1882 : //CATCH_REQUIRE(a[0][2] == a02);
1883 : //CATCH_REQUIRE(a[0][3] == a03);
1884 : //CATCH_REQUIRE(a[1][0] == a10);
1885 : //CATCH_REQUIRE(a[1][1] == a11);
1886 : //CATCH_REQUIRE(a[1][2] == a12);
1887 : //CATCH_REQUIRE(a[1][3] == a13);
1888 : //CATCH_REQUIRE(a[2][0] == a20);
1889 : //CATCH_REQUIRE(a[2][1] == a21);
1890 : //CATCH_REQUIRE(a[2][2] == a22);
1891 : //CATCH_REQUIRE(a[2][3] == a23);
1892 : //CATCH_REQUIRE(a[3][0] == a30);
1893 : //CATCH_REQUIRE(a[3][1] == a31);
1894 : //CATCH_REQUIRE(a[3][2] == a32);
1895 : //CATCH_REQUIRE(a[3][3] == a33);
1896 :
1897 1 : CATCH_REQUIRE(b[0][0] == b00);
1898 1 : CATCH_REQUIRE(b[0][1] == b01);
1899 1 : CATCH_REQUIRE(b[0][2] == b02);
1900 1 : CATCH_REQUIRE(b[0][3] == b03);
1901 1 : CATCH_REQUIRE(b[1][0] == b10);
1902 1 : CATCH_REQUIRE(b[1][1] == b11);
1903 1 : CATCH_REQUIRE(b[1][2] == b12);
1904 1 : CATCH_REQUIRE(b[1][3] == b13);
1905 1 : CATCH_REQUIRE(b[2][0] == b20);
1906 1 : CATCH_REQUIRE(b[2][1] == b21);
1907 1 : CATCH_REQUIRE(b[2][2] == b22);
1908 1 : CATCH_REQUIRE(b[2][3] == b23);
1909 1 : CATCH_REQUIRE(b[3][0] == b30);
1910 1 : CATCH_REQUIRE(b[3][1] == b31);
1911 1 : CATCH_REQUIRE(b[3][2] == b32);
1912 1 : CATCH_REQUIRE(b[3][3] == b33);
1913 :
1914 1 : CATCH_REQUIRE(a[0][0] == a00 - b00);
1915 1 : CATCH_REQUIRE(a[0][1] == a01 - b01);
1916 1 : CATCH_REQUIRE(a[0][2] == a02 - b02);
1917 1 : CATCH_REQUIRE(a[0][3] == a03 - b03);
1918 1 : CATCH_REQUIRE(a[1][0] == a10 - b10);
1919 1 : CATCH_REQUIRE(a[1][1] == a11 - b11);
1920 1 : CATCH_REQUIRE(a[1][2] == a12 - b12);
1921 1 : CATCH_REQUIRE(a[1][3] == a13 - b13);
1922 1 : CATCH_REQUIRE(a[2][0] == a20 - b20);
1923 1 : CATCH_REQUIRE(a[2][1] == a21 - b21);
1924 1 : CATCH_REQUIRE(a[2][2] == a22 - b22);
1925 1 : CATCH_REQUIRE(a[2][3] == a23 - b23);
1926 1 : CATCH_REQUIRE(a[3][0] == a30 - b30);
1927 1 : CATCH_REQUIRE(a[3][1] == a31 - b31);
1928 1 : CATCH_REQUIRE(a[3][2] == a32 - b32);
1929 1 : CATCH_REQUIRE(a[3][3] == a33 - b33);
1930 1 : }
1931 4 : CATCH_END_SECTION()
1932 8 : }
1933 8 : }
1934 :
1935 :
1936 8 : CATCH_TEST_CASE("matrix_util", "[matrix][math]")
1937 : {
1938 : // various ways to change the data order, minor, enlarge
1939 : //
1940 8 : CATCH_GIVEN("move data")
1941 : {
1942 5 : CATCH_START_SECTION("matrix: minor")
1943 : {
1944 : // setup A
1945 : //
1946 1 : snapdev::matrix<double> a(4, 4);
1947 :
1948 1 : CATCH_REQUIRE(a.rows() == 4);
1949 1 : CATCH_REQUIRE(a.columns() == 4);
1950 :
1951 1 : double const a00(frand());
1952 1 : double const a01(frand());
1953 1 : double const a02(frand());
1954 1 : double const a03(frand());
1955 1 : double const a10(frand());
1956 1 : double const a11(frand());
1957 1 : double const a12(frand());
1958 1 : double const a13(frand());
1959 1 : double const a20(frand());
1960 1 : double const a21(frand());
1961 1 : double const a22(frand());
1962 1 : double const a23(frand());
1963 1 : double const a30(frand());
1964 1 : double const a31(frand());
1965 1 : double const a32(frand());
1966 1 : double const a33(frand());
1967 :
1968 1 : a[0][0] = a00;
1969 1 : a[0][1] = a01;
1970 1 : a[0][2] = a02;
1971 1 : a[0][3] = a03;
1972 1 : a[1][0] = a10;
1973 1 : a[1][1] = a11;
1974 1 : a[1][2] = a12;
1975 1 : a[1][3] = a13;
1976 1 : a[2][0] = a20;
1977 1 : a[2][1] = a21;
1978 1 : a[2][2] = a22;
1979 1 : a[2][3] = a23;
1980 1 : a[3][0] = a30;
1981 1 : a[3][1] = a31;
1982 1 : a[3][2] = a32;
1983 1 : a[3][3] = a33;
1984 :
1985 1 : CATCH_REQUIRE(a[0][0] == a00);
1986 1 : CATCH_REQUIRE(a[0][1] == a01);
1987 1 : CATCH_REQUIRE(a[0][2] == a02);
1988 1 : CATCH_REQUIRE(a[0][3] == a03);
1989 1 : CATCH_REQUIRE(a[1][0] == a10);
1990 1 : CATCH_REQUIRE(a[1][1] == a11);
1991 1 : CATCH_REQUIRE(a[1][2] == a12);
1992 1 : CATCH_REQUIRE(a[1][3] == a13);
1993 1 : CATCH_REQUIRE(a[2][0] == a20);
1994 1 : CATCH_REQUIRE(a[2][1] == a21);
1995 1 : CATCH_REQUIRE(a[2][2] == a22);
1996 1 : CATCH_REQUIRE(a[2][3] == a23);
1997 1 : CATCH_REQUIRE(a[3][0] == a30);
1998 1 : CATCH_REQUIRE(a[3][1] == a31);
1999 1 : CATCH_REQUIRE(a[3][2] == a32);
2000 1 : CATCH_REQUIRE(a[3][3] == a33);
2001 :
2002 : {
2003 1 : snapdev::matrix<double> p = a.minor_matrix(0, 0);
2004 :
2005 1 : CATCH_REQUIRE(p.rows() == 3);
2006 1 : CATCH_REQUIRE(p.columns() == 3);
2007 :
2008 1 : CATCH_REQUIRE(p[0][0] == a11);
2009 1 : CATCH_REQUIRE(p[0][1] == a12);
2010 1 : CATCH_REQUIRE(p[0][2] == a13);
2011 1 : CATCH_REQUIRE(p[1][0] == a21);
2012 1 : CATCH_REQUIRE(p[1][1] == a22);
2013 1 : CATCH_REQUIRE(p[1][2] == a23);
2014 1 : CATCH_REQUIRE(p[2][0] == a31);
2015 1 : CATCH_REQUIRE(p[2][1] == a32);
2016 1 : CATCH_REQUIRE(p[2][2] == a33);
2017 1 : }
2018 :
2019 : {
2020 1 : snapdev::matrix<double> p = a.minor_matrix(0, 1);
2021 :
2022 1 : CATCH_REQUIRE(p.rows() == 3);
2023 1 : CATCH_REQUIRE(p.columns() == 3);
2024 :
2025 1 : CATCH_REQUIRE(p[0][0] == a10);
2026 1 : CATCH_REQUIRE(p[0][1] == a12);
2027 1 : CATCH_REQUIRE(p[0][2] == a13);
2028 1 : CATCH_REQUIRE(p[1][0] == a20);
2029 1 : CATCH_REQUIRE(p[1][1] == a22);
2030 1 : CATCH_REQUIRE(p[1][2] == a23);
2031 1 : CATCH_REQUIRE(p[2][0] == a30);
2032 1 : CATCH_REQUIRE(p[2][1] == a32);
2033 1 : CATCH_REQUIRE(p[2][2] == a33);
2034 1 : }
2035 :
2036 : {
2037 1 : snapdev::matrix<double> p = a.minor_matrix(0, 2);
2038 :
2039 1 : CATCH_REQUIRE(p.rows() == 3);
2040 1 : CATCH_REQUIRE(p.columns() == 3);
2041 :
2042 1 : CATCH_REQUIRE(p[0][0] == a10);
2043 1 : CATCH_REQUIRE(p[0][1] == a11);
2044 1 : CATCH_REQUIRE(p[0][2] == a13);
2045 1 : CATCH_REQUIRE(p[1][0] == a20);
2046 1 : CATCH_REQUIRE(p[1][1] == a21);
2047 1 : CATCH_REQUIRE(p[1][2] == a23);
2048 1 : CATCH_REQUIRE(p[2][0] == a30);
2049 1 : CATCH_REQUIRE(p[2][1] == a31);
2050 1 : CATCH_REQUIRE(p[2][2] == a33);
2051 1 : }
2052 :
2053 : {
2054 1 : snapdev::matrix<double> p = a.minor_matrix(0, 3);
2055 :
2056 1 : CATCH_REQUIRE(p.rows() == 3);
2057 1 : CATCH_REQUIRE(p.columns() == 3);
2058 :
2059 1 : CATCH_REQUIRE(p[0][0] == a10);
2060 1 : CATCH_REQUIRE(p[0][1] == a11);
2061 1 : CATCH_REQUIRE(p[0][2] == a12);
2062 1 : CATCH_REQUIRE(p[1][0] == a20);
2063 1 : CATCH_REQUIRE(p[1][1] == a21);
2064 1 : CATCH_REQUIRE(p[1][2] == a22);
2065 1 : CATCH_REQUIRE(p[2][0] == a30);
2066 1 : CATCH_REQUIRE(p[2][1] == a31);
2067 1 : CATCH_REQUIRE(p[2][2] == a32);
2068 1 : }
2069 :
2070 : {
2071 1 : snapdev::matrix<double> p = a.minor_matrix(2, 1);
2072 :
2073 1 : CATCH_REQUIRE(p.rows() == 3);
2074 1 : CATCH_REQUIRE(p.columns() == 3);
2075 :
2076 1 : CATCH_REQUIRE(p[0][0] == a00);
2077 1 : CATCH_REQUIRE(p[0][1] == a02);
2078 1 : CATCH_REQUIRE(p[0][2] == a03);
2079 1 : CATCH_REQUIRE(p[1][0] == a10);
2080 1 : CATCH_REQUIRE(p[1][1] == a12);
2081 1 : CATCH_REQUIRE(p[1][2] == a13);
2082 1 : CATCH_REQUIRE(p[2][0] == a30);
2083 1 : CATCH_REQUIRE(p[2][1] == a32);
2084 1 : CATCH_REQUIRE(p[2][2] == a33);
2085 1 : }
2086 1 : }
2087 5 : CATCH_END_SECTION()
2088 :
2089 5 : CATCH_START_SECTION("matrix: transpose 4x4")
2090 : {
2091 : // setup A
2092 : //
2093 1 : snapdev::matrix<double> a(4, 4);
2094 :
2095 1 : CATCH_REQUIRE(a.rows() == 4);
2096 1 : CATCH_REQUIRE(a.columns() == 4);
2097 :
2098 1 : double const a00(frand());
2099 1 : double const a01(frand());
2100 1 : double const a02(frand());
2101 1 : double const a03(frand());
2102 1 : double const a10(frand());
2103 1 : double const a11(frand());
2104 1 : double const a12(frand());
2105 1 : double const a13(frand());
2106 1 : double const a20(frand());
2107 1 : double const a21(frand());
2108 1 : double const a22(frand());
2109 1 : double const a23(frand());
2110 1 : double const a30(frand());
2111 1 : double const a31(frand());
2112 1 : double const a32(frand());
2113 1 : double const a33(frand());
2114 :
2115 1 : a[0][0] = a00;
2116 1 : a[0][1] = a01;
2117 1 : a[0][2] = a02;
2118 1 : a[0][3] = a03;
2119 1 : a[1][0] = a10;
2120 1 : a[1][1] = a11;
2121 1 : a[1][2] = a12;
2122 1 : a[1][3] = a13;
2123 1 : a[2][0] = a20;
2124 1 : a[2][1] = a21;
2125 1 : a[2][2] = a22;
2126 1 : a[2][3] = a23;
2127 1 : a[3][0] = a30;
2128 1 : a[3][1] = a31;
2129 1 : a[3][2] = a32;
2130 1 : a[3][3] = a33;
2131 :
2132 1 : CATCH_REQUIRE(a[0][0] == a00);
2133 1 : CATCH_REQUIRE(a[0][1] == a01);
2134 1 : CATCH_REQUIRE(a[0][2] == a02);
2135 1 : CATCH_REQUIRE(a[0][3] == a03);
2136 1 : CATCH_REQUIRE(a[1][0] == a10);
2137 1 : CATCH_REQUIRE(a[1][1] == a11);
2138 1 : CATCH_REQUIRE(a[1][2] == a12);
2139 1 : CATCH_REQUIRE(a[1][3] == a13);
2140 1 : CATCH_REQUIRE(a[2][0] == a20);
2141 1 : CATCH_REQUIRE(a[2][1] == a21);
2142 1 : CATCH_REQUIRE(a[2][2] == a22);
2143 1 : CATCH_REQUIRE(a[2][3] == a23);
2144 1 : CATCH_REQUIRE(a[3][0] == a30);
2145 1 : CATCH_REQUIRE(a[3][1] == a31);
2146 1 : CATCH_REQUIRE(a[3][2] == a32);
2147 1 : CATCH_REQUIRE(a[3][3] == a33);
2148 :
2149 1 : snapdev::matrix<double> t = a.transpose();
2150 :
2151 1 : CATCH_REQUIRE(t.rows() == 4);
2152 1 : CATCH_REQUIRE(t.columns() == 4);
2153 :
2154 1 : CATCH_REQUIRE(t[0][0] == a00);
2155 1 : CATCH_REQUIRE(t[0][1] == a10);
2156 1 : CATCH_REQUIRE(t[0][2] == a20);
2157 1 : CATCH_REQUIRE(t[0][3] == a30);
2158 1 : CATCH_REQUIRE(t[1][0] == a01);
2159 1 : CATCH_REQUIRE(t[1][1] == a11);
2160 1 : CATCH_REQUIRE(t[1][2] == a21);
2161 1 : CATCH_REQUIRE(t[1][3] == a31);
2162 1 : CATCH_REQUIRE(t[2][0] == a02);
2163 1 : CATCH_REQUIRE(t[2][1] == a12);
2164 1 : CATCH_REQUIRE(t[2][2] == a22);
2165 1 : CATCH_REQUIRE(t[2][3] == a32);
2166 1 : CATCH_REQUIRE(t[3][0] == a03);
2167 1 : CATCH_REQUIRE(t[3][1] == a13);
2168 1 : CATCH_REQUIRE(t[3][2] == a23);
2169 1 : CATCH_REQUIRE(t[3][3] == a33);
2170 1 : }
2171 5 : CATCH_END_SECTION()
2172 :
2173 5 : CATCH_START_SECTION("matrix: transpose 6x2")
2174 : {
2175 : // setup A
2176 : //
2177 1 : snapdev::matrix<double> a(6, 2);
2178 :
2179 1 : CATCH_REQUIRE(a.rows() == 6);
2180 1 : CATCH_REQUIRE(a.columns() == 2);
2181 :
2182 1 : double const a00(frand());
2183 1 : double const a01(frand());
2184 1 : double const a10(frand());
2185 1 : double const a11(frand());
2186 1 : double const a20(frand());
2187 1 : double const a21(frand());
2188 1 : double const a30(frand());
2189 1 : double const a31(frand());
2190 1 : double const a40(frand());
2191 1 : double const a41(frand());
2192 1 : double const a50(frand());
2193 1 : double const a51(frand());
2194 :
2195 1 : a[0][0] = a00;
2196 1 : a[0][1] = a01;
2197 1 : a[1][0] = a10;
2198 1 : a[1][1] = a11;
2199 1 : a[2][0] = a20;
2200 1 : a[2][1] = a21;
2201 1 : a[3][0] = a30;
2202 1 : a[3][1] = a31;
2203 1 : a[4][0] = a40;
2204 1 : a[4][1] = a41;
2205 1 : a[5][0] = a50;
2206 1 : a[5][1] = a51;
2207 :
2208 1 : CATCH_REQUIRE(a[0][0] == a00);
2209 1 : CATCH_REQUIRE(a[0][1] == a01);
2210 1 : CATCH_REQUIRE(a[1][0] == a10);
2211 1 : CATCH_REQUIRE(a[1][1] == a11);
2212 1 : CATCH_REQUIRE(a[2][0] == a20);
2213 1 : CATCH_REQUIRE(a[2][1] == a21);
2214 1 : CATCH_REQUIRE(a[3][0] == a30);
2215 1 : CATCH_REQUIRE(a[3][1] == a31);
2216 1 : CATCH_REQUIRE(a[4][0] == a40);
2217 1 : CATCH_REQUIRE(a[4][1] == a41);
2218 1 : CATCH_REQUIRE(a[5][0] == a50);
2219 1 : CATCH_REQUIRE(a[5][1] == a51);
2220 :
2221 1 : snapdev::matrix<double> t = a.transpose();
2222 :
2223 1 : CATCH_REQUIRE(t.rows() == 2);
2224 1 : CATCH_REQUIRE(t.columns() == 6);
2225 :
2226 : // original
2227 : // a00 a01
2228 : // a10 a11
2229 : // a20 a21
2230 : // a30 a31
2231 : // a40 a41
2232 : // a50 a51
2233 : //
2234 : // transposed
2235 : // a00 a10 a20 a30 a40 a50
2236 : // a01 a11 a21 a31 a41 a51
2237 : //
2238 :
2239 1 : CATCH_REQUIRE(t[0][0] == a00);
2240 1 : CATCH_REQUIRE(t[0][1] == a10);
2241 1 : CATCH_REQUIRE(t[0][2] == a20);
2242 1 : CATCH_REQUIRE(t[0][3] == a30);
2243 1 : CATCH_REQUIRE(t[0][4] == a40);
2244 1 : CATCH_REQUIRE(t[0][5] == a50);
2245 1 : CATCH_REQUIRE(t[1][0] == a01);
2246 1 : CATCH_REQUIRE(t[1][1] == a11);
2247 1 : CATCH_REQUIRE(t[1][2] == a21);
2248 1 : CATCH_REQUIRE(t[1][3] == a31);
2249 1 : CATCH_REQUIRE(t[1][4] == a41);
2250 1 : CATCH_REQUIRE(t[1][5] == a51);
2251 1 : }
2252 5 : CATCH_END_SECTION()
2253 :
2254 5 : CATCH_START_SECTION("matrix: adjugate 2x2")
2255 : {
2256 : // setup A
2257 : //
2258 1 : snapdev::matrix<double> a(2, 2);
2259 :
2260 1 : CATCH_REQUIRE(a.rows() == 2);
2261 1 : CATCH_REQUIRE(a.columns() == 2);
2262 :
2263 1 : double const a00(frand());
2264 1 : double const a01(frand());
2265 1 : double const a10(frand());
2266 1 : double const a11(frand());
2267 :
2268 1 : a[0][0] = a00;
2269 1 : a[0][1] = a01;
2270 1 : a[1][0] = a10;
2271 1 : a[1][1] = a11;
2272 :
2273 1 : CATCH_REQUIRE(a[0][0] == a00);
2274 1 : CATCH_REQUIRE(a[0][1] == a01);
2275 1 : CATCH_REQUIRE(a[1][0] == a10);
2276 1 : CATCH_REQUIRE(a[1][1] == a11);
2277 :
2278 1 : snapdev::matrix<double> m = a.adjugate();
2279 :
2280 1 : CATCH_REQUIRE(m.rows() == 2);
2281 1 : CATCH_REQUIRE(m.columns() == 2);
2282 :
2283 1 : CATCH_REQUIRE(m[0][0] == a11);
2284 1 : CATCH_REQUIRE(m[0][1] == -a01);
2285 1 : CATCH_REQUIRE(m[1][0] == -a10);
2286 1 : CATCH_REQUIRE(m[1][1] == a00);
2287 1 : }
2288 5 : CATCH_END_SECTION()
2289 :
2290 5 : CATCH_START_SECTION("matrix: adjugate 3x3")
2291 : {
2292 : // setup A
2293 : //
2294 1 : snapdev::matrix<double> a(3, 3);
2295 :
2296 1 : CATCH_REQUIRE(a.rows() == 3);
2297 1 : CATCH_REQUIRE(a.columns() == 3);
2298 :
2299 1 : double const a00(frand());
2300 1 : double const a01(frand());
2301 1 : double const a02(frand());
2302 1 : double const a10(frand());
2303 1 : double const a11(frand());
2304 1 : double const a12(frand());
2305 1 : double const a20(frand());
2306 1 : double const a21(frand());
2307 1 : double const a22(frand());
2308 :
2309 1 : a[0][0] = a00;
2310 1 : a[0][1] = a01;
2311 1 : a[0][2] = a02;
2312 1 : a[1][0] = a10;
2313 1 : a[1][1] = a11;
2314 1 : a[1][2] = a12;
2315 1 : a[2][0] = a20;
2316 1 : a[2][1] = a21;
2317 1 : a[2][2] = a22;
2318 :
2319 1 : CATCH_REQUIRE(a[0][0] == a00);
2320 1 : CATCH_REQUIRE(a[0][1] == a01);
2321 1 : CATCH_REQUIRE(a[0][2] == a02);
2322 1 : CATCH_REQUIRE(a[1][0] == a10);
2323 1 : CATCH_REQUIRE(a[1][1] == a11);
2324 1 : CATCH_REQUIRE(a[1][2] == a12);
2325 1 : CATCH_REQUIRE(a[2][0] == a20);
2326 1 : CATCH_REQUIRE(a[2][1] == a21);
2327 1 : CATCH_REQUIRE(a[2][2] == a22);
2328 :
2329 1 : snapdev::matrix<double> m = a.adjugate();
2330 :
2331 1 : CATCH_REQUIRE(m.rows() == 3);
2332 1 : CATCH_REQUIRE(m.columns() == 3);
2333 :
2334 : {
2335 1 : snapdev::matrix<double> p(a.minor_matrix(0, 0));
2336 1 : double const e(m[0][0] - p.determinant());
2337 1 : CATCH_REQUIRE(fabs(e) < 0.0001);
2338 1 : }
2339 :
2340 : {
2341 1 : snapdev::matrix<double> p(a.minor_matrix(0, 1));
2342 1 : double const e(m[1][0] + p.determinant());
2343 1 : CATCH_REQUIRE(fabs(e) < 0.0001);
2344 1 : }
2345 :
2346 : {
2347 1 : snapdev::matrix<double> p(a.minor_matrix(0, 2));
2348 1 : double const e(m[2][0] - p.determinant());
2349 1 : CATCH_REQUIRE(fabs(e) < 0.0001);
2350 1 : }
2351 :
2352 : {
2353 1 : snapdev::matrix<double> p(a.minor_matrix(1, 0));
2354 1 : double const e(m[0][1] + p.determinant());
2355 1 : CATCH_REQUIRE(fabs(e) < 0.0001);
2356 1 : }
2357 :
2358 : {
2359 1 : snapdev::matrix<double> p(a.minor_matrix(1, 1));
2360 1 : double const e(m[1][1] - p.determinant());
2361 1 : CATCH_REQUIRE(fabs(e) < 0.0001);
2362 1 : }
2363 :
2364 : {
2365 1 : snapdev::matrix<double> p(a.minor_matrix(1, 2));
2366 1 : double const e(m[2][1] + p.determinant());
2367 1 : CATCH_REQUIRE(fabs(e) < 0.0001);
2368 1 : }
2369 :
2370 : {
2371 1 : snapdev::matrix<double> p(a.minor_matrix(2, 0));
2372 1 : double const e(m[0][2] - p.determinant());
2373 1 : CATCH_REQUIRE(fabs(e) < 0.0001);
2374 1 : }
2375 :
2376 : {
2377 1 : snapdev::matrix<double> p(a.minor_matrix(2, 1));
2378 1 : double const e(m[1][2] + p.determinant());
2379 1 : CATCH_REQUIRE(fabs(e) < 0.0001);
2380 1 : }
2381 :
2382 : {
2383 1 : snapdev::matrix<double> p(a.minor_matrix(2, 2));
2384 1 : double const e(m[2][2] - p.determinant());
2385 1 : CATCH_REQUIRE(fabs(e) < 0.0001);
2386 1 : }
2387 1 : }
2388 5 : CATCH_END_SECTION()
2389 8 : }
2390 :
2391 : // calculations such as the determinant
2392 : //
2393 8 : CATCH_GIVEN("calculations")
2394 : {
2395 3 : CATCH_START_SECTION("matrix: 2x2 determinant")
2396 : {
2397 : // setup A
2398 : //
2399 1 : snapdev::matrix<double> a(2, 2);
2400 :
2401 1 : CATCH_REQUIRE(a.rows() == 2);
2402 1 : CATCH_REQUIRE(a.columns() == 2);
2403 :
2404 1 : double const a00(frand());
2405 1 : double const a01(frand());
2406 1 : double const a10(frand());
2407 1 : double const a11(frand());
2408 :
2409 1 : a[0][0] = a00;
2410 1 : a[0][1] = a01;
2411 1 : a[1][0] = a10;
2412 1 : a[1][1] = a11;
2413 :
2414 1 : CATCH_REQUIRE(a[0][0] == a00);
2415 1 : CATCH_REQUIRE(a[0][1] == a01);
2416 1 : CATCH_REQUIRE(a[1][0] == a10);
2417 1 : CATCH_REQUIRE(a[1][1] == a11);
2418 :
2419 1 : double d(a.determinant());
2420 :
2421 : //return f_vector[0 + 0 * 2] * f_vector[1 + 1 * 2]
2422 : // - f_vector[1 + 0 * 2] * f_vector[0 + 1 * 2];
2423 1 : double e = a00 * a11 - a10 * a01;
2424 :
2425 1 : CATCH_REQUIRE(fabs(d - e) < 0.0001);
2426 1 : }
2427 3 : CATCH_END_SECTION()
2428 :
2429 3 : CATCH_START_SECTION("matrix: 3x3 determinant with specific data")
2430 : {
2431 : // setup A
2432 : //
2433 1 : snapdev::matrix<double> a(3, 3);
2434 :
2435 1 : CATCH_REQUIRE(a.rows() == 3);
2436 1 : CATCH_REQUIRE(a.columns() == 3);
2437 :
2438 1 : a[0][0] = 5.0;
2439 1 : a[0][1] = -2.0;
2440 1 : a[0][2] = 1.0;
2441 1 : a[1][0] = 0.0;
2442 1 : a[1][1] = 3.0;
2443 1 : a[1][2] = -1.0;
2444 1 : a[2][0] = 2.0;
2445 1 : a[2][1] = 0.0;
2446 1 : a[2][2] = 7.0;
2447 :
2448 1 : double d(a.determinant());
2449 :
2450 : // we know that the exact answer for this one is 103
2451 : //
2452 1 : CATCH_REQUIRE(fabs(d - 103.0) < 0.0001);
2453 1 : }
2454 3 : CATCH_END_SECTION()
2455 :
2456 3 : CATCH_START_SECTION("matrix: 3x3 determinant with random data")
2457 : {
2458 : // verify 10 times (with different values)
2459 : //
2460 11 : for(int repeat(0); repeat < 10; ++repeat)
2461 : {
2462 : // setup A
2463 : //
2464 10 : snapdev::matrix<double> a(3, 3);
2465 :
2466 10 : CATCH_REQUIRE(a.rows() == 3);
2467 10 : CATCH_REQUIRE(a.columns() == 3);
2468 :
2469 10 : double const a00(frand());
2470 10 : double const a01(frand());
2471 10 : double const a02(frand());
2472 10 : double const a10(frand());
2473 10 : double const a11(frand());
2474 10 : double const a12(frand());
2475 10 : double const a20(frand());
2476 10 : double const a21(frand());
2477 10 : double const a22(frand());
2478 :
2479 10 : a[0][0] = a00;
2480 10 : a[0][1] = a01;
2481 10 : a[0][2] = a02;
2482 10 : a[1][0] = a10;
2483 10 : a[1][1] = a11;
2484 10 : a[1][2] = a12;
2485 10 : a[2][0] = a20;
2486 10 : a[2][1] = a21;
2487 10 : a[2][2] = a22;
2488 :
2489 10 : CATCH_REQUIRE(a[0][0] == a00);
2490 10 : CATCH_REQUIRE(a[0][1] == a01);
2491 10 : CATCH_REQUIRE(a[0][2] == a02);
2492 10 : CATCH_REQUIRE(a[1][0] == a10);
2493 10 : CATCH_REQUIRE(a[1][1] == a11);
2494 10 : CATCH_REQUIRE(a[1][2] == a12);
2495 10 : CATCH_REQUIRE(a[2][0] == a20);
2496 10 : CATCH_REQUIRE(a[2][1] == a21);
2497 10 : CATCH_REQUIRE(a[2][2] == a22);
2498 :
2499 10 : double d(a.determinant());
2500 :
2501 : //return f_vector[0 + 0 * 2] * f_vector[1 + 1 * 2]
2502 : // - f_vector[1 + 0 * 2] * f_vector[0 + 1 * 2];
2503 10 : double e = a00 * a11 * a22 + a01 * a12 * a20 + a02 * a10 * a21
2504 10 : - a00 * a12 * a21 - a01 * a10 * a22 - a02 * a11 * a20;
2505 :
2506 10 : CATCH_REQUIRE(fabs(d - e) < 0.0001);
2507 10 : }
2508 : }
2509 3 : CATCH_END_SECTION()
2510 8 : }
2511 8 : }
2512 :
2513 :
2514 8 : CATCH_TEST_CASE("matrix_multiplicative", "[matrix][math]")
2515 : {
2516 : // create two random 4x4 matrices and make sure the add works
2517 : //
2518 8 : CATCH_GIVEN("mul")
2519 : {
2520 4 : CATCH_START_SECTION("matrix: b=a*<scalar>")
2521 : {
2522 : // setup A
2523 : //
2524 1 : snapdev::matrix<double> a(4, 4);
2525 :
2526 1 : CATCH_REQUIRE(a.rows() == 4);
2527 1 : CATCH_REQUIRE(a.columns() == 4);
2528 :
2529 1 : double const a00(frand());
2530 1 : double const a01(frand());
2531 1 : double const a02(frand());
2532 1 : double const a03(frand());
2533 1 : double const a10(frand());
2534 1 : double const a11(frand());
2535 1 : double const a12(frand());
2536 1 : double const a13(frand());
2537 1 : double const a20(frand());
2538 1 : double const a21(frand());
2539 1 : double const a22(frand());
2540 1 : double const a23(frand());
2541 1 : double const a30(frand());
2542 1 : double const a31(frand());
2543 1 : double const a32(frand());
2544 1 : double const a33(frand());
2545 :
2546 1 : a[0][0] = a00;
2547 1 : a[0][1] = a01;
2548 1 : a[0][2] = a02;
2549 1 : a[0][3] = a03;
2550 1 : a[1][0] = a10;
2551 1 : a[1][1] = a11;
2552 1 : a[1][2] = a12;
2553 1 : a[1][3] = a13;
2554 1 : a[2][0] = a20;
2555 1 : a[2][1] = a21;
2556 1 : a[2][2] = a22;
2557 1 : a[2][3] = a23;
2558 1 : a[3][0] = a30;
2559 1 : a[3][1] = a31;
2560 1 : a[3][2] = a32;
2561 1 : a[3][3] = a33;
2562 :
2563 1 : CATCH_REQUIRE(a[0][0] == a00);
2564 1 : CATCH_REQUIRE(a[0][1] == a01);
2565 1 : CATCH_REQUIRE(a[0][2] == a02);
2566 1 : CATCH_REQUIRE(a[0][3] == a03);
2567 1 : CATCH_REQUIRE(a[1][0] == a10);
2568 1 : CATCH_REQUIRE(a[1][1] == a11);
2569 1 : CATCH_REQUIRE(a[1][2] == a12);
2570 1 : CATCH_REQUIRE(a[1][3] == a13);
2571 1 : CATCH_REQUIRE(a[2][0] == a20);
2572 1 : CATCH_REQUIRE(a[2][1] == a21);
2573 1 : CATCH_REQUIRE(a[2][2] == a22);
2574 1 : CATCH_REQUIRE(a[2][3] == a23);
2575 1 : CATCH_REQUIRE(a[3][0] == a30);
2576 1 : CATCH_REQUIRE(a[3][1] == a31);
2577 1 : CATCH_REQUIRE(a[3][2] == a32);
2578 1 : CATCH_REQUIRE(a[3][3] == a33);
2579 :
2580 : // setup B
2581 : //
2582 1 : snapdev::matrix<double> b(4, 4);
2583 :
2584 1 : CATCH_REQUIRE(b.rows() == 4);
2585 1 : CATCH_REQUIRE(b.columns() == 4);
2586 :
2587 1 : double const b00(frand());
2588 1 : double const b01(frand());
2589 1 : double const b02(frand());
2590 1 : double const b03(frand());
2591 1 : double const b10(frand());
2592 1 : double const b11(frand());
2593 1 : double const b12(frand());
2594 1 : double const b13(frand());
2595 1 : double const b20(frand());
2596 1 : double const b21(frand());
2597 1 : double const b22(frand());
2598 1 : double const b23(frand());
2599 1 : double const b30(frand());
2600 1 : double const b31(frand());
2601 1 : double const b32(frand());
2602 1 : double const b33(frand());
2603 :
2604 1 : b[0][0] = b00;
2605 1 : b[0][1] = b01;
2606 1 : b[0][2] = b02;
2607 1 : b[0][3] = b03;
2608 1 : b[1][0] = b10;
2609 1 : b[1][1] = b11;
2610 1 : b[1][2] = b12;
2611 1 : b[1][3] = b13;
2612 1 : b[2][0] = b20;
2613 1 : b[2][1] = b21;
2614 1 : b[2][2] = b22;
2615 1 : b[2][3] = b23;
2616 1 : b[3][0] = b30;
2617 1 : b[3][1] = b31;
2618 1 : b[3][2] = b32;
2619 1 : b[3][3] = b33;
2620 :
2621 1 : CATCH_REQUIRE(b[0][0] == b00);
2622 1 : CATCH_REQUIRE(b[0][1] == b01);
2623 1 : CATCH_REQUIRE(b[0][2] == b02);
2624 1 : CATCH_REQUIRE(b[0][3] == b03);
2625 1 : CATCH_REQUIRE(b[1][0] == b10);
2626 1 : CATCH_REQUIRE(b[1][1] == b11);
2627 1 : CATCH_REQUIRE(b[1][2] == b12);
2628 1 : CATCH_REQUIRE(b[1][3] == b13);
2629 1 : CATCH_REQUIRE(b[2][0] == b20);
2630 1 : CATCH_REQUIRE(b[2][1] == b21);
2631 1 : CATCH_REQUIRE(b[2][2] == b22);
2632 1 : CATCH_REQUIRE(b[2][3] == b23);
2633 1 : CATCH_REQUIRE(b[3][0] == b30);
2634 1 : CATCH_REQUIRE(b[3][1] == b31);
2635 1 : CATCH_REQUIRE(b[3][2] == b32);
2636 1 : CATCH_REQUIRE(b[3][3] == b33);
2637 :
2638 : // create a scalar for our test
2639 : //
2640 1 : double scalar(frand());
2641 1 : while(scalar == 0.0 || scalar == 1.0)
2642 : {
2643 0 : scalar = frand();
2644 : }
2645 :
2646 : // run operation B = A + <scalar>
2647 : //
2648 1 : b = a * scalar;
2649 :
2650 1 : CATCH_REQUIRE(a[0][0] == a00);
2651 1 : CATCH_REQUIRE(a[0][1] == a01);
2652 1 : CATCH_REQUIRE(a[0][2] == a02);
2653 1 : CATCH_REQUIRE(a[0][3] == a03);
2654 1 : CATCH_REQUIRE(a[1][0] == a10);
2655 1 : CATCH_REQUIRE(a[1][1] == a11);
2656 1 : CATCH_REQUIRE(a[1][2] == a12);
2657 1 : CATCH_REQUIRE(a[1][3] == a13);
2658 1 : CATCH_REQUIRE(a[2][0] == a20);
2659 1 : CATCH_REQUIRE(a[2][1] == a21);
2660 1 : CATCH_REQUIRE(a[2][2] == a22);
2661 1 : CATCH_REQUIRE(a[2][3] == a23);
2662 1 : CATCH_REQUIRE(a[3][0] == a30);
2663 1 : CATCH_REQUIRE(a[3][1] == a31);
2664 1 : CATCH_REQUIRE(a[3][2] == a32);
2665 1 : CATCH_REQUIRE(a[3][3] == a33);
2666 :
2667 : // this can't fail because we ensure scalar != 0 or 1
2668 : //
2669 1 : CATCH_REQUIRE_FALSE(b[0][0] == b00);
2670 1 : CATCH_REQUIRE_FALSE(b[0][1] == b01);
2671 1 : CATCH_REQUIRE_FALSE(b[0][2] == b02);
2672 1 : CATCH_REQUIRE_FALSE(b[0][3] == b03);
2673 1 : CATCH_REQUIRE_FALSE(b[1][0] == b10);
2674 1 : CATCH_REQUIRE_FALSE(b[1][1] == b11);
2675 1 : CATCH_REQUIRE_FALSE(b[1][2] == b12);
2676 1 : CATCH_REQUIRE_FALSE(b[1][3] == b13);
2677 1 : CATCH_REQUIRE_FALSE(b[2][0] == b20);
2678 1 : CATCH_REQUIRE_FALSE(b[2][1] == b21);
2679 1 : CATCH_REQUIRE_FALSE(b[2][2] == b22);
2680 1 : CATCH_REQUIRE_FALSE(b[2][3] == b23);
2681 1 : CATCH_REQUIRE_FALSE(b[3][0] == b30);
2682 1 : CATCH_REQUIRE_FALSE(b[3][1] == b31);
2683 1 : CATCH_REQUIRE_FALSE(b[3][2] == b32);
2684 1 : CATCH_REQUIRE_FALSE(b[3][3] == b33);
2685 :
2686 1 : CATCH_REQUIRE(b[0][0] == a00 * scalar);
2687 1 : CATCH_REQUIRE(b[0][1] == a01 * scalar);
2688 1 : CATCH_REQUIRE(b[0][2] == a02 * scalar);
2689 1 : CATCH_REQUIRE(b[0][3] == a03 * scalar);
2690 1 : CATCH_REQUIRE(b[1][0] == a10 * scalar);
2691 1 : CATCH_REQUIRE(b[1][1] == a11 * scalar);
2692 1 : CATCH_REQUIRE(b[1][2] == a12 * scalar);
2693 1 : CATCH_REQUIRE(b[1][3] == a13 * scalar);
2694 1 : CATCH_REQUIRE(b[2][0] == a20 * scalar);
2695 1 : CATCH_REQUIRE(b[2][1] == a21 * scalar);
2696 1 : CATCH_REQUIRE(b[2][2] == a22 * scalar);
2697 1 : CATCH_REQUIRE(b[2][3] == a23 * scalar);
2698 1 : CATCH_REQUIRE(b[3][0] == a30 * scalar);
2699 1 : CATCH_REQUIRE(b[3][1] == a31 * scalar);
2700 1 : CATCH_REQUIRE(b[3][2] == a32 * scalar);
2701 1 : CATCH_REQUIRE(b[3][3] == a33 * scalar);
2702 1 : }
2703 4 : CATCH_END_SECTION()
2704 :
2705 4 : CATCH_START_SECTION("matrix: a*=<scalar>")
2706 : {
2707 : // setup A
2708 : //
2709 1 : snapdev::matrix<double> a(4, 4);
2710 :
2711 1 : CATCH_REQUIRE(a.rows() == 4);
2712 1 : CATCH_REQUIRE(a.columns() == 4);
2713 :
2714 1 : double const a00(frand());
2715 1 : double const a01(frand());
2716 1 : double const a02(frand());
2717 1 : double const a03(frand());
2718 1 : double const a10(frand());
2719 1 : double const a11(frand());
2720 1 : double const a12(frand());
2721 1 : double const a13(frand());
2722 1 : double const a20(frand());
2723 1 : double const a21(frand());
2724 1 : double const a22(frand());
2725 1 : double const a23(frand());
2726 1 : double const a30(frand());
2727 1 : double const a31(frand());
2728 1 : double const a32(frand());
2729 1 : double const a33(frand());
2730 :
2731 1 : a[0][0] = a00;
2732 1 : a[0][1] = a01;
2733 1 : a[0][2] = a02;
2734 1 : a[0][3] = a03;
2735 1 : a[1][0] = a10;
2736 1 : a[1][1] = a11;
2737 1 : a[1][2] = a12;
2738 1 : a[1][3] = a13;
2739 1 : a[2][0] = a20;
2740 1 : a[2][1] = a21;
2741 1 : a[2][2] = a22;
2742 1 : a[2][3] = a23;
2743 1 : a[3][0] = a30;
2744 1 : a[3][1] = a31;
2745 1 : a[3][2] = a32;
2746 1 : a[3][3] = a33;
2747 :
2748 1 : CATCH_REQUIRE(a[0][0] == a00);
2749 1 : CATCH_REQUIRE(a[0][1] == a01);
2750 1 : CATCH_REQUIRE(a[0][2] == a02);
2751 1 : CATCH_REQUIRE(a[0][3] == a03);
2752 1 : CATCH_REQUIRE(a[1][0] == a10);
2753 1 : CATCH_REQUIRE(a[1][1] == a11);
2754 1 : CATCH_REQUIRE(a[1][2] == a12);
2755 1 : CATCH_REQUIRE(a[1][3] == a13);
2756 1 : CATCH_REQUIRE(a[2][0] == a20);
2757 1 : CATCH_REQUIRE(a[2][1] == a21);
2758 1 : CATCH_REQUIRE(a[2][2] == a22);
2759 1 : CATCH_REQUIRE(a[2][3] == a23);
2760 1 : CATCH_REQUIRE(a[3][0] == a30);
2761 1 : CATCH_REQUIRE(a[3][1] == a31);
2762 1 : CATCH_REQUIRE(a[3][2] == a32);
2763 1 : CATCH_REQUIRE(a[3][3] == a33);
2764 :
2765 : // create a scalar for our test
2766 : //
2767 1 : double scalar(frand());
2768 1 : while(scalar == 0.0 || scalar == 1.0)
2769 : {
2770 0 : scalar = frand();
2771 : }
2772 :
2773 : // run operation A *= <scalar>
2774 : //
2775 1 : a *= scalar;
2776 :
2777 : // this can't fail because we ensure scalar != 0 or 1
2778 : //
2779 1 : CATCH_REQUIRE_FALSE(a[0][0] == a00);
2780 1 : CATCH_REQUIRE_FALSE(a[0][1] == a01);
2781 1 : CATCH_REQUIRE_FALSE(a[0][2] == a02);
2782 1 : CATCH_REQUIRE_FALSE(a[0][3] == a03);
2783 1 : CATCH_REQUIRE_FALSE(a[1][0] == a10);
2784 1 : CATCH_REQUIRE_FALSE(a[1][1] == a11);
2785 1 : CATCH_REQUIRE_FALSE(a[1][2] == a12);
2786 1 : CATCH_REQUIRE_FALSE(a[1][3] == a13);
2787 1 : CATCH_REQUIRE_FALSE(a[2][0] == a20);
2788 1 : CATCH_REQUIRE_FALSE(a[2][1] == a21);
2789 1 : CATCH_REQUIRE_FALSE(a[2][2] == a22);
2790 1 : CATCH_REQUIRE_FALSE(a[2][3] == a23);
2791 1 : CATCH_REQUIRE_FALSE(a[3][0] == a30);
2792 1 : CATCH_REQUIRE_FALSE(a[3][1] == a31);
2793 1 : CATCH_REQUIRE_FALSE(a[3][2] == a32);
2794 1 : CATCH_REQUIRE_FALSE(a[3][3] == a33);
2795 :
2796 1 : CATCH_REQUIRE(a[0][0] == a00 * scalar);
2797 1 : CATCH_REQUIRE(a[0][1] == a01 * scalar);
2798 1 : CATCH_REQUIRE(a[0][2] == a02 * scalar);
2799 1 : CATCH_REQUIRE(a[0][3] == a03 * scalar);
2800 1 : CATCH_REQUIRE(a[1][0] == a10 * scalar);
2801 1 : CATCH_REQUIRE(a[1][1] == a11 * scalar);
2802 1 : CATCH_REQUIRE(a[1][2] == a12 * scalar);
2803 1 : CATCH_REQUIRE(a[1][3] == a13 * scalar);
2804 1 : CATCH_REQUIRE(a[2][0] == a20 * scalar);
2805 1 : CATCH_REQUIRE(a[2][1] == a21 * scalar);
2806 1 : CATCH_REQUIRE(a[2][2] == a22 * scalar);
2807 1 : CATCH_REQUIRE(a[2][3] == a23 * scalar);
2808 1 : CATCH_REQUIRE(a[3][0] == a30 * scalar);
2809 1 : CATCH_REQUIRE(a[3][1] == a31 * scalar);
2810 1 : CATCH_REQUIRE(a[3][2] == a32 * scalar);
2811 1 : CATCH_REQUIRE(a[3][3] == a33 * scalar);
2812 1 : }
2813 4 : CATCH_END_SECTION()
2814 :
2815 4 : CATCH_START_SECTION("matrix: c=a*b")
2816 : {
2817 : // setup A
2818 : //
2819 1 : snapdev::matrix<double> a(4, 4);
2820 :
2821 1 : CATCH_REQUIRE(a.rows() == 4);
2822 1 : CATCH_REQUIRE(a.columns() == 4);
2823 :
2824 1 : double const a00(frand());
2825 1 : double const a01(frand());
2826 1 : double const a02(frand());
2827 1 : double const a03(frand());
2828 1 : double const a10(frand());
2829 1 : double const a11(frand());
2830 1 : double const a12(frand());
2831 1 : double const a13(frand());
2832 1 : double const a20(frand());
2833 1 : double const a21(frand());
2834 1 : double const a22(frand());
2835 1 : double const a23(frand());
2836 1 : double const a30(frand());
2837 1 : double const a31(frand());
2838 1 : double const a32(frand());
2839 1 : double const a33(frand());
2840 :
2841 1 : a[0][0] = a00;
2842 1 : a[0][1] = a01;
2843 1 : a[0][2] = a02;
2844 1 : a[0][3] = a03;
2845 1 : a[1][0] = a10;
2846 1 : a[1][1] = a11;
2847 1 : a[1][2] = a12;
2848 1 : a[1][3] = a13;
2849 1 : a[2][0] = a20;
2850 1 : a[2][1] = a21;
2851 1 : a[2][2] = a22;
2852 1 : a[2][3] = a23;
2853 1 : a[3][0] = a30;
2854 1 : a[3][1] = a31;
2855 1 : a[3][2] = a32;
2856 1 : a[3][3] = a33;
2857 :
2858 1 : CATCH_REQUIRE(a[0][0] == a00);
2859 1 : CATCH_REQUIRE(a[0][1] == a01);
2860 1 : CATCH_REQUIRE(a[0][2] == a02);
2861 1 : CATCH_REQUIRE(a[0][3] == a03);
2862 1 : CATCH_REQUIRE(a[1][0] == a10);
2863 1 : CATCH_REQUIRE(a[1][1] == a11);
2864 1 : CATCH_REQUIRE(a[1][2] == a12);
2865 1 : CATCH_REQUIRE(a[1][3] == a13);
2866 1 : CATCH_REQUIRE(a[2][0] == a20);
2867 1 : CATCH_REQUIRE(a[2][1] == a21);
2868 1 : CATCH_REQUIRE(a[2][2] == a22);
2869 1 : CATCH_REQUIRE(a[2][3] == a23);
2870 1 : CATCH_REQUIRE(a[3][0] == a30);
2871 1 : CATCH_REQUIRE(a[3][1] == a31);
2872 1 : CATCH_REQUIRE(a[3][2] == a32);
2873 1 : CATCH_REQUIRE(a[3][3] == a33);
2874 :
2875 : // setup B
2876 : //
2877 1 : snapdev::matrix<double> b(4, 4);
2878 :
2879 1 : CATCH_REQUIRE(b.rows() == 4);
2880 1 : CATCH_REQUIRE(b.columns() == 4);
2881 :
2882 1 : double const b00(frand());
2883 1 : double const b01(frand());
2884 1 : double const b02(frand());
2885 1 : double const b03(frand());
2886 1 : double const b10(frand());
2887 1 : double const b11(frand());
2888 1 : double const b12(frand());
2889 1 : double const b13(frand());
2890 1 : double const b20(frand());
2891 1 : double const b21(frand());
2892 1 : double const b22(frand());
2893 1 : double const b23(frand());
2894 1 : double const b30(frand());
2895 1 : double const b31(frand());
2896 1 : double const b32(frand());
2897 1 : double const b33(frand());
2898 :
2899 1 : b[0][0] = b00;
2900 1 : b[0][1] = b01;
2901 1 : b[0][2] = b02;
2902 1 : b[0][3] = b03;
2903 1 : b[1][0] = b10;
2904 1 : b[1][1] = b11;
2905 1 : b[1][2] = b12;
2906 1 : b[1][3] = b13;
2907 1 : b[2][0] = b20;
2908 1 : b[2][1] = b21;
2909 1 : b[2][2] = b22;
2910 1 : b[2][3] = b23;
2911 1 : b[3][0] = b30;
2912 1 : b[3][1] = b31;
2913 1 : b[3][2] = b32;
2914 1 : b[3][3] = b33;
2915 :
2916 1 : CATCH_REQUIRE(b[0][0] == b00);
2917 1 : CATCH_REQUIRE(b[0][1] == b01);
2918 1 : CATCH_REQUIRE(b[0][2] == b02);
2919 1 : CATCH_REQUIRE(b[0][3] == b03);
2920 1 : CATCH_REQUIRE(b[1][0] == b10);
2921 1 : CATCH_REQUIRE(b[1][1] == b11);
2922 1 : CATCH_REQUIRE(b[1][2] == b12);
2923 1 : CATCH_REQUIRE(b[1][3] == b13);
2924 1 : CATCH_REQUIRE(b[2][0] == b20);
2925 1 : CATCH_REQUIRE(b[2][1] == b21);
2926 1 : CATCH_REQUIRE(b[2][2] == b22);
2927 1 : CATCH_REQUIRE(b[2][3] == b23);
2928 1 : CATCH_REQUIRE(b[3][0] == b30);
2929 1 : CATCH_REQUIRE(b[3][1] == b31);
2930 1 : CATCH_REQUIRE(b[3][2] == b32);
2931 1 : CATCH_REQUIRE(b[3][3] == b33);
2932 :
2933 : // setup C
2934 : //
2935 1 : snapdev::matrix<double> c(4, 4);
2936 :
2937 1 : CATCH_REQUIRE(c.rows() == 4);
2938 1 : CATCH_REQUIRE(c.columns() == 4);
2939 :
2940 1 : double const c00(frand());
2941 1 : double const c01(frand());
2942 1 : double const c02(frand());
2943 1 : double const c03(frand());
2944 1 : double const c10(frand());
2945 1 : double const c11(frand());
2946 1 : double const c12(frand());
2947 1 : double const c13(frand());
2948 1 : double const c20(frand());
2949 1 : double const c21(frand());
2950 1 : double const c22(frand());
2951 1 : double const c23(frand());
2952 1 : double const c30(frand());
2953 1 : double const c31(frand());
2954 1 : double const c32(frand());
2955 1 : double const c33(frand());
2956 :
2957 1 : c[0][0] = c00;
2958 1 : c[0][1] = c01;
2959 1 : c[0][2] = c02;
2960 1 : c[0][3] = c03;
2961 1 : c[1][0] = c10;
2962 1 : c[1][1] = c11;
2963 1 : c[1][2] = c12;
2964 1 : c[1][3] = c13;
2965 1 : c[2][0] = c20;
2966 1 : c[2][1] = c21;
2967 1 : c[2][2] = c22;
2968 1 : c[2][3] = c23;
2969 1 : c[3][0] = c30;
2970 1 : c[3][1] = c31;
2971 1 : c[3][2] = c32;
2972 1 : c[3][3] = c33;
2973 :
2974 1 : CATCH_REQUIRE(c[0][0] == c00);
2975 1 : CATCH_REQUIRE(c[0][1] == c01);
2976 1 : CATCH_REQUIRE(c[0][2] == c02);
2977 1 : CATCH_REQUIRE(c[0][3] == c03);
2978 1 : CATCH_REQUIRE(c[1][0] == c10);
2979 1 : CATCH_REQUIRE(c[1][1] == c11);
2980 1 : CATCH_REQUIRE(c[1][2] == c12);
2981 1 : CATCH_REQUIRE(c[1][3] == c13);
2982 1 : CATCH_REQUIRE(c[2][0] == c20);
2983 1 : CATCH_REQUIRE(c[2][1] == c21);
2984 1 : CATCH_REQUIRE(c[2][2] == c22);
2985 1 : CATCH_REQUIRE(c[2][3] == c23);
2986 1 : CATCH_REQUIRE(c[3][0] == c30);
2987 1 : CATCH_REQUIRE(c[3][1] == c31);
2988 1 : CATCH_REQUIRE(c[3][2] == c32);
2989 1 : CATCH_REQUIRE(c[3][3] == c33);
2990 :
2991 : // run operation C = A * B
2992 : //
2993 1 : c = a * b;
2994 :
2995 1 : CATCH_REQUIRE(a[0][0] == a00);
2996 1 : CATCH_REQUIRE(a[0][1] == a01);
2997 1 : CATCH_REQUIRE(a[0][2] == a02);
2998 1 : CATCH_REQUIRE(a[0][3] == a03);
2999 1 : CATCH_REQUIRE(a[1][0] == a10);
3000 1 : CATCH_REQUIRE(a[1][1] == a11);
3001 1 : CATCH_REQUIRE(a[1][2] == a12);
3002 1 : CATCH_REQUIRE(a[1][3] == a13);
3003 1 : CATCH_REQUIRE(a[2][0] == a20);
3004 1 : CATCH_REQUIRE(a[2][1] == a21);
3005 1 : CATCH_REQUIRE(a[2][2] == a22);
3006 1 : CATCH_REQUIRE(a[2][3] == a23);
3007 1 : CATCH_REQUIRE(a[3][0] == a30);
3008 1 : CATCH_REQUIRE(a[3][1] == a31);
3009 1 : CATCH_REQUIRE(a[3][2] == a32);
3010 1 : CATCH_REQUIRE(a[3][3] == a33);
3011 :
3012 1 : CATCH_REQUIRE(b[0][0] == b00);
3013 1 : CATCH_REQUIRE(b[0][1] == b01);
3014 1 : CATCH_REQUIRE(b[0][2] == b02);
3015 1 : CATCH_REQUIRE(b[0][3] == b03);
3016 1 : CATCH_REQUIRE(b[1][0] == b10);
3017 1 : CATCH_REQUIRE(b[1][1] == b11);
3018 1 : CATCH_REQUIRE(b[1][2] == b12);
3019 1 : CATCH_REQUIRE(b[1][3] == b13);
3020 1 : CATCH_REQUIRE(b[2][0] == b20);
3021 1 : CATCH_REQUIRE(b[2][1] == b21);
3022 1 : CATCH_REQUIRE(b[2][2] == b22);
3023 1 : CATCH_REQUIRE(b[2][3] == b23);
3024 1 : CATCH_REQUIRE(b[3][0] == b30);
3025 1 : CATCH_REQUIRE(b[3][1] == b31);
3026 1 : CATCH_REQUIRE(b[3][2] == b32);
3027 1 : CATCH_REQUIRE(b[3][3] == b33);
3028 :
3029 : // this could fail (albeit rather unlikely)
3030 : //CATCH_REQUIRE_FALSE(c[0][0] == c00);
3031 : //CATCH_REQUIRE_FALSE(c[0][1] == c01);
3032 : //CATCH_REQUIRE_FALSE(c[0][2] == c02);
3033 : //CATCH_REQUIRE_FALSE(c[0][3] == c03);
3034 : //CATCH_REQUIRE_FALSE(c[1][0] == c10);
3035 : //CATCH_REQUIRE_FALSE(c[1][1] == c11);
3036 : //CATCH_REQUIRE_FALSE(c[1][2] == c12);
3037 : //CATCH_REQUIRE_FALSE(c[1][3] == c13);
3038 : //CATCH_REQUIRE_FALSE(c[2][0] == c20);
3039 : //CATCH_REQUIRE_FALSE(c[2][1] == c21);
3040 : //CATCH_REQUIRE_FALSE(c[2][2] == c22);
3041 : //CATCH_REQUIRE_FALSE(c[2][3] == c23);
3042 : //CATCH_REQUIRE_FALSE(c[3][0] == c30);
3043 : //CATCH_REQUIRE_FALSE(c[3][1] == c31);
3044 : //CATCH_REQUIRE_FALSE(c[3][2] == c32);
3045 : //CATCH_REQUIRE_FALSE(c[3][3] == c33);
3046 :
3047 : //std::cout << "A = " << a << std::endl;
3048 : //std::cout << "B = " << b << std::endl;
3049 : //std::cout << "C = " << c << std::endl;
3050 :
3051 1 : CATCH_REQUIRE(c[0][0] == a00 * b00 + a01 * b10 + a02 * b20 + a03 * b30);
3052 1 : CATCH_REQUIRE(c[0][1] == a00 * b01 + a01 * b11 + a02 * b21 + a03 * b31);
3053 1 : CATCH_REQUIRE(c[0][2] == a00 * b02 + a01 * b12 + a02 * b22 + a03 * b32);
3054 1 : CATCH_REQUIRE(c[0][3] == a00 * b03 + a01 * b13 + a02 * b23 + a03 * b33);
3055 :
3056 1 : CATCH_REQUIRE(c[1][0] == a10 * b00 + a11 * b10 + a12 * b20 + a13 * b30);
3057 1 : CATCH_REQUIRE(c[1][1] == a10 * b01 + a11 * b11 + a12 * b21 + a13 * b31);
3058 1 : CATCH_REQUIRE(c[1][2] == a10 * b02 + a11 * b12 + a12 * b22 + a13 * b32);
3059 1 : CATCH_REQUIRE(c[1][3] == a10 * b03 + a11 * b13 + a12 * b23 + a13 * b33);
3060 :
3061 1 : CATCH_REQUIRE(c[2][0] == a20 * b00 + a21 * b10 + a22 * b20 + a23 * b30);
3062 1 : CATCH_REQUIRE(c[2][1] == a20 * b01 + a21 * b11 + a22 * b21 + a23 * b31);
3063 1 : CATCH_REQUIRE(c[2][2] == a20 * b02 + a21 * b12 + a22 * b22 + a23 * b32);
3064 1 : CATCH_REQUIRE(c[2][3] == a20 * b03 + a21 * b13 + a22 * b23 + a23 * b33);
3065 :
3066 1 : CATCH_REQUIRE(c[3][0] == a30 * b00 + a31 * b10 + a32 * b20 + a33 * b30);
3067 1 : CATCH_REQUIRE(c[3][1] == a30 * b01 + a31 * b11 + a32 * b21 + a33 * b31);
3068 1 : CATCH_REQUIRE(c[3][2] == a30 * b02 + a31 * b12 + a32 * b22 + a33 * b32);
3069 1 : CATCH_REQUIRE(c[3][3] == a30 * b03 + a31 * b13 + a32 * b23 + a33 * b33);
3070 :
3071 : // if swapped it would equal this instead
3072 : //CATCH_REQUIRE(c[0][0] == a00 * b00 + a10 * b01 + a20 * b02 + a30 * b03);
3073 : //CATCH_REQUIRE(c[0][1] == a01 * b00 + a11 * b01 + a21 * b02 + a31 * b03);
3074 : //CATCH_REQUIRE(c[0][2] == a02 * b00 + a12 * b01 + a22 * b02 + a32 * b03);
3075 : //CATCH_REQUIRE(c[0][3] == a03 * b00 + a13 * b01 + a23 * b02 + a33 * b03);
3076 :
3077 : //CATCH_REQUIRE(c[1][0] == a00 * b10 + a10 * b11 + a20 * b12 + a30 * b13);
3078 : //CATCH_REQUIRE(c[1][1] == a01 * b10 + a11 * b11 + a21 * b12 + a31 * b13);
3079 : //CATCH_REQUIRE(c[1][2] == a02 * b10 + a12 * b11 + a22 * b12 + a32 * b13);
3080 : //CATCH_REQUIRE(c[1][3] == a03 * b10 + a13 * b11 + a23 * b12 + a33 * b13);
3081 :
3082 : //CATCH_REQUIRE(c[2][0] == a00 * b20 + a10 * b21 + a20 * b22 + a30 * b23);
3083 : //CATCH_REQUIRE(c[2][1] == a01 * b20 + a11 * b21 + a21 * b22 + a31 * b23);
3084 : //CATCH_REQUIRE(c[2][2] == a02 * b20 + a12 * b21 + a22 * b22 + a32 * b23);
3085 : //CATCH_REQUIRE(c[2][3] == a03 * b20 + a13 * b21 + a23 * b22 + a33 * b23);
3086 :
3087 : //CATCH_REQUIRE(c[3][0] == a00 * b30 + a10 * b31 + a20 * b32 + a30 * b33);
3088 : //CATCH_REQUIRE(c[3][1] == a01 * b30 + a11 * b31 + a21 * b32 + a31 * b33);
3089 : //CATCH_REQUIRE(c[3][2] == a02 * b30 + a12 * b31 + a22 * b32 + a32 * b33);
3090 : //CATCH_REQUIRE(c[3][3] == a03 * b30 + a13 * b31 + a23 * b32 + a33 * b33);
3091 1 : }
3092 4 : CATCH_END_SECTION()
3093 :
3094 4 : CATCH_START_SECTION("matrix: a*=b")
3095 : {
3096 : // setup A
3097 : //
3098 1 : snapdev::matrix<double> a(4, 4);
3099 :
3100 1 : CATCH_REQUIRE(a.rows() == 4);
3101 1 : CATCH_REQUIRE(a.columns() == 4);
3102 :
3103 1 : double const a00(frand());
3104 1 : double const a01(frand());
3105 1 : double const a02(frand());
3106 1 : double const a03(frand());
3107 1 : double const a10(frand());
3108 1 : double const a11(frand());
3109 1 : double const a12(frand());
3110 1 : double const a13(frand());
3111 1 : double const a20(frand());
3112 1 : double const a21(frand());
3113 1 : double const a22(frand());
3114 1 : double const a23(frand());
3115 1 : double const a30(frand());
3116 1 : double const a31(frand());
3117 1 : double const a32(frand());
3118 1 : double const a33(frand());
3119 :
3120 1 : a[0][0] = a00;
3121 1 : a[0][1] = a01;
3122 1 : a[0][2] = a02;
3123 1 : a[0][3] = a03;
3124 1 : a[1][0] = a10;
3125 1 : a[1][1] = a11;
3126 1 : a[1][2] = a12;
3127 1 : a[1][3] = a13;
3128 1 : a[2][0] = a20;
3129 1 : a[2][1] = a21;
3130 1 : a[2][2] = a22;
3131 1 : a[2][3] = a23;
3132 1 : a[3][0] = a30;
3133 1 : a[3][1] = a31;
3134 1 : a[3][2] = a32;
3135 1 : a[3][3] = a33;
3136 :
3137 1 : CATCH_REQUIRE(a[0][0] == a00);
3138 1 : CATCH_REQUIRE(a[0][1] == a01);
3139 1 : CATCH_REQUIRE(a[0][2] == a02);
3140 1 : CATCH_REQUIRE(a[0][3] == a03);
3141 1 : CATCH_REQUIRE(a[1][0] == a10);
3142 1 : CATCH_REQUIRE(a[1][1] == a11);
3143 1 : CATCH_REQUIRE(a[1][2] == a12);
3144 1 : CATCH_REQUIRE(a[1][3] == a13);
3145 1 : CATCH_REQUIRE(a[2][0] == a20);
3146 1 : CATCH_REQUIRE(a[2][1] == a21);
3147 1 : CATCH_REQUIRE(a[2][2] == a22);
3148 1 : CATCH_REQUIRE(a[2][3] == a23);
3149 1 : CATCH_REQUIRE(a[3][0] == a30);
3150 1 : CATCH_REQUIRE(a[3][1] == a31);
3151 1 : CATCH_REQUIRE(a[3][2] == a32);
3152 1 : CATCH_REQUIRE(a[3][3] == a33);
3153 :
3154 : // setup B
3155 : //
3156 1 : snapdev::matrix<double> b(4, 4);
3157 :
3158 1 : CATCH_REQUIRE(b.rows() == 4);
3159 1 : CATCH_REQUIRE(b.columns() == 4);
3160 :
3161 1 : double const b00(frand());
3162 1 : double const b01(frand());
3163 1 : double const b02(frand());
3164 1 : double const b03(frand());
3165 1 : double const b10(frand());
3166 1 : double const b11(frand());
3167 1 : double const b12(frand());
3168 1 : double const b13(frand());
3169 1 : double const b20(frand());
3170 1 : double const b21(frand());
3171 1 : double const b22(frand());
3172 1 : double const b23(frand());
3173 1 : double const b30(frand());
3174 1 : double const b31(frand());
3175 1 : double const b32(frand());
3176 1 : double const b33(frand());
3177 :
3178 1 : b[0][0] = b00;
3179 1 : b[0][1] = b01;
3180 1 : b[0][2] = b02;
3181 1 : b[0][3] = b03;
3182 1 : b[1][0] = b10;
3183 1 : b[1][1] = b11;
3184 1 : b[1][2] = b12;
3185 1 : b[1][3] = b13;
3186 1 : b[2][0] = b20;
3187 1 : b[2][1] = b21;
3188 1 : b[2][2] = b22;
3189 1 : b[2][3] = b23;
3190 1 : b[3][0] = b30;
3191 1 : b[3][1] = b31;
3192 1 : b[3][2] = b32;
3193 1 : b[3][3] = b33;
3194 :
3195 1 : CATCH_REQUIRE(b[0][0] == b00);
3196 1 : CATCH_REQUIRE(b[0][1] == b01);
3197 1 : CATCH_REQUIRE(b[0][2] == b02);
3198 1 : CATCH_REQUIRE(b[0][3] == b03);
3199 1 : CATCH_REQUIRE(b[1][0] == b10);
3200 1 : CATCH_REQUIRE(b[1][1] == b11);
3201 1 : CATCH_REQUIRE(b[1][2] == b12);
3202 1 : CATCH_REQUIRE(b[1][3] == b13);
3203 1 : CATCH_REQUIRE(b[2][0] == b20);
3204 1 : CATCH_REQUIRE(b[2][1] == b21);
3205 1 : CATCH_REQUIRE(b[2][2] == b22);
3206 1 : CATCH_REQUIRE(b[2][3] == b23);
3207 1 : CATCH_REQUIRE(b[3][0] == b30);
3208 1 : CATCH_REQUIRE(b[3][1] == b31);
3209 1 : CATCH_REQUIRE(b[3][2] == b32);
3210 1 : CATCH_REQUIRE(b[3][3] == b33);
3211 :
3212 : // run operation A *= B
3213 : //
3214 1 : a *= b;
3215 :
3216 : // this could fail
3217 : //
3218 : //CATCH_REQUIRE_FALSE(a[0][0] == a00);
3219 : //CATCH_REQUIRE_FALSE(a[0][1] == a01);
3220 : //CATCH_REQUIRE_FALSE(a[0][2] == a02);
3221 : //CATCH_REQUIRE_FALSE(a[0][3] == a03);
3222 : //CATCH_REQUIRE_FALSE(a[1][0] == a10);
3223 : //CATCH_REQUIRE_FALSE(a[1][1] == a11);
3224 : //CATCH_REQUIRE_FALSE(a[1][2] == a12);
3225 : //CATCH_REQUIRE_FALSE(a[1][3] == a13);
3226 : //CATCH_REQUIRE_FALSE(a[2][0] == a20);
3227 : //CATCH_REQUIRE_FALSE(a[2][1] == a21);
3228 : //CATCH_REQUIRE_FALSE(a[2][2] == a22);
3229 : //CATCH_REQUIRE_FALSE(a[2][3] == a23);
3230 : //CATCH_REQUIRE_FALSE(a[3][0] == a30);
3231 : //CATCH_REQUIRE_FALSE(a[3][1] == a31);
3232 : //CATCH_REQUIRE_FALSE(a[3][2] == a32);
3233 : //CATCH_REQUIRE_FALSE(a[3][3] == a33);
3234 :
3235 1 : CATCH_REQUIRE(b[0][0] == b00);
3236 1 : CATCH_REQUIRE(b[0][1] == b01);
3237 1 : CATCH_REQUIRE(b[0][2] == b02);
3238 1 : CATCH_REQUIRE(b[0][3] == b03);
3239 1 : CATCH_REQUIRE(b[1][0] == b10);
3240 1 : CATCH_REQUIRE(b[1][1] == b11);
3241 1 : CATCH_REQUIRE(b[1][2] == b12);
3242 1 : CATCH_REQUIRE(b[1][3] == b13);
3243 1 : CATCH_REQUIRE(b[2][0] == b20);
3244 1 : CATCH_REQUIRE(b[2][1] == b21);
3245 1 : CATCH_REQUIRE(b[2][2] == b22);
3246 1 : CATCH_REQUIRE(b[2][3] == b23);
3247 1 : CATCH_REQUIRE(b[3][0] == b30);
3248 1 : CATCH_REQUIRE(b[3][1] == b31);
3249 1 : CATCH_REQUIRE(b[3][2] == b32);
3250 1 : CATCH_REQUIRE(b[3][3] == b33);
3251 :
3252 1 : CATCH_REQUIRE(a[0][0] == a00 * b00 + a01 * b10 + a02 * b20 + a03 * b30);
3253 1 : CATCH_REQUIRE(a[0][1] == a00 * b01 + a01 * b11 + a02 * b21 + a03 * b31);
3254 1 : CATCH_REQUIRE(a[0][2] == a00 * b02 + a01 * b12 + a02 * b22 + a03 * b32);
3255 1 : CATCH_REQUIRE(a[0][3] == a00 * b03 + a01 * b13 + a02 * b23 + a03 * b33);
3256 :
3257 1 : CATCH_REQUIRE(a[1][0] == a10 * b00 + a11 * b10 + a12 * b20 + a13 * b30);
3258 1 : CATCH_REQUIRE(a[1][1] == a10 * b01 + a11 * b11 + a12 * b21 + a13 * b31);
3259 1 : CATCH_REQUIRE(a[1][2] == a10 * b02 + a11 * b12 + a12 * b22 + a13 * b32);
3260 1 : CATCH_REQUIRE(a[1][3] == a10 * b03 + a11 * b13 + a12 * b23 + a13 * b33);
3261 :
3262 1 : CATCH_REQUIRE(a[2][0] == a20 * b00 + a21 * b10 + a22 * b20 + a23 * b30);
3263 1 : CATCH_REQUIRE(a[2][1] == a20 * b01 + a21 * b11 + a22 * b21 + a23 * b31);
3264 1 : CATCH_REQUIRE(a[2][2] == a20 * b02 + a21 * b12 + a22 * b22 + a23 * b32);
3265 1 : CATCH_REQUIRE(a[2][3] == a20 * b03 + a21 * b13 + a22 * b23 + a23 * b33);
3266 :
3267 1 : CATCH_REQUIRE(a[3][0] == a30 * b00 + a31 * b10 + a32 * b20 + a33 * b30);
3268 1 : CATCH_REQUIRE(a[3][1] == a30 * b01 + a31 * b11 + a32 * b21 + a33 * b31);
3269 1 : CATCH_REQUIRE(a[3][2] == a30 * b02 + a31 * b12 + a32 * b22 + a33 * b32);
3270 1 : CATCH_REQUIRE(a[3][3] == a30 * b03 + a31 * b13 + a32 * b23 + a33 * b33);
3271 :
3272 : //CATCH_REQUIRE(a[0][0] == a00 * b00 + a10 * b01 + a20 * b02 + a30 * b03);
3273 : //CATCH_REQUIRE(a[0][1] == a01 * b00 + a11 * b01 + a21 * b02 + a31 * b03);
3274 : //CATCH_REQUIRE(a[0][2] == a02 * b00 + a12 * b01 + a22 * b02 + a32 * b03);
3275 : //CATCH_REQUIRE(a[0][3] == a03 * b00 + a13 * b01 + a23 * b02 + a33 * b03);
3276 :
3277 : //CATCH_REQUIRE(a[1][0] == a00 * b10 + a10 * b11 + a20 * b12 + a30 * b13);
3278 : //CATCH_REQUIRE(a[1][1] == a01 * b10 + a11 * b11 + a21 * b12 + a31 * b13);
3279 : //CATCH_REQUIRE(a[1][2] == a02 * b10 + a12 * b11 + a22 * b12 + a32 * b13);
3280 : //CATCH_REQUIRE(a[1][3] == a03 * b10 + a13 * b11 + a23 * b12 + a33 * b13);
3281 :
3282 : //CATCH_REQUIRE(a[2][0] == a00 * b20 + a10 * b21 + a20 * b22 + a30 * b23);
3283 : //CATCH_REQUIRE(a[2][1] == a01 * b20 + a11 * b21 + a21 * b22 + a31 * b23);
3284 : //CATCH_REQUIRE(a[2][2] == a02 * b20 + a12 * b21 + a22 * b22 + a32 * b23);
3285 : //CATCH_REQUIRE(a[2][3] == a03 * b20 + a13 * b21 + a23 * b22 + a33 * b23);
3286 :
3287 : //CATCH_REQUIRE(a[3][0] == a00 * b30 + a10 * b31 + a20 * b32 + a30 * b33);
3288 : //CATCH_REQUIRE(a[3][1] == a01 * b30 + a11 * b31 + a21 * b32 + a31 * b33);
3289 : //CATCH_REQUIRE(a[3][2] == a02 * b30 + a12 * b31 + a22 * b32 + a32 * b33);
3290 : //CATCH_REQUIRE(a[3][3] == a03 * b30 + a13 * b31 + a23 * b32 + a33 * b33);
3291 1 : }
3292 4 : CATCH_END_SECTION()
3293 8 : }
3294 :
3295 : // create two random 4x4 matrices and make sure the add works
3296 : //
3297 8 : CATCH_GIVEN("div")
3298 : {
3299 4 : CATCH_START_SECTION("matrix: b=a/<scalar>")
3300 : {
3301 : // setup A
3302 : //
3303 1 : snapdev::matrix<double> a(4, 4);
3304 :
3305 1 : CATCH_REQUIRE(a.rows() == 4);
3306 1 : CATCH_REQUIRE(a.columns() == 4);
3307 :
3308 1 : double const a00(frand());
3309 1 : double const a01(frand());
3310 1 : double const a02(frand());
3311 1 : double const a03(frand());
3312 1 : double const a10(frand());
3313 1 : double const a11(frand());
3314 1 : double const a12(frand());
3315 1 : double const a13(frand());
3316 1 : double const a20(frand());
3317 1 : double const a21(frand());
3318 1 : double const a22(frand());
3319 1 : double const a23(frand());
3320 1 : double const a30(frand());
3321 1 : double const a31(frand());
3322 1 : double const a32(frand());
3323 1 : double const a33(frand());
3324 :
3325 1 : a[0][0] = a00;
3326 1 : a[0][1] = a01;
3327 1 : a[0][2] = a02;
3328 1 : a[0][3] = a03;
3329 1 : a[1][0] = a10;
3330 1 : a[1][1] = a11;
3331 1 : a[1][2] = a12;
3332 1 : a[1][3] = a13;
3333 1 : a[2][0] = a20;
3334 1 : a[2][1] = a21;
3335 1 : a[2][2] = a22;
3336 1 : a[2][3] = a23;
3337 1 : a[3][0] = a30;
3338 1 : a[3][1] = a31;
3339 1 : a[3][2] = a32;
3340 1 : a[3][3] = a33;
3341 :
3342 1 : CATCH_REQUIRE(a[0][0] == a00);
3343 1 : CATCH_REQUIRE(a[0][1] == a01);
3344 1 : CATCH_REQUIRE(a[0][2] == a02);
3345 1 : CATCH_REQUIRE(a[0][3] == a03);
3346 1 : CATCH_REQUIRE(a[1][0] == a10);
3347 1 : CATCH_REQUIRE(a[1][1] == a11);
3348 1 : CATCH_REQUIRE(a[1][2] == a12);
3349 1 : CATCH_REQUIRE(a[1][3] == a13);
3350 1 : CATCH_REQUIRE(a[2][0] == a20);
3351 1 : CATCH_REQUIRE(a[2][1] == a21);
3352 1 : CATCH_REQUIRE(a[2][2] == a22);
3353 1 : CATCH_REQUIRE(a[2][3] == a23);
3354 1 : CATCH_REQUIRE(a[3][0] == a30);
3355 1 : CATCH_REQUIRE(a[3][1] == a31);
3356 1 : CATCH_REQUIRE(a[3][2] == a32);
3357 1 : CATCH_REQUIRE(a[3][3] == a33);
3358 :
3359 : // setup B
3360 : //
3361 1 : snapdev::matrix<double> b(4, 4);
3362 :
3363 1 : CATCH_REQUIRE(b.rows() == 4);
3364 1 : CATCH_REQUIRE(b.columns() == 4);
3365 :
3366 1 : double const b00(frand());
3367 1 : double const b01(frand());
3368 1 : double const b02(frand());
3369 1 : double const b03(frand());
3370 1 : double const b10(frand());
3371 1 : double const b11(frand());
3372 1 : double const b12(frand());
3373 1 : double const b13(frand());
3374 1 : double const b20(frand());
3375 1 : double const b21(frand());
3376 1 : double const b22(frand());
3377 1 : double const b23(frand());
3378 1 : double const b30(frand());
3379 1 : double const b31(frand());
3380 1 : double const b32(frand());
3381 1 : double const b33(frand());
3382 :
3383 1 : b[0][0] = b00;
3384 1 : b[0][1] = b01;
3385 1 : b[0][2] = b02;
3386 1 : b[0][3] = b03;
3387 1 : b[1][0] = b10;
3388 1 : b[1][1] = b11;
3389 1 : b[1][2] = b12;
3390 1 : b[1][3] = b13;
3391 1 : b[2][0] = b20;
3392 1 : b[2][1] = b21;
3393 1 : b[2][2] = b22;
3394 1 : b[2][3] = b23;
3395 1 : b[3][0] = b30;
3396 1 : b[3][1] = b31;
3397 1 : b[3][2] = b32;
3398 1 : b[3][3] = b33;
3399 :
3400 1 : CATCH_REQUIRE(b[0][0] == b00);
3401 1 : CATCH_REQUIRE(b[0][1] == b01);
3402 1 : CATCH_REQUIRE(b[0][2] == b02);
3403 1 : CATCH_REQUIRE(b[0][3] == b03);
3404 1 : CATCH_REQUIRE(b[1][0] == b10);
3405 1 : CATCH_REQUIRE(b[1][1] == b11);
3406 1 : CATCH_REQUIRE(b[1][2] == b12);
3407 1 : CATCH_REQUIRE(b[1][3] == b13);
3408 1 : CATCH_REQUIRE(b[2][0] == b20);
3409 1 : CATCH_REQUIRE(b[2][1] == b21);
3410 1 : CATCH_REQUIRE(b[2][2] == b22);
3411 1 : CATCH_REQUIRE(b[2][3] == b23);
3412 1 : CATCH_REQUIRE(b[3][0] == b30);
3413 1 : CATCH_REQUIRE(b[3][1] == b31);
3414 1 : CATCH_REQUIRE(b[3][2] == b32);
3415 1 : CATCH_REQUIRE(b[3][3] == b33);
3416 :
3417 : // create a scalar for our test
3418 : //
3419 1 : double scalar(frand());
3420 1 : while(scalar == 0.0 || scalar == 1.0)
3421 : {
3422 0 : scalar = frand();
3423 : }
3424 :
3425 : // run operation B = A / <scalar>
3426 : //
3427 1 : b = a / scalar;
3428 :
3429 1 : CATCH_REQUIRE(a[0][0] == a00);
3430 1 : CATCH_REQUIRE(a[0][1] == a01);
3431 1 : CATCH_REQUIRE(a[0][2] == a02);
3432 1 : CATCH_REQUIRE(a[0][3] == a03);
3433 1 : CATCH_REQUIRE(a[1][0] == a10);
3434 1 : CATCH_REQUIRE(a[1][1] == a11);
3435 1 : CATCH_REQUIRE(a[1][2] == a12);
3436 1 : CATCH_REQUIRE(a[1][3] == a13);
3437 1 : CATCH_REQUIRE(a[2][0] == a20);
3438 1 : CATCH_REQUIRE(a[2][1] == a21);
3439 1 : CATCH_REQUIRE(a[2][2] == a22);
3440 1 : CATCH_REQUIRE(a[2][3] == a23);
3441 1 : CATCH_REQUIRE(a[3][0] == a30);
3442 1 : CATCH_REQUIRE(a[3][1] == a31);
3443 1 : CATCH_REQUIRE(a[3][2] == a32);
3444 1 : CATCH_REQUIRE(a[3][3] == a33);
3445 :
3446 : // this can't fail because we ensure scalar != 0 or 1
3447 : //
3448 1 : CATCH_REQUIRE_FALSE(b[0][0] == b00);
3449 1 : CATCH_REQUIRE_FALSE(b[0][1] == b01);
3450 1 : CATCH_REQUIRE_FALSE(b[0][2] == b02);
3451 1 : CATCH_REQUIRE_FALSE(b[0][3] == b03);
3452 1 : CATCH_REQUIRE_FALSE(b[1][0] == b10);
3453 1 : CATCH_REQUIRE_FALSE(b[1][1] == b11);
3454 1 : CATCH_REQUIRE_FALSE(b[1][2] == b12);
3455 1 : CATCH_REQUIRE_FALSE(b[1][3] == b13);
3456 1 : CATCH_REQUIRE_FALSE(b[2][0] == b20);
3457 1 : CATCH_REQUIRE_FALSE(b[2][1] == b21);
3458 1 : CATCH_REQUIRE_FALSE(b[2][2] == b22);
3459 1 : CATCH_REQUIRE_FALSE(b[2][3] == b23);
3460 1 : CATCH_REQUIRE_FALSE(b[3][0] == b30);
3461 1 : CATCH_REQUIRE_FALSE(b[3][1] == b31);
3462 1 : CATCH_REQUIRE_FALSE(b[3][2] == b32);
3463 1 : CATCH_REQUIRE_FALSE(b[3][3] == b33);
3464 :
3465 1 : CATCH_REQUIRE(b[0][0] == a00 / scalar);
3466 1 : CATCH_REQUIRE(b[0][1] == a01 / scalar);
3467 1 : CATCH_REQUIRE(b[0][2] == a02 / scalar);
3468 1 : CATCH_REQUIRE(b[0][3] == a03 / scalar);
3469 1 : CATCH_REQUIRE(b[1][0] == a10 / scalar);
3470 1 : CATCH_REQUIRE(b[1][1] == a11 / scalar);
3471 1 : CATCH_REQUIRE(b[1][2] == a12 / scalar);
3472 1 : CATCH_REQUIRE(b[1][3] == a13 / scalar);
3473 1 : CATCH_REQUIRE(b[2][0] == a20 / scalar);
3474 1 : CATCH_REQUIRE(b[2][1] == a21 / scalar);
3475 1 : CATCH_REQUIRE(b[2][2] == a22 / scalar);
3476 1 : CATCH_REQUIRE(b[2][3] == a23 / scalar);
3477 1 : CATCH_REQUIRE(b[3][0] == a30 / scalar);
3478 1 : CATCH_REQUIRE(b[3][1] == a31 / scalar);
3479 1 : CATCH_REQUIRE(b[3][2] == a32 / scalar);
3480 1 : CATCH_REQUIRE(b[3][3] == a33 / scalar);
3481 1 : }
3482 4 : CATCH_END_SECTION()
3483 :
3484 4 : CATCH_START_SECTION("matrix: a/=<scalar>")
3485 : {
3486 : // setup A
3487 : //
3488 1 : snapdev::matrix<double> a(4, 4);
3489 :
3490 1 : CATCH_REQUIRE(a.rows() == 4);
3491 1 : CATCH_REQUIRE(a.columns() == 4);
3492 :
3493 1 : double const a00(frand());
3494 1 : double const a01(frand());
3495 1 : double const a02(frand());
3496 1 : double const a03(frand());
3497 1 : double const a10(frand());
3498 1 : double const a11(frand());
3499 1 : double const a12(frand());
3500 1 : double const a13(frand());
3501 1 : double const a20(frand());
3502 1 : double const a21(frand());
3503 1 : double const a22(frand());
3504 1 : double const a23(frand());
3505 1 : double const a30(frand());
3506 1 : double const a31(frand());
3507 1 : double const a32(frand());
3508 1 : double const a33(frand());
3509 :
3510 1 : a[0][0] = a00;
3511 1 : a[0][1] = a01;
3512 1 : a[0][2] = a02;
3513 1 : a[0][3] = a03;
3514 1 : a[1][0] = a10;
3515 1 : a[1][1] = a11;
3516 1 : a[1][2] = a12;
3517 1 : a[1][3] = a13;
3518 1 : a[2][0] = a20;
3519 1 : a[2][1] = a21;
3520 1 : a[2][2] = a22;
3521 1 : a[2][3] = a23;
3522 1 : a[3][0] = a30;
3523 1 : a[3][1] = a31;
3524 1 : a[3][2] = a32;
3525 1 : a[3][3] = a33;
3526 :
3527 1 : CATCH_REQUIRE(a[0][0] == a00);
3528 1 : CATCH_REQUIRE(a[0][1] == a01);
3529 1 : CATCH_REQUIRE(a[0][2] == a02);
3530 1 : CATCH_REQUIRE(a[0][3] == a03);
3531 1 : CATCH_REQUIRE(a[1][0] == a10);
3532 1 : CATCH_REQUIRE(a[1][1] == a11);
3533 1 : CATCH_REQUIRE(a[1][2] == a12);
3534 1 : CATCH_REQUIRE(a[1][3] == a13);
3535 1 : CATCH_REQUIRE(a[2][0] == a20);
3536 1 : CATCH_REQUIRE(a[2][1] == a21);
3537 1 : CATCH_REQUIRE(a[2][2] == a22);
3538 1 : CATCH_REQUIRE(a[2][3] == a23);
3539 1 : CATCH_REQUIRE(a[3][0] == a30);
3540 1 : CATCH_REQUIRE(a[3][1] == a31);
3541 1 : CATCH_REQUIRE(a[3][2] == a32);
3542 1 : CATCH_REQUIRE(a[3][3] == a33);
3543 :
3544 : // create a scalar for our test
3545 : //
3546 1 : double scalar(frand());
3547 1 : while(scalar == 0.0)
3548 : {
3549 0 : scalar = frand();
3550 : }
3551 :
3552 : // run operation A /= <scalar>
3553 : //
3554 1 : a /= scalar;
3555 :
3556 : // this can't fail because we ensure scalar != 0 or 1
3557 : //
3558 1 : CATCH_REQUIRE_FALSE(a[0][0] == a00);
3559 1 : CATCH_REQUIRE_FALSE(a[0][1] == a01);
3560 1 : CATCH_REQUIRE_FALSE(a[0][2] == a02);
3561 1 : CATCH_REQUIRE_FALSE(a[0][3] == a03);
3562 1 : CATCH_REQUIRE_FALSE(a[1][0] == a10);
3563 1 : CATCH_REQUIRE_FALSE(a[1][1] == a11);
3564 1 : CATCH_REQUIRE_FALSE(a[1][2] == a12);
3565 1 : CATCH_REQUIRE_FALSE(a[1][3] == a13);
3566 1 : CATCH_REQUIRE_FALSE(a[2][0] == a20);
3567 1 : CATCH_REQUIRE_FALSE(a[2][1] == a21);
3568 1 : CATCH_REQUIRE_FALSE(a[2][2] == a22);
3569 1 : CATCH_REQUIRE_FALSE(a[2][3] == a23);
3570 1 : CATCH_REQUIRE_FALSE(a[3][0] == a30);
3571 1 : CATCH_REQUIRE_FALSE(a[3][1] == a31);
3572 1 : CATCH_REQUIRE_FALSE(a[3][2] == a32);
3573 1 : CATCH_REQUIRE_FALSE(a[3][3] == a33);
3574 :
3575 1 : CATCH_REQUIRE(a[0][0] == a00 / scalar);
3576 1 : CATCH_REQUIRE(a[0][1] == a01 / scalar);
3577 1 : CATCH_REQUIRE(a[0][2] == a02 / scalar);
3578 1 : CATCH_REQUIRE(a[0][3] == a03 / scalar);
3579 1 : CATCH_REQUIRE(a[1][0] == a10 / scalar);
3580 1 : CATCH_REQUIRE(a[1][1] == a11 / scalar);
3581 1 : CATCH_REQUIRE(a[1][2] == a12 / scalar);
3582 1 : CATCH_REQUIRE(a[1][3] == a13 / scalar);
3583 1 : CATCH_REQUIRE(a[2][0] == a20 / scalar);
3584 1 : CATCH_REQUIRE(a[2][1] == a21 / scalar);
3585 1 : CATCH_REQUIRE(a[2][2] == a22 / scalar);
3586 1 : CATCH_REQUIRE(a[2][3] == a23 / scalar);
3587 1 : CATCH_REQUIRE(a[3][0] == a30 / scalar);
3588 1 : CATCH_REQUIRE(a[3][1] == a31 / scalar);
3589 1 : CATCH_REQUIRE(a[3][2] == a32 / scalar);
3590 1 : CATCH_REQUIRE(a[3][3] == a33 / scalar);
3591 1 : }
3592 4 : CATCH_END_SECTION()
3593 :
3594 4 : CATCH_START_SECTION("matrix: c=a/b")
3595 : {
3596 : // setup A
3597 : //
3598 1 : snapdev::matrix<double> a(4, 4);
3599 :
3600 1 : CATCH_REQUIRE(a.rows() == 4);
3601 1 : CATCH_REQUIRE(a.columns() == 4);
3602 :
3603 1 : double const a00(frand());
3604 1 : double const a01(frand());
3605 1 : double const a02(frand());
3606 1 : double const a03(frand());
3607 1 : double const a10(frand());
3608 1 : double const a11(frand());
3609 1 : double const a12(frand());
3610 1 : double const a13(frand());
3611 1 : double const a20(frand());
3612 1 : double const a21(frand());
3613 1 : double const a22(frand());
3614 1 : double const a23(frand());
3615 1 : double const a30(frand());
3616 1 : double const a31(frand());
3617 1 : double const a32(frand());
3618 1 : double const a33(frand());
3619 :
3620 1 : a[0][0] = a00;
3621 1 : a[0][1] = a01;
3622 1 : a[0][2] = a02;
3623 1 : a[0][3] = a03;
3624 1 : a[1][0] = a10;
3625 1 : a[1][1] = a11;
3626 1 : a[1][2] = a12;
3627 1 : a[1][3] = a13;
3628 1 : a[2][0] = a20;
3629 1 : a[2][1] = a21;
3630 1 : a[2][2] = a22;
3631 1 : a[2][3] = a23;
3632 1 : a[3][0] = a30;
3633 1 : a[3][1] = a31;
3634 1 : a[3][2] = a32;
3635 1 : a[3][3] = a33;
3636 :
3637 1 : CATCH_REQUIRE(a[0][0] == a00);
3638 1 : CATCH_REQUIRE(a[0][1] == a01);
3639 1 : CATCH_REQUIRE(a[0][2] == a02);
3640 1 : CATCH_REQUIRE(a[0][3] == a03);
3641 1 : CATCH_REQUIRE(a[1][0] == a10);
3642 1 : CATCH_REQUIRE(a[1][1] == a11);
3643 1 : CATCH_REQUIRE(a[1][2] == a12);
3644 1 : CATCH_REQUIRE(a[1][3] == a13);
3645 1 : CATCH_REQUIRE(a[2][0] == a20);
3646 1 : CATCH_REQUIRE(a[2][1] == a21);
3647 1 : CATCH_REQUIRE(a[2][2] == a22);
3648 1 : CATCH_REQUIRE(a[2][3] == a23);
3649 1 : CATCH_REQUIRE(a[3][0] == a30);
3650 1 : CATCH_REQUIRE(a[3][1] == a31);
3651 1 : CATCH_REQUIRE(a[3][2] == a32);
3652 1 : CATCH_REQUIRE(a[3][3] == a33);
3653 :
3654 : // setup B
3655 : //
3656 1 : snapdev::matrix<double> b(4, 4);
3657 :
3658 1 : CATCH_REQUIRE(b.rows() == 4);
3659 1 : CATCH_REQUIRE(b.columns() == 4);
3660 :
3661 : double b00;
3662 : double b01;
3663 : double b02;
3664 : double b03;
3665 : double b10;
3666 : double b11;
3667 : double b12;
3668 : double b13;
3669 : double b20;
3670 : double b21;
3671 : double b22;
3672 : double b23;
3673 : double b30;
3674 : double b31;
3675 : double b32;
3676 : double b33;
3677 :
3678 : // create an inversible matrix (most are with random numbers)
3679 1 : snapdev::matrix<double> t(4, 4);
3680 : do
3681 : {
3682 1 : b00 = frand();
3683 1 : b01 = frand();
3684 1 : b02 = frand();
3685 1 : b03 = frand();
3686 1 : b10 = frand();
3687 1 : b11 = frand();
3688 1 : b12 = frand();
3689 1 : b13 = frand();
3690 1 : b20 = frand();
3691 1 : b21 = frand();
3692 1 : b22 = frand();
3693 1 : b23 = frand();
3694 1 : b30 = frand();
3695 1 : b31 = frand();
3696 1 : b32 = frand();
3697 1 : b33 = frand();
3698 :
3699 1 : b[0][0] = b00;
3700 1 : b[0][1] = b01;
3701 1 : b[0][2] = b02;
3702 1 : b[0][3] = b03;
3703 1 : b[1][0] = b10;
3704 1 : b[1][1] = b11;
3705 1 : b[1][2] = b12;
3706 1 : b[1][3] = b13;
3707 1 : b[2][0] = b20;
3708 1 : b[2][1] = b21;
3709 1 : b[2][2] = b22;
3710 1 : b[2][3] = b23;
3711 1 : b[3][0] = b30;
3712 1 : b[3][1] = b31;
3713 1 : b[3][2] = b32;
3714 1 : b[3][3] = b33;
3715 :
3716 1 : CATCH_REQUIRE(b[0][0] == b00);
3717 1 : CATCH_REQUIRE(b[0][1] == b01);
3718 1 : CATCH_REQUIRE(b[0][2] == b02);
3719 1 : CATCH_REQUIRE(b[0][3] == b03);
3720 1 : CATCH_REQUIRE(b[1][0] == b10);
3721 1 : CATCH_REQUIRE(b[1][1] == b11);
3722 1 : CATCH_REQUIRE(b[1][2] == b12);
3723 1 : CATCH_REQUIRE(b[1][3] == b13);
3724 1 : CATCH_REQUIRE(b[2][0] == b20);
3725 1 : CATCH_REQUIRE(b[2][1] == b21);
3726 1 : CATCH_REQUIRE(b[2][2] == b22);
3727 1 : CATCH_REQUIRE(b[2][3] == b23);
3728 1 : CATCH_REQUIRE(b[3][0] == b30);
3729 1 : CATCH_REQUIRE(b[3][1] == b31);
3730 1 : CATCH_REQUIRE(b[3][2] == b32);
3731 1 : CATCH_REQUIRE(b[3][3] == b33);
3732 :
3733 1 : t = b;
3734 : }
3735 1 : while(!t.inverse());
3736 :
3737 : // setup C
3738 : //
3739 1 : snapdev::matrix<double> c(4, 4);
3740 :
3741 1 : CATCH_REQUIRE(c.rows() == 4);
3742 1 : CATCH_REQUIRE(c.columns() == 4);
3743 :
3744 1 : double const c00(frand());
3745 1 : double const c01(frand());
3746 1 : double const c02(frand());
3747 1 : double const c03(frand());
3748 1 : double const c10(frand());
3749 1 : double const c11(frand());
3750 1 : double const c12(frand());
3751 1 : double const c13(frand());
3752 1 : double const c20(frand());
3753 1 : double const c21(frand());
3754 1 : double const c22(frand());
3755 1 : double const c23(frand());
3756 1 : double const c30(frand());
3757 1 : double const c31(frand());
3758 1 : double const c32(frand());
3759 1 : double const c33(frand());
3760 :
3761 1 : c[0][0] = c00;
3762 1 : c[0][1] = c01;
3763 1 : c[0][2] = c02;
3764 1 : c[0][3] = c03;
3765 1 : c[1][0] = c10;
3766 1 : c[1][1] = c11;
3767 1 : c[1][2] = c12;
3768 1 : c[1][3] = c13;
3769 1 : c[2][0] = c20;
3770 1 : c[2][1] = c21;
3771 1 : c[2][2] = c22;
3772 1 : c[2][3] = c23;
3773 1 : c[3][0] = c30;
3774 1 : c[3][1] = c31;
3775 1 : c[3][2] = c32;
3776 1 : c[3][3] = c33;
3777 :
3778 1 : CATCH_REQUIRE(c[0][0] == c00);
3779 1 : CATCH_REQUIRE(c[0][1] == c01);
3780 1 : CATCH_REQUIRE(c[0][2] == c02);
3781 1 : CATCH_REQUIRE(c[0][3] == c03);
3782 1 : CATCH_REQUIRE(c[1][0] == c10);
3783 1 : CATCH_REQUIRE(c[1][1] == c11);
3784 1 : CATCH_REQUIRE(c[1][2] == c12);
3785 1 : CATCH_REQUIRE(c[1][3] == c13);
3786 1 : CATCH_REQUIRE(c[2][0] == c20);
3787 1 : CATCH_REQUIRE(c[2][1] == c21);
3788 1 : CATCH_REQUIRE(c[2][2] == c22);
3789 1 : CATCH_REQUIRE(c[2][3] == c23);
3790 1 : CATCH_REQUIRE(c[3][0] == c30);
3791 1 : CATCH_REQUIRE(c[3][1] == c31);
3792 1 : CATCH_REQUIRE(c[3][2] == c32);
3793 1 : CATCH_REQUIRE(c[3][3] == c33);
3794 :
3795 : // run operation C = A / B
3796 : //
3797 1 : c = a / b;
3798 :
3799 1 : CATCH_REQUIRE(a[0][0] == a00);
3800 1 : CATCH_REQUIRE(a[0][1] == a01);
3801 1 : CATCH_REQUIRE(a[0][2] == a02);
3802 1 : CATCH_REQUIRE(a[0][3] == a03);
3803 1 : CATCH_REQUIRE(a[1][0] == a10);
3804 1 : CATCH_REQUIRE(a[1][1] == a11);
3805 1 : CATCH_REQUIRE(a[1][2] == a12);
3806 1 : CATCH_REQUIRE(a[1][3] == a13);
3807 1 : CATCH_REQUIRE(a[2][0] == a20);
3808 1 : CATCH_REQUIRE(a[2][1] == a21);
3809 1 : CATCH_REQUIRE(a[2][2] == a22);
3810 1 : CATCH_REQUIRE(a[2][3] == a23);
3811 1 : CATCH_REQUIRE(a[3][0] == a30);
3812 1 : CATCH_REQUIRE(a[3][1] == a31);
3813 1 : CATCH_REQUIRE(a[3][2] == a32);
3814 1 : CATCH_REQUIRE(a[3][3] == a33);
3815 :
3816 1 : CATCH_REQUIRE(b[0][0] == b00);
3817 1 : CATCH_REQUIRE(b[0][1] == b01);
3818 1 : CATCH_REQUIRE(b[0][2] == b02);
3819 1 : CATCH_REQUIRE(b[0][3] == b03);
3820 1 : CATCH_REQUIRE(b[1][0] == b10);
3821 1 : CATCH_REQUIRE(b[1][1] == b11);
3822 1 : CATCH_REQUIRE(b[1][2] == b12);
3823 1 : CATCH_REQUIRE(b[1][3] == b13);
3824 1 : CATCH_REQUIRE(b[2][0] == b20);
3825 1 : CATCH_REQUIRE(b[2][1] == b21);
3826 1 : CATCH_REQUIRE(b[2][2] == b22);
3827 1 : CATCH_REQUIRE(b[2][3] == b23);
3828 1 : CATCH_REQUIRE(b[3][0] == b30);
3829 1 : CATCH_REQUIRE(b[3][1] == b31);
3830 1 : CATCH_REQUIRE(b[3][2] == b32);
3831 1 : CATCH_REQUIRE(b[3][3] == b33);
3832 :
3833 : // this could fail
3834 : //CATCH_REQUIRE_FALSE(c[0][0] == c00);
3835 : //CATCH_REQUIRE_FALSE(c[0][1] == c01);
3836 : //CATCH_REQUIRE_FALSE(c[0][2] == c02);
3837 : //CATCH_REQUIRE_FALSE(c[0][3] == c03);
3838 : //CATCH_REQUIRE_FALSE(c[1][0] == c10);
3839 : //CATCH_REQUIRE_FALSE(c[1][1] == c11);
3840 : //CATCH_REQUIRE_FALSE(c[1][2] == c12);
3841 : //CATCH_REQUIRE_FALSE(c[1][3] == c13);
3842 : //CATCH_REQUIRE_FALSE(c[2][0] == c20);
3843 : //CATCH_REQUIRE_FALSE(c[2][1] == c21);
3844 : //CATCH_REQUIRE_FALSE(c[2][2] == c22);
3845 : //CATCH_REQUIRE_FALSE(c[2][3] == c23);
3846 : //CATCH_REQUIRE_FALSE(c[3][0] == c30);
3847 : //CATCH_REQUIRE_FALSE(c[3][1] == c31);
3848 : //CATCH_REQUIRE_FALSE(c[3][2] == c32);
3849 : //CATCH_REQUIRE_FALSE(c[3][3] == c33);
3850 :
3851 : // first level verification, which is exactly what the
3852 : // division function does, it does not mean that the
3853 : // division works properly, though
3854 : //
3855 1 : snapdev::matrix<double> r(4, 4);
3856 1 : r = a * t;
3857 :
3858 1 : CATCH_REQUIRE(c[0][0] == r[0][0]);
3859 1 : CATCH_REQUIRE(c[0][1] == r[0][1]);
3860 1 : CATCH_REQUIRE(c[0][2] == r[0][2]);
3861 1 : CATCH_REQUIRE(c[0][3] == r[0][3]);
3862 1 : CATCH_REQUIRE(c[1][0] == r[1][0]);
3863 1 : CATCH_REQUIRE(c[1][1] == r[1][1]);
3864 1 : CATCH_REQUIRE(c[1][2] == r[1][2]);
3865 1 : CATCH_REQUIRE(c[1][3] == r[1][3]);
3866 1 : CATCH_REQUIRE(c[2][0] == r[2][0]);
3867 1 : CATCH_REQUIRE(c[2][1] == r[2][1]);
3868 1 : CATCH_REQUIRE(c[2][2] == r[2][2]);
3869 1 : CATCH_REQUIRE(c[2][3] == r[2][3]);
3870 1 : CATCH_REQUIRE(c[3][0] == r[3][0]);
3871 1 : CATCH_REQUIRE(c[3][1] == r[3][1]);
3872 1 : CATCH_REQUIRE(c[3][2] == r[3][2]);
3873 1 : CATCH_REQUIRE(c[3][3] == r[3][3]);
3874 :
3875 1 : double const determinant(b.determinant());
3876 1 : snapdev::matrix<double> adjugate(b.adjugate());
3877 :
3878 1 : snapdev::matrix<double> inv(adjugate * (1.0 / determinant));
3879 1 : snapdev::matrix<double> div(a * inv);
3880 :
3881 1 : CATCH_REQUIRE(fabs(c[0][0] - div[0][0]) < 0.0001);
3882 1 : CATCH_REQUIRE(fabs(c[0][1] - div[0][1]) < 0.0001);
3883 1 : CATCH_REQUIRE(fabs(c[0][2] - div[0][2]) < 0.0001);
3884 1 : CATCH_REQUIRE(fabs(c[0][3] - div[0][3]) < 0.0001);
3885 1 : CATCH_REQUIRE(fabs(c[1][0] - div[1][0]) < 0.0001);
3886 1 : CATCH_REQUIRE(fabs(c[1][1] - div[1][1]) < 0.0001);
3887 1 : CATCH_REQUIRE(fabs(c[1][2] - div[1][2]) < 0.0001);
3888 1 : CATCH_REQUIRE(fabs(c[1][3] - div[1][3]) < 0.0001);
3889 1 : CATCH_REQUIRE(fabs(c[2][0] - div[2][0]) < 0.0001);
3890 1 : CATCH_REQUIRE(fabs(c[2][1] - div[2][1]) < 0.0001);
3891 1 : CATCH_REQUIRE(fabs(c[2][2] - div[2][2]) < 0.0001);
3892 1 : CATCH_REQUIRE(fabs(c[2][3] - div[2][3]) < 0.0001);
3893 1 : CATCH_REQUIRE(fabs(c[3][0] - div[3][0]) < 0.0001);
3894 1 : CATCH_REQUIRE(fabs(c[3][1] - div[3][1]) < 0.0001);
3895 1 : CATCH_REQUIRE(fabs(c[3][2] - div[3][2]) < 0.0001);
3896 1 : CATCH_REQUIRE(fabs(c[3][3] - div[3][3]) < 0.0001);
3897 1 : }
3898 4 : CATCH_END_SECTION()
3899 :
3900 4 : CATCH_START_SECTION("matrix: a/=b")
3901 : {
3902 : // setup A
3903 : //
3904 1 : snapdev::matrix<double> a(4, 4);
3905 :
3906 1 : CATCH_REQUIRE(a.rows() == 4);
3907 1 : CATCH_REQUIRE(a.columns() == 4);
3908 :
3909 1 : double const a00(frand());
3910 1 : double const a01(frand());
3911 1 : double const a02(frand());
3912 1 : double const a03(frand());
3913 1 : double const a10(frand());
3914 1 : double const a11(frand());
3915 1 : double const a12(frand());
3916 1 : double const a13(frand());
3917 1 : double const a20(frand());
3918 1 : double const a21(frand());
3919 1 : double const a22(frand());
3920 1 : double const a23(frand());
3921 1 : double const a30(frand());
3922 1 : double const a31(frand());
3923 1 : double const a32(frand());
3924 1 : double const a33(frand());
3925 :
3926 1 : a[0][0] = a00;
3927 1 : a[0][1] = a01;
3928 1 : a[0][2] = a02;
3929 1 : a[0][3] = a03;
3930 1 : a[1][0] = a10;
3931 1 : a[1][1] = a11;
3932 1 : a[1][2] = a12;
3933 1 : a[1][3] = a13;
3934 1 : a[2][0] = a20;
3935 1 : a[2][1] = a21;
3936 1 : a[2][2] = a22;
3937 1 : a[2][3] = a23;
3938 1 : a[3][0] = a30;
3939 1 : a[3][1] = a31;
3940 1 : a[3][2] = a32;
3941 1 : a[3][3] = a33;
3942 :
3943 1 : CATCH_REQUIRE(a[0][0] == a00);
3944 1 : CATCH_REQUIRE(a[0][1] == a01);
3945 1 : CATCH_REQUIRE(a[0][2] == a02);
3946 1 : CATCH_REQUIRE(a[0][3] == a03);
3947 1 : CATCH_REQUIRE(a[1][0] == a10);
3948 1 : CATCH_REQUIRE(a[1][1] == a11);
3949 1 : CATCH_REQUIRE(a[1][2] == a12);
3950 1 : CATCH_REQUIRE(a[1][3] == a13);
3951 1 : CATCH_REQUIRE(a[2][0] == a20);
3952 1 : CATCH_REQUIRE(a[2][1] == a21);
3953 1 : CATCH_REQUIRE(a[2][2] == a22);
3954 1 : CATCH_REQUIRE(a[2][3] == a23);
3955 1 : CATCH_REQUIRE(a[3][0] == a30);
3956 1 : CATCH_REQUIRE(a[3][1] == a31);
3957 1 : CATCH_REQUIRE(a[3][2] == a32);
3958 1 : CATCH_REQUIRE(a[3][3] == a33);
3959 :
3960 : // setup B
3961 : //
3962 1 : snapdev::matrix<double> b(4, 4);
3963 :
3964 1 : CATCH_REQUIRE(b.rows() == 4);
3965 1 : CATCH_REQUIRE(b.columns() == 4);
3966 :
3967 1 : double const b00(frand());
3968 1 : double const b01(frand());
3969 1 : double const b02(frand());
3970 1 : double const b03(frand());
3971 1 : double const b10(frand());
3972 1 : double const b11(frand());
3973 1 : double const b12(frand());
3974 1 : double const b13(frand());
3975 1 : double const b20(frand());
3976 1 : double const b21(frand());
3977 1 : double const b22(frand());
3978 1 : double const b23(frand());
3979 1 : double const b30(frand());
3980 1 : double const b31(frand());
3981 1 : double const b32(frand());
3982 1 : double const b33(frand());
3983 :
3984 1 : b[0][0] = b00;
3985 1 : b[0][1] = b01;
3986 1 : b[0][2] = b02;
3987 1 : b[0][3] = b03;
3988 1 : b[1][0] = b10;
3989 1 : b[1][1] = b11;
3990 1 : b[1][2] = b12;
3991 1 : b[1][3] = b13;
3992 1 : b[2][0] = b20;
3993 1 : b[2][1] = b21;
3994 1 : b[2][2] = b22;
3995 1 : b[2][3] = b23;
3996 1 : b[3][0] = b30;
3997 1 : b[3][1] = b31;
3998 1 : b[3][2] = b32;
3999 1 : b[3][3] = b33;
4000 :
4001 1 : CATCH_REQUIRE(b[0][0] == b00);
4002 1 : CATCH_REQUIRE(b[0][1] == b01);
4003 1 : CATCH_REQUIRE(b[0][2] == b02);
4004 1 : CATCH_REQUIRE(b[0][3] == b03);
4005 1 : CATCH_REQUIRE(b[1][0] == b10);
4006 1 : CATCH_REQUIRE(b[1][1] == b11);
4007 1 : CATCH_REQUIRE(b[1][2] == b12);
4008 1 : CATCH_REQUIRE(b[1][3] == b13);
4009 1 : CATCH_REQUIRE(b[2][0] == b20);
4010 1 : CATCH_REQUIRE(b[2][1] == b21);
4011 1 : CATCH_REQUIRE(b[2][2] == b22);
4012 1 : CATCH_REQUIRE(b[2][3] == b23);
4013 1 : CATCH_REQUIRE(b[3][0] == b30);
4014 1 : CATCH_REQUIRE(b[3][1] == b31);
4015 1 : CATCH_REQUIRE(b[3][2] == b32);
4016 1 : CATCH_REQUIRE(b[3][3] == b33);
4017 :
4018 1 : double const determinant(b.determinant());
4019 1 : snapdev::matrix<double> adjugate(b.adjugate());
4020 :
4021 1 : snapdev::matrix<double> inv(adjugate * (1.0 / determinant));
4022 1 : snapdev::matrix<double> div(a * inv);
4023 :
4024 : // run operation A /= B
4025 : //
4026 1 : a /= b;
4027 :
4028 : // this could fail if one of bXX is 0.0
4029 : //
4030 : //CATCH_REQUIRE(a[0][0] == a00);
4031 : //CATCH_REQUIRE(a[0][1] == a01);
4032 : //CATCH_REQUIRE(a[0][2] == a02);
4033 : //CATCH_REQUIRE(a[0][3] == a03);
4034 : //CATCH_REQUIRE(a[1][0] == a10);
4035 : //CATCH_REQUIRE(a[1][1] == a11);
4036 : //CATCH_REQUIRE(a[1][2] == a12);
4037 : //CATCH_REQUIRE(a[1][3] == a13);
4038 : //CATCH_REQUIRE(a[2][0] == a20);
4039 : //CATCH_REQUIRE(a[2][1] == a21);
4040 : //CATCH_REQUIRE(a[2][2] == a22);
4041 : //CATCH_REQUIRE(a[2][3] == a23);
4042 : //CATCH_REQUIRE(a[3][0] == a30);
4043 : //CATCH_REQUIRE(a[3][1] == a31);
4044 : //CATCH_REQUIRE(a[3][2] == a32);
4045 : //CATCH_REQUIRE(a[3][3] == a33);
4046 :
4047 1 : CATCH_REQUIRE(b[0][0] == b00);
4048 1 : CATCH_REQUIRE(b[0][1] == b01);
4049 1 : CATCH_REQUIRE(b[0][2] == b02);
4050 1 : CATCH_REQUIRE(b[0][3] == b03);
4051 1 : CATCH_REQUIRE(b[1][0] == b10);
4052 1 : CATCH_REQUIRE(b[1][1] == b11);
4053 1 : CATCH_REQUIRE(b[1][2] == b12);
4054 1 : CATCH_REQUIRE(b[1][3] == b13);
4055 1 : CATCH_REQUIRE(b[2][0] == b20);
4056 1 : CATCH_REQUIRE(b[2][1] == b21);
4057 1 : CATCH_REQUIRE(b[2][2] == b22);
4058 1 : CATCH_REQUIRE(b[2][3] == b23);
4059 1 : CATCH_REQUIRE(b[3][0] == b30);
4060 1 : CATCH_REQUIRE(b[3][1] == b31);
4061 1 : CATCH_REQUIRE(b[3][2] == b32);
4062 1 : CATCH_REQUIRE(b[3][3] == b33);
4063 :
4064 1 : CATCH_REQUIRE(fabs(a[0][0] - div[0][0]) < 0.0001);
4065 1 : CATCH_REQUIRE(fabs(a[0][1] - div[0][1]) < 0.0001);
4066 1 : CATCH_REQUIRE(fabs(a[0][2] - div[0][2]) < 0.0001);
4067 1 : CATCH_REQUIRE(fabs(a[0][3] - div[0][3]) < 0.0001);
4068 1 : CATCH_REQUIRE(fabs(a[1][0] - div[1][0]) < 0.0001);
4069 1 : CATCH_REQUIRE(fabs(a[1][1] - div[1][1]) < 0.0001);
4070 1 : CATCH_REQUIRE(fabs(a[1][2] - div[1][2]) < 0.0001);
4071 1 : CATCH_REQUIRE(fabs(a[1][3] - div[1][3]) < 0.0001);
4072 1 : CATCH_REQUIRE(fabs(a[2][0] - div[2][0]) < 0.0001);
4073 1 : CATCH_REQUIRE(fabs(a[2][1] - div[2][1]) < 0.0001);
4074 1 : CATCH_REQUIRE(fabs(a[2][2] - div[2][2]) < 0.0001);
4075 1 : CATCH_REQUIRE(fabs(a[2][3] - div[2][3]) < 0.0001);
4076 1 : CATCH_REQUIRE(fabs(a[3][0] - div[3][0]) < 0.0001);
4077 1 : CATCH_REQUIRE(fabs(a[3][1] - div[3][1]) < 0.0001);
4078 1 : CATCH_REQUIRE(fabs(a[3][2] - div[3][2]) < 0.0001);
4079 1 : CATCH_REQUIRE(fabs(a[3][3] - div[3][3]) < 0.0001);
4080 1 : }
4081 4 : CATCH_END_SECTION()
4082 8 : }
4083 8 : }
4084 :
4085 :
4086 3 : CATCH_TEST_CASE("matrix_color", "[matrix][math]")
4087 : {
4088 : // generate a few brightness matrices
4089 : //
4090 3 : CATCH_GIVEN("brightness")
4091 : {
4092 1 : CATCH_START_SECTION("matrix: b=a*<brightness> (a is identity)")
4093 : {
4094 : // setup A
4095 : //
4096 1 : snapdev::matrix<double> a(4, 4);
4097 :
4098 1 : CATCH_REQUIRE(a.rows() == 4);
4099 1 : CATCH_REQUIRE(a.columns() == 4);
4100 :
4101 : // setup B
4102 : //
4103 1 : snapdev::matrix<double> b(4, 4);
4104 :
4105 1 : CATCH_REQUIRE(b.rows() == 4);
4106 1 : CATCH_REQUIRE(b.columns() == 4);
4107 :
4108 1001 : for(int idx(0); idx < 1000; ++idx)
4109 : {
4110 1000 : double const brightness(static_cast<double>(idx) / 1000.0);
4111 1000 : b = a.brightness(brightness);
4112 :
4113 1000 : CATCH_REQUIRE(fabs(b[0][0] - brightness) < 0.0001);
4114 1000 : CATCH_REQUIRE(fabs(b[0][1] - 0.0) < 0.0001);
4115 1000 : CATCH_REQUIRE(fabs(b[0][2] - 0.0) < 0.0001);
4116 1000 : CATCH_REQUIRE(fabs(b[0][3] - 0.0) < 0.0001);
4117 1000 : CATCH_REQUIRE(fabs(b[1][0] - 0.0) < 0.0001);
4118 1000 : CATCH_REQUIRE(fabs(b[1][1] - brightness) < 0.0001);
4119 1000 : CATCH_REQUIRE(fabs(b[1][2] - 0.0) < 0.0001);
4120 1000 : CATCH_REQUIRE(fabs(b[1][3] - 0.0) < 0.0001);
4121 1000 : CATCH_REQUIRE(fabs(b[2][0] - 0.0) < 0.0001);
4122 1000 : CATCH_REQUIRE(fabs(b[2][1] - 0.0) < 0.0001);
4123 1000 : CATCH_REQUIRE(fabs(b[2][2] - brightness) < 0.0001);
4124 1000 : CATCH_REQUIRE(fabs(b[2][3] - 0.0) < 0.0001);
4125 1000 : CATCH_REQUIRE(fabs(b[3][0] - 0.0) < 0.0001);
4126 1000 : CATCH_REQUIRE(fabs(b[3][1] - 0.0) < 0.0001);
4127 1000 : CATCH_REQUIRE(fabs(b[3][2] - 0.0) < 0.0001);
4128 1000 : CATCH_REQUIRE(fabs(b[3][3] - 1.0) < 0.0001);
4129 :
4130 : // * brightness_r & 0 & 0 & 0
4131 : // * \\ 0 & brightness_g & 0 & 0
4132 : // * \\ 0 & 0 & brightness_b & 0
4133 : // * \\ 0 & 0 & 0 & 1
4134 :
4135 : }
4136 1 : }
4137 1 : CATCH_END_SECTION()
4138 3 : }
4139 :
4140 : // generate a few saturation matrices
4141 : //
4142 3 : CATCH_GIVEN("saturation")
4143 : {
4144 1 : CATCH_START_SECTION("matrix: b=a*<saturation> (a is identity)")
4145 : {
4146 : // setup A
4147 : //
4148 1 : snapdev::matrix<double> a(4, 4);
4149 :
4150 1 : CATCH_REQUIRE(a.rows() == 4);
4151 1 : CATCH_REQUIRE(a.columns() == 4);
4152 :
4153 : // setup B
4154 : //
4155 1 : snapdev::matrix<double> b(4, 4);
4156 :
4157 1 : CATCH_REQUIRE(b.rows() == 4);
4158 1 : CATCH_REQUIRE(b.columns() == 4);
4159 :
4160 1 : double luma[6][3];
4161 :
4162 1 : luma[0][0] = a.HDTV_LUMA_RED ;
4163 1 : luma[0][1] = a.HDTV_LUMA_GREEN;
4164 1 : luma[0][2] = a.HDTV_LUMA_BLUE ;
4165 :
4166 1 : luma[1][0] = a.LED_LUMA_RED ;
4167 1 : luma[1][1] = a.LED_LUMA_GREEN;
4168 1 : luma[1][2] = a.LED_LUMA_BLUE ;
4169 :
4170 1 : luma[2][0] = a.CRT_LUMA_RED ;
4171 1 : luma[2][1] = a.CRT_LUMA_GREEN;
4172 1 : luma[2][2] = a.CRT_LUMA_BLUE ;
4173 :
4174 1 : luma[3][0] = a.NTSC_LUMA_RED ;
4175 1 : luma[3][1] = a.NTSC_LUMA_GREEN;
4176 1 : luma[3][2] = a.NTSC_LUMA_BLUE ;
4177 :
4178 1 : luma[4][0] = a.AVERAGE_LUMA_RED ;
4179 1 : luma[4][1] = a.AVERAGE_LUMA_GREEN;
4180 1 : luma[4][2] = a.AVERAGE_LUMA_BLUE ;
4181 :
4182 1 : luma[5][0] = 0.2; // to test any dynamic version if such exists for this algorithm
4183 1 : luma[5][1] = 0.7;
4184 1 : luma[5][2] = 0.1;
4185 :
4186 7 : for(int l(0); l < 6; ++l)
4187 : {
4188 : // set the luma weights
4189 : //
4190 6 : a.set_luma_vector(luma[l][0], luma[l][1], luma[l][2]);
4191 :
4192 6006 : for(int idx(0); idx < 1000; ++idx)
4193 : {
4194 6000 : double const s(static_cast<double>(idx) / 1000.0);
4195 6000 : b = a.saturation(s);
4196 :
4197 6000 : CATCH_REQUIRE(fabs(b[0][0] - (luma[l][0] * (1.0 - s) + s)) < 0.0001);
4198 6000 : CATCH_REQUIRE(fabs(b[0][1] - (luma[l][0] * (1.0 - s))) < 0.0001);
4199 6000 : CATCH_REQUIRE(fabs(b[0][2] - (luma[l][0] * (1.0 - s))) < 0.0001);
4200 6000 : CATCH_REQUIRE(fabs(b[0][3] - (0.0)) < 0.0001);
4201 6000 : CATCH_REQUIRE(fabs(b[1][0] - (luma[l][1] * (1.0 - s))) < 0.0001);
4202 6000 : CATCH_REQUIRE(fabs(b[1][1] - (luma[l][1] * (1.0 - s) + s)) < 0.0001);
4203 6000 : CATCH_REQUIRE(fabs(b[1][2] - (luma[l][1] * (1.0 - s))) < 0.0001);
4204 6000 : CATCH_REQUIRE(fabs(b[1][3] - (0.0)) < 0.0001);
4205 6000 : CATCH_REQUIRE(fabs(b[2][0] - (luma[l][2] * (1.0 - s))) < 0.0001);
4206 6000 : CATCH_REQUIRE(fabs(b[2][1] - (luma[l][2] * (1.0 - s))) < 0.0001);
4207 6000 : CATCH_REQUIRE(fabs(b[2][2] - (luma[l][2] * (1.0 - s) + s)) < 0.0001);
4208 6000 : CATCH_REQUIRE(fabs(b[2][3] - (0.0)) < 0.0001);
4209 6000 : CATCH_REQUIRE(fabs(b[3][0] - (0.0)) < 0.0001);
4210 6000 : CATCH_REQUIRE(fabs(b[3][1] - (0.0)) < 0.0001);
4211 6000 : CATCH_REQUIRE(fabs(b[3][2] - (0.0)) < 0.0001);
4212 6000 : CATCH_REQUIRE(fabs(b[3][3] - (1.0)) < 0.0001);
4213 : }
4214 : }
4215 1 : }
4216 1 : CATCH_END_SECTION()
4217 3 : }
4218 :
4219 : // generate a few brightness matrices
4220 : //
4221 3 : CATCH_GIVEN("hue")
4222 : {
4223 1 : CATCH_START_SECTION("matrix: b=a*<hue> (a is identity)")
4224 : {
4225 : // setup A
4226 : //
4227 1 : snapdev::matrix<double> a(4, 4);
4228 :
4229 1 : CATCH_REQUIRE(a.rows() == 4);
4230 1 : CATCH_REQUIRE(a.columns() == 4);
4231 :
4232 : // setup B
4233 : //
4234 1 : snapdev::matrix<double> b(4, 4);
4235 :
4236 1 : CATCH_REQUIRE(b.rows() == 4);
4237 1 : CATCH_REQUIRE(b.columns() == 4);
4238 :
4239 1 : double luma[6][3];
4240 :
4241 1 : luma[0][0] = a.HDTV_LUMA_RED ;
4242 1 : luma[0][1] = a.HDTV_LUMA_GREEN;
4243 1 : luma[0][2] = a.HDTV_LUMA_BLUE ;
4244 :
4245 1 : luma[1][0] = a.LED_LUMA_RED ;
4246 1 : luma[1][1] = a.LED_LUMA_GREEN;
4247 1 : luma[1][2] = a.LED_LUMA_BLUE ;
4248 :
4249 1 : luma[2][0] = a.CRT_LUMA_RED ;
4250 1 : luma[2][1] = a.CRT_LUMA_GREEN;
4251 1 : luma[2][2] = a.CRT_LUMA_BLUE ;
4252 :
4253 1 : luma[3][0] = a.NTSC_LUMA_RED ;
4254 1 : luma[3][1] = a.NTSC_LUMA_GREEN;
4255 1 : luma[3][2] = a.NTSC_LUMA_BLUE ;
4256 :
4257 1 : luma[4][0] = a.AVERAGE_LUMA_RED ;
4258 1 : luma[4][1] = a.AVERAGE_LUMA_GREEN;
4259 1 : luma[4][2] = a.AVERAGE_LUMA_BLUE ;
4260 :
4261 1 : luma[5][0] = 0.2; // to test any dynamic version if such exists for this algorithm
4262 1 : luma[5][1] = 0.7;
4263 1 : luma[5][2] = 0.1;
4264 :
4265 1 : std::vector<std::string> last_hue_matrix{
4266 : "hdtv",
4267 : "led",
4268 : "crt",
4269 : "ntsc",
4270 : "average",
4271 : "dynamic"
4272 10 : };
4273 :
4274 7 : for(int l(0); l < 6; ++l)
4275 : {
4276 : // set the luma weights
4277 : //
4278 6 : a.set_luma_vector(luma[l][0], luma[l][1], luma[l][2]);
4279 :
4280 : //std::cerr << "*** working on luma " << l << "\n";
4281 6006 : for(int idx(0); idx < 1000; ++idx)
4282 : {
4283 : // hue is an angle from 0 to 360 (albeit in radian)
4284 : //
4285 6000 : double const hue(static_cast<double>(idx) / 1000.0 * M_PI * 2.0);
4286 6000 : b = a.hue(hue);
4287 :
4288 : #ifdef _DEBUG
4289 6000 : CATCH_REQUIRE(a.get_last_hue_matrix() == last_hue_matrix[l]);
4290 : #endif
4291 :
4292 : // this one requires us to work!
4293 : // here we calculate the hue matrix by hand instead of using
4294 : // pre-calculated numbers; that way we know we have it right
4295 : //
4296 :
4297 : // $R_r$ (red rotation)
4298 : //
4299 6000 : snapdev::matrix<double> r_r(4, 4);
4300 6000 : r_r[1][1] = 1.0 / sqrt(2.0);
4301 6000 : r_r[1][2] = 1.0 / sqrt(2.0);
4302 6000 : r_r[2][1] = -1.0 / sqrt(2.0);
4303 6000 : r_r[2][2] = 1.0 / sqrt(2.0);
4304 :
4305 : // $R_g$ (green rotation)
4306 : //
4307 6000 : snapdev::matrix<double> r_g(4, 4);
4308 6000 : r_g[0][0] = sqrt(2.0) / sqrt(3.0);
4309 6000 : r_g[0][2] = 1.0 / sqrt(3.0);
4310 6000 : r_g[2][0] = -1.0 / sqrt(3.0);
4311 6000 : r_g[2][2] = sqrt(2.0) / sqrt(3.0);
4312 :
4313 : // $R_rg = R_r R_g$ (red and green rotations combined)
4314 : //
4315 6000 : snapdev::matrix<double> r_rg(r_r * r_g);
4316 :
4317 : // $w$ (luma vector, a.k.a. color weights to calcalute the perfect grayscale for your monitor)
4318 : //
4319 6000 : snapdev::matrix<double> w(a.get_luma_vector());
4320 : //w[0][0] = a.RED_WEIGHT;
4321 : //w[1][0] = a.GREEN_WEIGHT;
4322 : //w[2][0] = a.BLUE_WEIGHT;
4323 : //w[3][0] = 0;
4324 :
4325 : // sm (skew matrix)
4326 : //
4327 6000 : snapdev::matrix<double> sm(r_rg * w);
4328 :
4329 : // s (scaling with skew matrix)
4330 : //
4331 6000 : snapdev::matrix<double> s(4, 4);
4332 6000 : s[0][2] = sm[0][0] / sm[2][0];
4333 6000 : s[1][2] = sm[1][0] / sm[2][0];
4334 :
4335 : // p (skewed red & green rotation)
4336 : //
4337 6000 : snapdev::matrix<double> p(r_rg);
4338 6000 : p *= s;
4339 :
4340 : // r_b (blue rotation)
4341 : //
4342 6000 : snapdev::matrix<double> r_b(4, 4);
4343 6000 : double const rot_cos(cos(hue));
4344 6000 : double const rot_sin(sin(hue));
4345 6000 : r_b[0][0] = rot_cos;
4346 6000 : r_b[0][1] = rot_sin;
4347 6000 : r_b[1][0] = -rot_sin;
4348 6000 : r_b[1][1] = rot_cos;
4349 :
4350 : // hue_matrix (the resulting matrix depending on angle)
4351 : //
4352 6000 : snapdev::matrix<double> hue_matrix(p * r_b / p);
4353 :
4354 : //std::cout << "Got b = " << b << "\n";
4355 : //std::cout << "Got our hue matrix = " << hue_matrix << "\n";
4356 :
4357 6000 : CATCH_REQUIRE(fabs(b[0][0] - hue_matrix[0][0]) < 0.0001);
4358 6000 : CATCH_REQUIRE(fabs(b[0][1] - hue_matrix[0][1]) < 0.0001);
4359 6000 : CATCH_REQUIRE(fabs(b[0][2] - hue_matrix[0][2]) < 0.0001);
4360 6000 : CATCH_REQUIRE(fabs(b[0][3] - hue_matrix[0][3]) < 0.0001);
4361 6000 : CATCH_REQUIRE(fabs(b[1][0] - hue_matrix[1][0]) < 0.0001);
4362 6000 : CATCH_REQUIRE(fabs(b[1][1] - hue_matrix[1][1]) < 0.0001);
4363 6000 : CATCH_REQUIRE(fabs(b[1][2] - hue_matrix[1][2]) < 0.0001);
4364 6000 : CATCH_REQUIRE(fabs(b[1][3] - hue_matrix[1][3]) < 0.0001);
4365 6000 : CATCH_REQUIRE(fabs(b[2][0] - hue_matrix[2][0]) < 0.0001);
4366 6000 : CATCH_REQUIRE(fabs(b[2][1] - hue_matrix[2][1]) < 0.0001);
4367 6000 : CATCH_REQUIRE(fabs(b[2][2] - hue_matrix[2][2]) < 0.0001);
4368 6000 : CATCH_REQUIRE(fabs(b[2][3] - hue_matrix[2][3]) < 0.0001);
4369 6000 : CATCH_REQUIRE(fabs(b[3][0] - hue_matrix[3][0]) < 0.0001);
4370 6000 : CATCH_REQUIRE(fabs(b[3][1] - hue_matrix[3][1]) < 0.0001);
4371 6000 : CATCH_REQUIRE(fabs(b[3][2] - hue_matrix[3][2]) < 0.0001);
4372 6000 : CATCH_REQUIRE(fabs(b[3][3] - hue_matrix[3][3]) < 0.0001);
4373 :
4374 : // * $$
4375 : // * R_r =
4376 : // * \begin{bmatrix}
4377 : // * 1 & 0 & 0 & 0
4378 : // * \\ 0 & {1 \over \sqrt 2 } & {1 \over \sqrt 2 } & 0
4379 : // * \\ 0 & -{1 \over \sqrt 2 } & {1 \over \sqrt 2 } & 0
4380 : // * \\ 0 & 0 & 0 & 1
4381 : // * \end{bmatrix}
4382 : // * $$
4383 : // *
4384 : // * $$
4385 : // * R_g =
4386 : // * \begin{bmatrix}
4387 : // * 1 & 0 & 0 & 0
4388 : // * \\ 0 & {\sqrt 2 \over \sqrt 3} & -{1 \over \sqrt 3} & 0
4389 : // * \\ 0 & {1 \over \sqrt 3} & {\sqrt 2 \over \sqrt 3} & 0
4390 : // * \\ 0 & 0 & 0 & 1
4391 : // * \end{bmatrix}
4392 : // * $$
4393 : // *
4394 : // * $$
4395 : // * R_{rg} = R_r R_g
4396 : // * $$
4397 : // *
4398 : // * $$
4399 : // * \begin{bmatrix}
4400 : // * L_r
4401 : // * \\ L_g
4402 : // * \\ L_b
4403 : // * \\ 0
4404 : // * \end{bmatrix}
4405 : // * =
4406 : // * R_{rg}
4407 : // * \begin{bmatrix}
4408 : // * W_r
4409 : // * \\ W_g
4410 : // * \\ W_b
4411 : // * \\ 0
4412 : // * \end{bmatrix}
4413 : // * $$
4414 : // *
4415 : // * Now we define a rotation matrix for the blue axis. This one uses
4416 : // * a variable as the angle. This is the actual rotation angle which
4417 : // * can go from 0 to $2 \pi$:
4418 : // *
4419 : // * $$
4420 : // * R_b =
4421 : // * \begin{bmatrix}
4422 : // * cos( \theta ) & sin( \theta ) & 0 & 0
4423 : // * \\ -sin( \theta ) & cos( \theta ) & 0 & 0
4424 : // * \\ 0 & 0 & 1 & 0
4425 : // * \\ 0 & 0 & 0 & 1
4426 : // * \end{bmatrix}
4427 : // * $$
4428 : // *
4429 : // * Now we have all the matrices to caculate the hue rotation
4430 : // * of all the components of an image:
4431 : // *
4432 : // * $$
4433 : // * H = R_{rg} S R_b S^{-1} R_{rg}^{-1}
4434 : // * $$
4435 : // *
4436 : // * $H$ can then be used as in:
4437 : // *
4438 : // * $$
4439 : // * \begin{bmatrix}
4440 : // * R'
4441 : // * \\ G'
4442 : // * \\ B'
4443 : // * \end{bmatrix}
4444 : // * =
4445 : // * H
4446 : // * \begin{bmatrix}
4447 : // * R
4448 : // * \\ G
4449 : // * \\ B
4450 : // * \end{bmatrix}
4451 : // * $$
4452 : // *
4453 6000 : }
4454 : }
4455 1 : }
4456 1 : CATCH_END_SECTION()
4457 3 : }
4458 6 : }
4459 :
4460 :
4461 : // vim: ts=4 sw=4 et
|
