libs/multiprecision/test/test_eigen_interop.cpp

   1 //  Copyright 2012 John Maddock. Distributed under the Boost
   2 //  Software License, Version 1.0. (See accompanying file
   3 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
   4
   5 #include <boost/multiprecision/cpp_int.hpp>
   6 #include <boost/multiprecision/cpp_dec_float.hpp>
   7 #include <boost/multiprecision/cpp_bin_float.hpp>
   8 #include <iostream>
   9 #include <Eigen/Dense>
  10
  11 #include <boost/multiprecision/mpfr.hpp>
  12 #include <boost/multiprecision/logged_adaptor.hpp>
  13 #include <boost/multiprecision/gmp.hpp>
  14 #include <boost/multiprecision/mpc.hpp>
  15
  16 #include "test.hpp"
  17
  18 using namespace Eigen;
  19
  20 namespace Eigen {
  21 template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates>
  22 struct NumTraits<boost::multiprecision::number<Backend, ExpressionTemplates> >
  23 {
  24    typedef boost::multiprecision::number<Backend, ExpressionTemplates>                          self_type;
  25    typedef typename boost::multiprecision::scalar_result_from_possible_complex<self_type>::type Real;
  26    typedef self_type                                                                            NonInteger; // Not correct but we can't do much better??
  27    typedef double                                                                               Literal;
  28    typedef self_type                                                                            Nested;
  29    enum
  30    {
  31       IsComplex             = boost::multiprecision::number_category<self_type>::value == boost::multiprecision::number_kind_complex,
  32       IsInteger             = boost::multiprecision::number_category<self_type>::value == boost::multiprecision::number_kind_integer,
  33       ReadCost              = 1,
  34       AddCost               = 4,
  35       MulCost               = 8,
  36       IsSigned              = std::numeric_limits<self_type>::is_specialized ? std::numeric_limits<self_type>::is_signed : true,
  37       RequireInitialization = 1,
  38    };
  39    static Real epsilon()
  40    {
  41       return std::numeric_limits<Real>::epsilon();
  42    }
  43    static Real dummy_precision()
  44    {
  45       return sqrt(epsilon());
  46    }
  47    static Real highest()
  48    {
  49       return (std::numeric_limits<Real>::max)();
  50    }
  51    static Real lowest()
  52    {
  53       return (std::numeric_limits<Real>::min)();
  54    }
  55    static int digits10_imp(const boost::mpl::true_&)
  56    {
  57       return std::numeric_limits<Real>::digits10;
  58    }
  59    template <bool B>
  60    static int digits10_imp(const boost::mpl::bool_<B>&)
  61    {
  62       return Real::default_precision();
  63    }
  64    static int digits10()
  65    {
  66       return digits10_imp(boost::mpl::bool_ < std::numeric_limits<Real>::digits10 && (std::numeric_limits<Real>::digits10 != INT_MAX) ? true : false > ());
  67    }
  68 };
  69
  70 #define BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(A)                                                                                                                                                                       \
  71    template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, typename BinaryOp>                                                                                              \
  72    struct ScalarBinaryOpTraits<boost::multiprecision::number<Backend, ExpressionTemplates>, A, BinaryOp>                                                                                                           \
  73    {                                                                                                                                                                                                               \
  74       static_assert(boost::multiprecision::is_compatible_arithmetic_type<A, boost::multiprecision::number<Backend, ExpressionTemplates> >::value, "Interoperability with this arithmetic type is not supported."); \
  75       typedef boost::multiprecision::number<Backend, ExpressionTemplates> ReturnType;                                                                                                                              \
  76    };                                                                                                                                                                                                              \
  77    template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, typename BinaryOp>                                                                                              \
  78    struct ScalarBinaryOpTraits<A, boost::multiprecision::number<Backend, ExpressionTemplates>, BinaryOp>                                                                                                           \
  79    {                                                                                                                                                                                                               \
  80       static_assert(boost::multiprecision::is_compatible_arithmetic_type<A, boost::multiprecision::number<Backend, ExpressionTemplates> >::value, "Interoperability with this arithmetic type is not supported."); \
  81       typedef boost::multiprecision::number<Backend, ExpressionTemplates> ReturnType;                                                                                                                              \
  82    };
  83
  84 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(float)
  85 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(double)
  86 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(long double)
  87 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(char)
  88 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned char)
  89 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(signed char)
  90 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(short)
  91 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned short)
  92 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(int)
  93 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned int)
  94 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(long)
  95 BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned long)
  96 #if 0
  97    template<class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, class Backend2, boost::multiprecision::expression_template_option ExpressionTemplates2, typename BinaryOp>
  98    struct ScalarBinaryOpTraits<boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::number<Backend2, ExpressionTemplates2>, BinaryOp>
  99    {
 100       static_assert(
 101          boost::multiprecision::is_compatible_arithmetic_type<boost::multiprecision::number<Backend2, ExpressionTemplates2>, boost::multiprecision::number<Backend, ExpressionTemplates> >::value
 102          || boost::multiprecision::is_compatible_arithmetic_type<boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::number<Backend2, ExpressionTemplates2> >::value, "Interoperability with this arithmetic type is not supported.");
 103       typedef typename boost::mpl::if_c<boost::is_convertible<boost::multiprecision::number<Backend2, ExpressionTemplates2>, boost::multiprecision::number<Backend, ExpressionTemplates> >::value,
 104          boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::number<Backend2, ExpressionTemplates2> >::type ReturnType;
 105    };
 106
 107    template<unsigned D, typename BinaryOp>
 108    struct ScalarBinaryOpTraits<boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<D>, boost::multiprecision::et_on>, boost::multiprecision::mpfr_float, BinaryOp>
 109    {
 110       typedef boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<D>, boost::multiprecision::et_on> ReturnType;
 111    };
 112
 113    template<typename BinaryOp>
 114    struct ScalarBinaryOpTraits<boost::multiprecision::mpfr_float, boost::multiprecision::mpc_complex, BinaryOp>
 115    {
 116       typedef boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<0>, boost::multiprecision::et_on> ReturnType;
 117    };
 118
 119    template<class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, typename BinaryOp>
 120    struct ScalarBinaryOpTraits<boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::number<Backend, ExpressionTemplates>, BinaryOp>
 121    {
 122       typedef boost::multiprecision::number<Backend, ExpressionTemplates> ReturnType;
 123    };
 124 #endif
 125 template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, class tag, class Arg1, class Arg2, class Arg3, class Arg4, typename BinaryOp>
 126 struct ScalarBinaryOpTraits<boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, BinaryOp>
 127 {
 128    static_assert(boost::is_convertible<typename boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type, boost::multiprecision::number<Backend, ExpressionTemplates> >::value, "Interoperability with this arithmetic type is not supported.");
 129    typedef boost::multiprecision::number<Backend, ExpressionTemplates> ReturnType;
 130 };
 131
 132 template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, typename BinaryOp>
 133 struct ScalarBinaryOpTraits<boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, boost::multiprecision::number<Backend, ExpressionTemplates>, BinaryOp>
 134 {
 135    static_assert(boost::is_convertible<typename boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type, boost::multiprecision::number<Backend, ExpressionTemplates> >::value, "Interoperability with this arithmetic type is not supported.");
 136    typedef boost::multiprecision::number<Backend, ExpressionTemplates> ReturnType;
 137 };
 138
 139 } // namespace Eigen
 140
 141 template <class T>
 142 struct related_number
 143 {
 144    typedef T type;
 145 };
 146 /*
 147 template <>
 148 struct related_number<boost::multiprecision::cpp_bin_float_50>
 149 {
 150    typedef boost::multiprecision::cpp_bin_float_quad type;
 151 };
 152 template <>
 153 struct related_number<boost::multiprecision::cpp_dec_float_100>
 154 {
 155    typedef boost::multiprecision::cpp_dec_float_50 type;
 156 };*/
 157
 158 template <class Num>
 159 void example1()
 160 {
 161    // expected results first:
 162    Matrix<Num, 2, 2> r1, r2;
 163    r1 << 3, 5, 4, 8;
 164    r2 << -1, -1, 2, 0;
 165    Matrix<Num, 3, 1> r3;
 166    r3 << -1, -4, -6;
 167
 168    Matrix<Num, 2, 2> a;
 169    a << 1, 2, 3, 4;
 170    Matrix<Num, Dynamic, Dynamic> b(2, 2);
 171    b << 2, 3, 1, 4;
 172    std::cout << "a + b =\n"
 173              << a + b << std::endl;
 174    BOOST_CHECK_EQUAL(a + b, r1);
 175    std::cout << "a - b =\n"
 176              << a - b << std::endl;
 177    BOOST_CHECK_EQUAL(a - b, r2);
 178    std::cout << "Doing a += b;" << std::endl;
 179    a += b;
 180    std::cout << "Now a =\n"
 181              << a << std::endl;
 182    Matrix<Num, 3, 1> v(1, 2, 3);
 183    Matrix<Num, 3, 1> w(1, 0, 0);
 184    std::cout << "-v + w - v =\n"
 185              << -v + w - v << std::endl;
 186    BOOST_CHECK_EQUAL(-v + w - v, r3);
 187 }
 188
 189 template <class Num>
 190 void example2()
 191 {
 192    Matrix<Num, 2, 2> a;
 193    a << 1, 2, 3, 4;
 194    Matrix<Num, 3, 1> v(1, 2, 3);
 195    std::cout << "a * 2.5 =\n"
 196              << a * 2.5 << std::endl;
 197    std::cout << "0.1 * v =\n"
 198              << 0.1 * v << std::endl;
 199    std::cout << "Doing v *= 2;" << std::endl;
 200    v *= 2;
 201    std::cout << "Now v =\n"
 202              << v << std::endl;
 203    Num n(4);
 204    std::cout << "Doing v *= Num;" << std::endl;
 205    v *= n;
 206    std::cout << "Now v =\n"
 207              << v << std::endl;
 208    typedef typename related_number<Num>::type related_type;
 209    related_type                               r(6);
 210    std::cout << "Doing v *= RelatedType;" << std::endl;
 211    v *= r;
 212    std::cout << "Now v =\n"
 213              << v << std::endl;
 214    std::cout << "RelatedType * v =\n"
 215              << r * v << std::endl;
 216    std::cout << "Doing v *= RelatedType^2;" << std::endl;
 217    v *= r * r;
 218    std::cout << "Now v =\n"
 219              << v << std::endl;
 220    std::cout << "RelatedType^2 * v =\n"
 221              << r * r * v << std::endl;
 222
 223    static_assert(boost::is_same<typename Eigen::ScalarBinaryOpTraits<Num, related_type, Eigen::internal::scalar_product_op<Num, related_type> >::ReturnType, Num>::value, "Incorrect type.");
 224 }
 225
 226 template <class Num>
 227 void example3()
 228 {
 229    using namespace std;
 230    Matrix<Num, Dynamic, Dynamic> a = Matrix<Num, Dynamic, Dynamic>::Random(2, 2);
 231    cout << "Here is the matrix a\n"
 232         << a << endl;
 233    cout << "Here is the matrix a^T\n"
 234         << a.transpose() << endl;
 235    cout << "Here is the conjugate of a\n"
 236         << a.conjugate() << endl;
 237    cout << "Here is the matrix a^*\n"
 238         << a.adjoint() << endl;
 239 }
 240
 241 template <class Num>
 242 void example4()
 243 {
 244    Matrix<Num, 2, 2> mat;
 245    mat << 1, 2,
 246        3, 4;
 247    Matrix<Num, 2, 1> u(-1, 1), v(2, 0);
 248    std::cout << "Here is mat*mat:\n"
 249              << mat * mat << std::endl;
 250    std::cout << "Here is mat*u:\n"
 251              << mat * u << std::endl;
 252    std::cout << "Here is u^T*mat:\n"
 253              << u.transpose() * mat << std::endl;
 254    std::cout << "Here is u^T*v:\n"
 255              << u.transpose() * v << std::endl;
 256    std::cout << "Here is u*v^T:\n"
 257              << u * v.transpose() << std::endl;
 258    std::cout << "Let's multiply mat by itself" << std::endl;
 259    mat = mat * mat;
 260    std::cout << "Now mat is mat:\n"
 261              << mat << std::endl;
 262 }
 263
 264 template <class Num>
 265 void example5()
 266 {
 267    using namespace std;
 268    Matrix<Num, 3, 1> v(1, 2, 3);
 269    Matrix<Num, 3, 1> w(0, 1, 2);
 270    cout << "Dot product: " << v.dot(w) << endl;
 271    Num dp = v.adjoint() * w; // automatic conversion of the inner product to a scalar
 272    cout << "Dot product via a matrix product: " << dp << endl;
 273    cout << "Cross product:\n"
 274         << v.cross(w) << endl;
 275 }
 276
 277 template <class Num>
 278 void example6()
 279 {
 280    using namespace std;
 281    Matrix<Num, 2, 2> mat;
 282    mat << 1, 2,
 283        3, 4;
 284    cout << "Here is mat.sum():       " << mat.sum() << endl;
 285    cout << "Here is mat.prod():      " << mat.prod() << endl;
 286    cout << "Here is mat.mean():      " << mat.mean() << endl;
 287    cout << "Here is mat.minCoeff():  " << mat.minCoeff() << endl;
 288    cout << "Here is mat.maxCoeff():  " << mat.maxCoeff() << endl;
 289    cout << "Here is mat.trace():     " << mat.trace() << endl;
 290 }
 291
 292 template <class Num>
 293 void example7()
 294 {
 295    using namespace std;
 296
 297    Array<Num, Dynamic, Dynamic> m(2, 2);
 298
 299    // assign some values coefficient by coefficient
 300    m(0, 0) = 1.0;
 301    m(0, 1) = 2.0;
 302    m(1, 0) = 3.0;
 303    m(1, 1) = m(0, 1) + m(1, 0);
 304
 305    // print values to standard output
 306    cout << m << endl
 307         << endl;
 308
 309    // using the comma-initializer is also allowed
 310    m << 1.0, 2.0,
 311        3.0, 4.0;
 312
 313    // print values to standard output
 314    cout << m << endl;
 315 }
 316
 317 template <class Num>
 318 void example8()
 319 {
 320    using namespace std;
 321    Array<Num, Dynamic, Dynamic> a(3, 3);
 322    Array<Num, Dynamic, Dynamic> b(3, 3);
 323    a << 1, 2, 3,
 324        4, 5, 6,
 325        7, 8, 9;
 326    b << 1, 2, 3,
 327        1, 2, 3,
 328        1, 2, 3;
 329
 330    // Adding two arrays
 331    cout << "a + b = " << endl
 332         << a + b << endl
 333         << endl;
 334    // Subtracting a scalar from an array
 335    cout << "a - 2 = " << endl
 336         << a - 2 << endl;
 337 }
 338
 339 template <class Num>
 340 void example9()
 341 {
 342    using namespace std;
 343    Array<Num, Dynamic, Dynamic> a(2, 2);
 344    Array<Num, Dynamic, Dynamic> b(2, 2);
 345    a << 1, 2,
 346        3, 4;
 347    b << 5, 6,
 348        7, 8;
 349    cout << "a * b = " << endl
 350         << a * b << endl;
 351 }
 352
 353 template <class Num>
 354 void example10()
 355 {
 356    using namespace std;
 357    Array<Num, Dynamic, 1> a = Array<Num, Dynamic, 1>::Random(5);
 358    a *= 2;
 359    cout << "a =" << endl
 360         << a << endl;
 361    cout << "a.abs() =" << endl
 362         << a.abs() << endl;
 363    cout << "a.abs().sqrt() =" << endl
 364         << a.abs().sqrt() << endl;
 365    cout << "a.min(a.abs().sqrt()) =" << endl
 366         << a.std::min)(a.abs().sqrt()) << endl;
 367 }
 368
 369 template <class Num>
 370 void example11()
 371 {
 372    using namespace std;
 373    Matrix<Num, Dynamic, Dynamic> m(2, 2);
 374    Matrix<Num, Dynamic, Dynamic> n(2, 2);
 375    Matrix<Num, Dynamic, Dynamic> result(2, 2);
 376    m << 1, 2,
 377        3, 4;
 378    n << 5, 6,
 379        7, 8;
 380    result = m * n;
 381    cout << "-- Matrix m*n: --" << endl
 382         << result << endl
 383         << endl;
 384    result = m.array() * n.array();
 385    cout << "-- Array m*n: --" << endl
 386         << result << endl
 387         << endl;
 388    result = m.cwiseProduct(n);
 389    cout << "-- With cwiseProduct: --" << endl
 390         << result << endl
 391         << endl;
 392    result = m.array() + 4;
 393    cout << "-- Array m + 4: --" << endl
 394         << result << endl
 395         << endl;
 396 }
 397
 398 template <class Num>
 399 void example12()
 400 {
 401    using namespace std;
 402    Matrix<Num, Dynamic, Dynamic> m(2, 2);
 403    Matrix<Num, Dynamic, Dynamic> n(2, 2);
 404    Matrix<Num, Dynamic, Dynamic> result(2, 2);
 405    m << 1, 2,
 406        3, 4;
 407    n << 5, 6,
 408        7, 8;
 409
 410    result = (m.array() + 4).matrix() * m;
 411    cout << "-- Combination 1: --" << endl
 412         << result << endl
 413         << endl;
 414    result = (m.array() * n.array()).matrix() * m;
 415    cout << "-- Combination 2: --" << endl
 416         << result << endl
 417         << endl;
 418 }
 419
 420 template <class Num>
 421 void example13()
 422 {
 423    using namespace std;
 424    Matrix<Num, Dynamic, Dynamic> m(4, 4);
 425    m << 1, 2, 3, 4,
 426        5, 6, 7, 8,
 427        9, 10, 11, 12,
 428        13, 14, 15, 16;
 429    cout << "Block in the middle" << endl;
 430    cout << m.template block<2, 2>(1, 1) << endl
 431         << endl;
 432    for (int i = 1; i <= 3; ++i)
 433    {
 434       cout << "Block of size " << i << "x" << i << endl;
 435       cout << m.block(0, 0, i, i) << endl
 436            << endl;
 437    }
 438 }
 439
 440 template <class Num>
 441 void example14()
 442 {
 443    using namespace std;
 444    Array<Num, 2, 2> m;
 445    m << 1, 2,
 446        3, 4;
 447    Array<Num, 4, 4> a = Array<Num, 4, 4>::Constant(0.6);
 448    cout << "Here is the array a:" << endl
 449         << a << endl
 450         << endl;
 451    a.template block<2, 2>(1, 1) = m;
 452    cout << "Here is now a with m copied into its central 2x2 block:" << endl
 453         << a << endl
 454         << endl;
 455    a.block(0, 0, 2, 3) = a.block(2, 1, 2, 3);
 456    cout << "Here is now a with bottom-right 2x3 block copied into top-left 2x2 block:" << endl
 457         << a << endl
 458         << endl;
 459 }
 460
 461 template <class Num>
 462 void example15()
 463 {
 464    using namespace std;
 465    Eigen::Matrix<Num, Dynamic, Dynamic> m(3, 3);
 466    m << 1, 2, 3,
 467        4, 5, 6,
 468        7, 8, 9;
 469    cout << "Here is the matrix m:" << endl
 470         << m << endl;
 471    cout << "2nd Row: " << m.row(1) << endl;
 472    m.col(2) += 3 * m.col(0);
 473    cout << "After adding 3 times the first column into the third column, the matrix m is:\n";
 474    cout << m << endl;
 475 }
 476
 477 template <class Num>
 478 void example16()
 479 {
 480    using namespace std;
 481    Matrix<Num, 4, 4> m;
 482    m << 1, 2, 3, 4,
 483        5, 6, 7, 8,
 484        9, 10, 11, 12,
 485        13, 14, 15, 16;
 486    cout << "m.leftCols(2) =" << endl
 487         << m.leftCols(2) << endl
 488         << endl;
 489    cout << "m.bottomRows<2>() =" << endl
 490         << m.template bottomRows<2>() << endl
 491         << endl;
 492    m.topLeftCorner(1, 3) = m.bottomRightCorner(3, 1).transpose();
 493    cout << "After assignment, m = " << endl
 494         << m << endl;
 495 }
 496
 497 template <class Num>
 498 void example17()
 499 {
 500    using namespace std;
 501    Array<Num, Dynamic, 1> v(6);
 502    v << 1, 2, 3, 4, 5, 6;
 503    cout << "v.head(3) =" << endl
 504         << v.head(3) << endl
 505         << endl;
 506    cout << "v.tail<3>() = " << endl
 507         << v.template tail<3>() << endl
 508         << endl;
 509    v.segment(1, 4) *= 2;
 510    cout << "after 'v.segment(1,4) *= 2', v =" << endl
 511         << v << endl;
 512 }
 513
 514 template <class Num>
 515 void example18()
 516 {
 517    using namespace std;
 518    Matrix<Num, 2, 2> mat;
 519    mat << 1, 2,
 520        3, 4;
 521    cout << "Here is mat.sum():       " << mat.sum() << endl;
 522    cout << "Here is mat.prod():      " << mat.prod() << endl;
 523    cout << "Here is mat.mean():      " << mat.mean() << endl;
 524    cout << "Here is mat.minCoeff():  " << mat.minCoeff() << endl;
 525    cout << "Here is mat.maxCoeff():  " << mat.maxCoeff() << endl;
 526    cout << "Here is mat.trace():     " << mat.trace() << endl;
 527
 528    BOOST_CHECK_EQUAL(mat.sum(), 10);
 529    BOOST_CHECK_EQUAL(mat.prod(), 24);
 530    BOOST_CHECK_EQUAL(mat.mean(), Num(5) / 2);
 531    BOOST_CHECK_EQUAL(mat.minCoeff(), 1);
 532    BOOST_CHECK_EQUAL(mat.maxCoeff(), 4);
 533    BOOST_CHECK_EQUAL(mat.trace(), 5);
 534 }
 535
 536 template <class Num>
 537 void example18a()
 538 {
 539    using namespace std;
 540    Matrix<Num, 2, 2> mat;
 541    mat << 1, 2,
 542        3, 4;
 543    cout << "Here is mat.sum():       " << mat.sum() << endl;
 544    cout << "Here is mat.prod():      " << mat.prod() << endl;
 545    cout << "Here is mat.mean():      " << mat.mean() << endl;
 546    //cout << "Here is mat.minCoeff():  " << mat.minCoeff() << endl;
 547    //cout << "Here is mat.maxCoeff():  " << mat.maxCoeff() << endl;
 548    cout << "Here is mat.trace():     " << mat.trace() << endl;
 549 }
 550
 551 template <class Num>
 552 void example19()
 553 {
 554    using namespace std;
 555    Matrix<Num, Dynamic, 1>       v(2);
 556    Matrix<Num, Dynamic, Dynamic> m(2, 2), n(2, 2);
 557
 558    v << -1,
 559        2;
 560
 561    m << 1, -2,
 562        -3, 4;
 563    cout << "v.squaredNorm() = " << v.squaredNorm() << endl;
 564    cout << "v.norm() = " << v.norm() << endl;
 565    cout << "v.lpNorm<1>() = " << v.template lpNorm<1>() << endl;
 566    cout << "v.lpNorm<Infinity>() = " << v.template lpNorm<Infinity>() << endl;
 567    cout << endl;
 568    cout << "m.squaredNorm() = " << m.squaredNorm() << endl;
 569    cout << "m.norm() = " << m.norm() << endl;
 570    cout << "m.lpNorm<1>() = " << m.template lpNorm<1>() << endl;
 571    cout << "m.lpNorm<Infinity>() = " << m.template lpNorm<Infinity>() << endl;
 572 }
 573
 574 template <class Num>
 575 void example20()
 576 {
 577    using namespace std;
 578    Matrix<Num, 3, 3> A;
 579    Matrix<Num, 3, 1> b;
 580    A << 1, 2, 3, 4, 5, 6, 7, 8, 10;
 581    b << 3, 3, 4;
 582    cout << "Here is the matrix A:\n"
 583         << A << endl;
 584    cout << "Here is the vector b:\n"
 585         << b << endl;
 586    Matrix<Num, 3, 1> x = A.colPivHouseholderQr().solve(b);
 587    cout << "The solution is:\n"
 588         << x << endl;
 589 }
 590
 591 template <class Num>
 592 void example21()
 593 {
 594    using namespace std;
 595    Matrix<Num, 2, 2> A, b;
 596    A << 2, -1, -1, 3;
 597    b << 1, 2, 3, 1;
 598    cout << "Here is the matrix A:\n"
 599         << A << endl;
 600    cout << "Here is the right hand side b:\n"
 601         << b << endl;
 602    Matrix<Num, 2, 2> x = A.ldlt().solve(b);
 603    cout << "The solution is:\n"
 604         << x << endl;
 605 }
 606
 607 template <class Num>
 608 void example22()
 609 {
 610    using namespace std;
 611    Matrix<Num, Dynamic, Dynamic> A              = Matrix<Num, Dynamic, Dynamic>::Random(100, 100);
 612    Matrix<Num, Dynamic, Dynamic> b              = Matrix<Num, Dynamic, Dynamic>::Random(100, 50);
 613    Matrix<Num, Dynamic, Dynamic> x              = A.fullPivLu().solve(b);
 614    Matrix<Num, Dynamic, Dynamic> axmb           = A * x - b;
 615    double                        relative_error = static_cast<double>(abs(axmb.norm() / b.norm())); // norm() is L2 norm
 616    cout << "norm1 = " << axmb.norm() << endl;
 617    cout << "norm2 = " << b.norm() << endl;
 618    cout << "The relative error is:\n"
 619         << relative_error << endl;
 620 }
 621
 622 template <class Num>
 623 void example23()
 624 {
 625    using namespace std;
 626    Matrix<Num, 2, 2> A;
 627    A << 1, 2, 2, 3;
 628    cout << "Here is the matrix A:\n"
 629         << A << endl;
 630    SelfAdjointEigenSolver<Matrix<Num, 2, 2> > eigensolver(A);
 631    if (eigensolver.info() != Success)
 632    {
 633       std::cout << "Eigenvalue solver failed!" << endl;
 634    }
 635    else
 636    {
 637       cout << "The eigenvalues of A are:\n"
 638            << eigensolver.eigenvalues() << endl;
 639       cout << "Here's a matrix whose columns are eigenvectors of A \n"
 640            << "corresponding to these eigenvalues:\n"
 641            << eigensolver.eigenvectors() << endl;
 642    }
 643 }
 644
 645 template <class Num>
 646 void example24()
 647 {
 648    using namespace std;
 649    Matrix<Num, 3, 3> A;
 650    A << 1, 2, 1,
 651        2, 1, 0,
 652        -1, 1, 2;
 653    cout << "Here is the matrix A:\n"
 654         << A << endl;
 655    cout << "The determinant of A is " << A.determinant() << endl;
 656    cout << "The inverse of A is:\n"
 657         << A.inverse() << endl;
 658 }
 659
 660 template <class Num>
 661 void test_integer_type()
 662 {
 663    example1<Num>();
 664    //example2<Num>();
 665    example18<Num>();
 666 }
 667
 668 template <class Num>
 669 void test_float_type()
 670 {
 671    std::cout << "Epsilon    = " << Eigen::NumTraits<Num>::epsilon() << std::endl;
 672    std::cout << "Dummy Prec = " << Eigen::NumTraits<Num>::dummy_precision() << std::endl;
 673    std::cout << "Highest    = " << Eigen::NumTraits<Num>::highest() << std::endl;
 674    std::cout << "Lowest     = " << Eigen::NumTraits<Num>::lowest() << std::endl;
 675    std::cout << "Digits10   = " << Eigen::NumTraits<Num>::digits10() << std::endl;
 676
 677    example1<Num>();
 678    example2<Num>();
 679    example4<Num>();
 680    example5<Num>();
 681    example6<Num>();
 682    example7<Num>();
 683    example8<Num>();
 684    example9<Num>();
 685    example10<Num>();
 686    example11<Num>();
 687    example12<Num>();
 688    example13<Num>();
 689    example14<Num>();
 690    example15<Num>();
 691    example16<Num>();
 692    example17<Num>();
 693    example18<Num>();
 694    example19<Num>();
 695    example20<Num>();
 696    example21<Num>();
 697    example22<Num>();
 698    example23<Num>();
 699    example24<Num>();
 700 }
 701
 702 template <class Num>
 703 void test_complex_type()
 704 {
 705    std::cout << "Epsilon    = " << Eigen::NumTraits<Num>::epsilon() << std::endl;
 706    std::cout << "Dummy Prec = " << Eigen::NumTraits<Num>::dummy_precision() << std::endl;
 707    std::cout << "Highest    = " << Eigen::NumTraits<Num>::highest() << std::endl;
 708    std::cout << "Lowest     = " << Eigen::NumTraits<Num>::lowest() << std::endl;
 709    std::cout << "Digits10   = " << Eigen::NumTraits<Num>::digits10() << std::endl;
 710
 711    example1<Num>();
 712    example2<Num>();
 713    example3<Num>();
 714    example4<Num>();
 715    example5<Num>();
 716    example7<Num>();
 717    example8<Num>();
 718    example9<Num>();
 719    example11<Num>();
 720    example12<Num>();
 721    example13<Num>();
 722    example14<Num>();
 723    example15<Num>();
 724    example16<Num>();
 725    example17<Num>();
 726    example18a<Num>();
 727    example19<Num>();
 728    example20<Num>();
 729    example21<Num>();
 730    example22<Num>();
 731    // example23<Num>();  //requires comparisons.
 732    example24<Num>();
 733 }
 734
 735 namespace boost {
 736 namespace multiprecision {
 737
 738 template <unsigned D>
 739 inline void log_postfix_event(const mpc_complex_backend<D>& val, const char* event_description)
 740 {
 741    if (mpfr_nan_p(mpc_realref(val.data())))
 742    {
 743       std::cout << "Found a NaN! " << event_description << std::endl;
 744    }
 745 }
 746
 747 }
 748 } // namespace boost::multiprecision
 749
 750 int main()
 751 {
 752    using namespace boost::multiprecision;
 753    test_complex_type<mpc_complex>();
 754 #if 0
 755    test_integer_type<int>();
 756
 757    test_float_type<double>();
 758    test_complex_type<std::complex<double> >();
 759
 760    test_float_type<boost::multiprecision::cpp_dec_float_100>();
 761    test_float_type<boost::multiprecision::cpp_bin_float_50>();
 762    test_float_type<boost::multiprecision::mpfr_float>();
 763
 764    test_integer_type<boost::multiprecision::int256_t>();
 765    test_integer_type<boost::multiprecision::cpp_int>();
 766    test_integer_type<boost::multiprecision::cpp_rational>();
 767    test_integer_type<boost::multiprecision::mpz_int>();
 768    test_integer_type<boost::multiprecision::mpq_rational>();
 769
 770 #endif
 771    return 0;
 772 }