1 ///////////////////////////////////////////////////////////////////////////////
2 // Copyright 2018 John Maddock. Distributed under the Boost
3 // Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #include <boost/multiprecision/eigen.hpp>
7 #include <iostream>
8 #include <Eigen/Dense>
9 #include "test.hpp"
10
11 using namespace Eigen;
12
13 template <class T>
14 struct related_number
15 {
16 typedef T type;
17 };
18
19 template <class Num>
example1()20 void example1()
21 {
22 // expected results first:
23 Matrix<Num, 2, 2> r1, r2;
24 r1 << 3, 5, 4, 8;
25 r2 << -1, -1, 2, 0;
26 Matrix<Num, 3, 1> r3;
27 r3 << -1, -4, -6;
28
29 Matrix<Num, 2, 2> a;
30 a << 1, 2, 3, 4;
31 Matrix<Num, Dynamic, Dynamic> b(2, 2);
32 b << 2, 3, 1, 4;
33 std::cout << "a + b =\n"
34 << a + b << std::endl;
35 BOOST_CHECK_EQUAL(a + b, r1);
36 std::cout << "a - b =\n"
37 << a - b << std::endl;
38 BOOST_CHECK_EQUAL(a - b, r2);
39 std::cout << "Doing a += b;" << std::endl;
40 a += b;
41 std::cout << "Now a =\n"
42 << a << std::endl;
43 Matrix<Num, 3, 1> v(1, 2, 3);
44 Matrix<Num, 3, 1> w(1, 0, 0);
45 std::cout << "-v + w - v =\n"
46 << -v + w - v << std::endl;
47 BOOST_CHECK_EQUAL(-v + w - v, r3);
48 }
49
50 template <class Num>
example2()51 void example2()
52 {
53 Matrix<Num, 2, 2> a;
54 a << 1, 2, 3, 4;
55 Matrix<Num, 3, 1> v(1, 2, 3);
56 std::cout << "a * 2.5 =\n"
57 << a * 2.5 << std::endl;
58 std::cout << "0.1 * v =\n"
59 << 0.1 * v << std::endl;
60 std::cout << "Doing v *= 2;" << std::endl;
61 v *= 2;
62 std::cout << "Now v =\n"
63 << v << std::endl;
64 Num n(4);
65 std::cout << "Doing v *= Num;" << std::endl;
66 v *= n;
67 std::cout << "Now v =\n"
68 << v << std::endl;
69 typedef typename related_number<Num>::type related_type;
70 related_type r(6);
71 std::cout << "Doing v *= RelatedType;" << std::endl;
72 v *= r;
73 std::cout << "Now v =\n"
74 << v << std::endl;
75 std::cout << "RelatedType * v =\n"
76 << r * v << std::endl;
77 std::cout << "Doing v *= RelatedType^2;" << std::endl;
78 v *= r * r;
79 std::cout << "Now v =\n"
80 << v << std::endl;
81 std::cout << "RelatedType^2 * v =\n"
82 << r * r * v << std::endl;
83 }
84
85 template <class Num>
example3()86 void example3()
87 {
88 using namespace std;
89 Matrix<Num, Dynamic, Dynamic> a = Matrix<Num, Dynamic, Dynamic>::Random(2, 2);
90 cout << "Here is the matrix a\n"
91 << a << endl;
92 cout << "Here is the matrix a^T\n"
93 << a.transpose() << endl;
94 cout << "Here is the conjugate of a\n"
95 << a.conjugate() << endl;
96 cout << "Here is the matrix a^*\n"
97 << a.adjoint() << endl;
98 }
99
100 template <class Num>
example4()101 void example4()
102 {
103 Matrix<Num, 2, 2> mat;
104 mat << 1, 2,
105 3, 4;
106 Matrix<Num, 2, 1> u(-1, 1), v(2, 0);
107 std::cout << "Here is mat*mat:\n"
108 << mat * mat << std::endl;
109 std::cout << "Here is mat*u:\n"
110 << mat * u << std::endl;
111 std::cout << "Here is u^T*mat:\n"
112 << u.transpose() * mat << std::endl;
113 std::cout << "Here is u^T*v:\n"
114 << u.transpose() * v << std::endl;
115 std::cout << "Here is u*v^T:\n"
116 << u * v.transpose() << std::endl;
117 std::cout << "Let's multiply mat by itself" << std::endl;
118 mat = mat * mat;
119 std::cout << "Now mat is mat:\n"
120 << mat << std::endl;
121 }
122
123 template <class Num>
example5()124 void example5()
125 {
126 using namespace std;
127 Matrix<Num, 3, 1> v(1, 2, 3);
128 Matrix<Num, 3, 1> w(0, 1, 2);
129 cout << "Dot product: " << v.dot(w) << endl;
130 Num dp = v.adjoint() * w; // automatic conversion of the inner product to a scalar
131 cout << "Dot product via a matrix product: " << dp << endl;
132 cout << "Cross product:\n"
133 << v.cross(w) << endl;
134 }
135
136 template <class Num>
example6()137 void example6()
138 {
139 using namespace std;
140 Matrix<Num, 2, 2> mat;
141 mat << 1, 2,
142 3, 4;
143 cout << "Here is mat.sum(): " << mat.sum() << endl;
144 cout << "Here is mat.prod(): " << mat.prod() << endl;
145 cout << "Here is mat.mean(): " << mat.mean() << endl;
146 cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl;
147 cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl;
148 cout << "Here is mat.trace(): " << mat.trace() << endl;
149 }
150
151 template <class Num>
example7()152 void example7()
153 {
154 using namespace std;
155
156 Array<Num, Dynamic, Dynamic> m(2, 2);
157
158 // assign some values coefficient by coefficient
159 m(0, 0) = 1.0;
160 m(0, 1) = 2.0;
161 m(1, 0) = 3.0;
162 m(1, 1) = m(0, 1) + m(1, 0);
163
164 // print values to standard output
165 cout << m << endl
166 << endl;
167
168 // using the comma-initializer is also allowed
169 m << 1.0, 2.0,
170 3.0, 4.0;
171
172 // print values to standard output
173 cout << m << endl;
174 }
175
176 template <class Num>
example8()177 void example8()
178 {
179 using namespace std;
180 Array<Num, Dynamic, Dynamic> a(3, 3);
181 Array<Num, Dynamic, Dynamic> b(3, 3);
182 a << 1, 2, 3,
183 4, 5, 6,
184 7, 8, 9;
185 b << 1, 2, 3,
186 1, 2, 3,
187 1, 2, 3;
188
189 // Adding two arrays
190 cout << "a + b = " << endl
191 << a + b << endl
192 << endl;
193 // Subtracting a scalar from an array
194 cout << "a - 2 = " << endl
195 << a - 2 << endl;
196 }
197
198 template <class Num>
example9()199 void example9()
200 {
201 using namespace std;
202 Array<Num, Dynamic, Dynamic> a(2, 2);
203 Array<Num, Dynamic, Dynamic> b(2, 2);
204 a << 1, 2,
205 3, 4;
206 b << 5, 6,
207 7, 8;
208 cout << "a * b = " << endl
209 << a * b << endl;
210 }
211
212 template <class Num>
example10()213 void example10()
214 {
215 using namespace std;
216 Array<Num, Dynamic, 1> a = Array<Num, Dynamic, 1>::Random(5);
217 a *= 2;
218 cout << "a =" << endl
219 << a << endl;
220 cout << "a.abs() =" << endl
221 << a.abs() << endl;
222 cout << "a.abs().sqrt() =" << endl
223 << a.abs().sqrt() << endl;
224 cout << "a.min(a.abs().sqrt()) =" << endl
225 << (a.min)(a.abs().sqrt()) << endl;
226 }
227
228 template <class Num>
example11()229 void example11()
230 {
231 using namespace std;
232 Matrix<Num, Dynamic, Dynamic> m(2, 2);
233 Matrix<Num, Dynamic, Dynamic> n(2, 2);
234 Matrix<Num, Dynamic, Dynamic> result(2, 2);
235 m << 1, 2,
236 3, 4;
237 n << 5, 6,
238 7, 8;
239 result = m * n;
240 cout << "-- Matrix m*n: --" << endl
241 << result << endl
242 << endl;
243 result = m.array() * n.array();
244 cout << "-- Array m*n: --" << endl
245 << result << endl
246 << endl;
247 result = m.cwiseProduct(n);
248 cout << "-- With cwiseProduct: --" << endl
249 << result << endl
250 << endl;
251 result = m.array() + 4;
252 cout << "-- Array m + 4: --" << endl
253 << result << endl
254 << endl;
255 }
256
257 template <class Num>
example12()258 void example12()
259 {
260 using namespace std;
261 Matrix<Num, Dynamic, Dynamic> m(2, 2);
262 Matrix<Num, Dynamic, Dynamic> n(2, 2);
263 Matrix<Num, Dynamic, Dynamic> result(2, 2);
264 m << 1, 2,
265 3, 4;
266 n << 5, 6,
267 7, 8;
268
269 result = (m.array() + 4).matrix() * m;
270 cout << "-- Combination 1: --" << endl
271 << result << endl
272 << endl;
273 result = (m.array() * n.array()).matrix() * m;
274 cout << "-- Combination 2: --" << endl
275 << result << endl
276 << endl;
277 }
278
279 template <class Num>
example13()280 void example13()
281 {
282 using namespace std;
283 Matrix<Num, Dynamic, Dynamic> m(4, 4);
284 m << 1, 2, 3, 4,
285 5, 6, 7, 8,
286 9, 10, 11, 12,
287 13, 14, 15, 16;
288 cout << "Block in the middle" << endl;
289 cout << m.template block<2, 2>(1, 1) << endl
290 << endl;
291 for (int i = 1; i <= 3; ++i)
292 {
293 cout << "Block of size " << i << "x" << i << endl;
294 cout << m.block(0, 0, i, i) << endl
295 << endl;
296 }
297 }
298
299 template <class Num>
example14()300 void example14()
301 {
302 using namespace std;
303 Array<Num, 2, 2> m;
304 m << 1, 2,
305 3, 4;
306 Array<Num, 4, 4> a = Array<Num, 4, 4>::Constant(0.6);
307 cout << "Here is the array a:" << endl
308 << a << endl
309 << endl;
310 a.template block<2, 2>(1, 1) = m;
311 cout << "Here is now a with m copied into its central 2x2 block:" << endl
312 << a << endl
313 << endl;
314 a.block(0, 0, 2, 3) = a.block(2, 1, 2, 3);
315 cout << "Here is now a with bottom-right 2x3 block copied into top-left 2x2 block:" << endl
316 << a << endl
317 << endl;
318 }
319
320 template <class Num>
example15()321 void example15()
322 {
323 using namespace std;
324 Eigen::Matrix<Num, Dynamic, Dynamic> m(3, 3);
325 m << 1, 2, 3,
326 4, 5, 6,
327 7, 8, 9;
328 cout << "Here is the matrix m:" << endl
329 << m << endl;
330 cout << "2nd Row: " << m.row(1) << endl;
331 m.col(2) += 3 * m.col(0);
332 cout << "After adding 3 times the first column into the third column, the matrix m is:\n";
333 cout << m << endl;
334 }
335
336 template <class Num>
example16()337 void example16()
338 {
339 using namespace std;
340 Matrix<Num, 4, 4> m;
341 m << 1, 2, 3, 4,
342 5, 6, 7, 8,
343 9, 10, 11, 12,
344 13, 14, 15, 16;
345 cout << "m.leftCols(2) =" << endl
346 << m.leftCols(2) << endl
347 << endl;
348 cout << "m.bottomRows<2>() =" << endl
349 << m.template bottomRows<2>() << endl
350 << endl;
351 m.topLeftCorner(1, 3) = m.bottomRightCorner(3, 1).transpose();
352 cout << "After assignment, m = " << endl
353 << m << endl;
354 }
355
356 template <class Num>
example17()357 void example17()
358 {
359 using namespace std;
360 Array<Num, Dynamic, 1> v(6);
361 v << 1, 2, 3, 4, 5, 6;
362 cout << "v.head(3) =" << endl
363 << v.head(3) << endl
364 << endl;
365 cout << "v.tail<3>() = " << endl
366 << v.template tail<3>() << endl
367 << endl;
368 v.segment(1, 4) *= 2;
369 cout << "after 'v.segment(1,4) *= 2', v =" << endl
370 << v << endl;
371 }
372
373 template <class Num>
example18()374 void example18()
375 {
376 using namespace std;
377 Matrix<Num, 2, 2> mat;
378 mat << 1, 2,
379 3, 4;
380 cout << "Here is mat.sum(): " << mat.sum() << endl;
381 cout << "Here is mat.prod(): " << mat.prod() << endl;
382 cout << "Here is mat.mean(): " << mat.mean() << endl;
383 cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl;
384 cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl;
385 cout << "Here is mat.trace(): " << mat.trace() << endl;
386
387 BOOST_CHECK_EQUAL(mat.sum(), 10);
388 BOOST_CHECK_EQUAL(mat.prod(), 24);
389 BOOST_CHECK_EQUAL(mat.mean(), Num(5) / 2);
390 BOOST_CHECK_EQUAL(mat.minCoeff(), 1);
391 BOOST_CHECK_EQUAL(mat.maxCoeff(), 4);
392 BOOST_CHECK_EQUAL(mat.trace(), 5);
393 }
394
395 template <class Num>
example18a()396 void example18a()
397 {
398 using namespace std;
399 Matrix<Num, 2, 2> mat;
400 mat << 1, 2,
401 3, 4;
402 cout << "Here is mat.sum(): " << mat.sum() << endl;
403 cout << "Here is mat.prod(): " << mat.prod() << endl;
404 cout << "Here is mat.mean(): " << mat.mean() << endl;
405 //cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl;
406 //cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl;
407 cout << "Here is mat.trace(): " << mat.trace() << endl;
408 }
409
410 template <class Num>
example19()411 void example19()
412 {
413 using namespace std;
414 Matrix<Num, Dynamic, 1> v(2);
415 Matrix<Num, Dynamic, Dynamic> m(2, 2), n(2, 2);
416
417 v << -1,
418 2;
419
420 m << 1, -2,
421 -3, 4;
422 cout << "v.squaredNorm() = " << v.squaredNorm() << endl;
423 cout << "v.norm() = " << v.norm() << endl;
424 cout << "v.lpNorm<1>() = " << v.template lpNorm<1>() << endl;
425 cout << "v.lpNorm<Infinity>() = " << v.template lpNorm<Infinity>() << endl;
426 cout << endl;
427 cout << "m.squaredNorm() = " << m.squaredNorm() << endl;
428 cout << "m.norm() = " << m.norm() << endl;
429 cout << "m.lpNorm<1>() = " << m.template lpNorm<1>() << endl;
430 cout << "m.lpNorm<Infinity>() = " << m.template lpNorm<Infinity>() << endl;
431 }
432
433 template <class Num>
example20()434 void example20()
435 {
436 using namespace std;
437 Matrix<Num, 3, 3> A;
438 Matrix<Num, 3, 1> b;
439 A << 1, 2, 3, 4, 5, 6, 7, 8, 10;
440 b << 3, 3, 4;
441 cout << "Here is the matrix A:\n"
442 << A << endl;
443 cout << "Here is the vector b:\n"
444 << b << endl;
445 Matrix<Num, 3, 1> x = A.colPivHouseholderQr().solve(b);
446 cout << "The solution is:\n"
447 << x << endl;
448 }
449
450 template <class Num>
example21()451 void example21()
452 {
453 using namespace std;
454 Matrix<Num, 2, 2> A, b;
455 A << 2, -1, -1, 3;
456 b << 1, 2, 3, 1;
457 cout << "Here is the matrix A:\n"
458 << A << endl;
459 cout << "Here is the right hand side b:\n"
460 << b << endl;
461 Matrix<Num, 2, 2> x = A.ldlt().solve(b);
462 cout << "The solution is:\n"
463 << x << endl;
464 }
465
466 template <class Num>
example22()467 void example22()
468 {
469 using namespace std;
470 Matrix<Num, Dynamic, Dynamic> A = Matrix<Num, Dynamic, Dynamic>::Random(100, 100);
471 Matrix<Num, Dynamic, Dynamic> b = Matrix<Num, Dynamic, Dynamic>::Random(100, 50);
472 Matrix<Num, Dynamic, Dynamic> x = A.fullPivLu().solve(b);
473 Matrix<Num, Dynamic, Dynamic> axmb = A * x - b;
474 double relative_error = static_cast<double>(abs(axmb.norm() / b.norm())); // norm() is L2 norm
475 cout << "norm1 = " << axmb.norm() << endl;
476 cout << "norm2 = " << b.norm() << endl;
477 cout << "The relative error is:\n"
478 << relative_error << endl;
479 }
480
481 template <class Num>
example23()482 void example23()
483 {
484 using namespace std;
485 Matrix<Num, 2, 2> A;
486 A << 1, 2, 2, 3;
487 cout << "Here is the matrix A:\n"
488 << A << endl;
489 SelfAdjointEigenSolver<Matrix<Num, 2, 2> > eigensolver(A);
490 if (eigensolver.info() != Success)
491 {
492 std::cout << "Eigenvalue solver failed!" << endl;
493 }
494 else
495 {
496 cout << "The eigenvalues of A are:\n"
497 << eigensolver.eigenvalues() << endl;
498 cout << "Here's a matrix whose columns are eigenvectors of A \n"
499 << "corresponding to these eigenvalues:\n"
500 << eigensolver.eigenvectors() << endl;
501 }
502 }
503
504 template <class Num>
example24()505 void example24()
506 {
507 using namespace std;
508 Matrix<Num, 3, 3> A;
509 A << 1, 2, 1,
510 2, 1, 0,
511 -1, 1, 2;
512 cout << "Here is the matrix A:\n"
513 << A << endl;
514 cout << "The determinant of A is " << A.determinant() << endl;
515 cout << "The inverse of A is:\n"
516 << A.inverse() << endl;
517 }
518
519 template <class Num>
test_integer_type()520 void test_integer_type()
521 {
522 example1<Num>();
523 //example2<Num>();
524 example18<Num>();
525 }
526
527 template <class Num>
test_float_type()528 void test_float_type()
529 {
530 std::cout << "Epsilon = " << Eigen::NumTraits<Num>::epsilon() << std::endl;
531 std::cout << "Dummy Prec = " << Eigen::NumTraits<Num>::dummy_precision() << std::endl;
532 std::cout << "Highest = " << Eigen::NumTraits<Num>::highest() << std::endl;
533 std::cout << "Lowest = " << Eigen::NumTraits<Num>::lowest() << std::endl;
534 std::cout << "Digits10 = " << Eigen::NumTraits<Num>::digits10() << std::endl;
535
536 example1<Num>();
537 example2<Num>();
538 example4<Num>();
539 example5<Num>();
540 example6<Num>();
541 example7<Num>();
542 example8<Num>();
543 example9<Num>();
544 example10<Num>();
545 example11<Num>();
546 example12<Num>();
547 example13<Num>();
548 example14<Num>();
549 example15<Num>();
550 example16<Num>();
551 example17<Num>();
552 /*
553 example18<Num>();
554 example19<Num>();
555 example20<Num>();
556 example21<Num>();
557 example22<Num>();
558 example23<Num>();
559 example24<Num>();
560 */
561 }
562
563 template <class Num>
test_float_type_2()564 void test_float_type_2()
565 {
566 std::cout << "Epsilon = " << Eigen::NumTraits<Num>::epsilon() << std::endl;
567 std::cout << "Dummy Prec = " << Eigen::NumTraits<Num>::dummy_precision() << std::endl;
568 std::cout << "Highest = " << Eigen::NumTraits<Num>::highest() << std::endl;
569 std::cout << "Lowest = " << Eigen::NumTraits<Num>::lowest() << std::endl;
570 std::cout << "Digits10 = " << Eigen::NumTraits<Num>::digits10() << std::endl;
571
572 example18<Num>();
573 example19<Num>();
574 example20<Num>();
575 example21<Num>();
576
577 //example22<Num>();
578 //example23<Num>();
579 //example24<Num>();
580 }
581
582 template <class Num>
test_float_type_3()583 void test_float_type_3()
584 {
585 std::cout << "Epsilon = " << Eigen::NumTraits<Num>::epsilon() << std::endl;
586 std::cout << "Dummy Prec = " << Eigen::NumTraits<Num>::dummy_precision() << std::endl;
587 std::cout << "Highest = " << Eigen::NumTraits<Num>::highest() << std::endl;
588 std::cout << "Lowest = " << Eigen::NumTraits<Num>::lowest() << std::endl;
589 std::cout << "Digits10 = " << Eigen::NumTraits<Num>::digits10() << std::endl;
590
591 example22<Num>();
592 example23<Num>();
593 example24<Num>();
594 }
595
596 template <class Num>
test_complex_type()597 void test_complex_type()
598 {
599 std::cout << "Epsilon = " << Eigen::NumTraits<Num>::epsilon() << std::endl;
600 std::cout << "Dummy Prec = " << Eigen::NumTraits<Num>::dummy_precision() << std::endl;
601 std::cout << "Highest = " << Eigen::NumTraits<Num>::highest() << std::endl;
602 std::cout << "Lowest = " << Eigen::NumTraits<Num>::lowest() << std::endl;
603 std::cout << "Digits10 = " << Eigen::NumTraits<Num>::digits10() << std::endl;
604
605 example1<Num>();
606 example2<Num>();
607 example3<Num>();
608 example4<Num>();
609 example5<Num>();
610 example7<Num>();
611 example8<Num>();
612 example9<Num>();
613 example11<Num>();
614 example12<Num>();
615 example13<Num>();
616 example14<Num>();
617 example15<Num>();
618 example16<Num>();
619 example17<Num>();
620 example18a<Num>();
621 example19<Num>();
622 example20<Num>();
623 example21<Num>();
624 example22<Num>();
625 // example23<Num>(); //requires comparisons.
626 example24<Num>();
627 }
628