New scoring functor to select the pivot.

This is can be useful for non-floating point scalars, where choosing the biggest element is generally not the best choice.

1 files changed, 61 insertions, 0 deletions

diff --git a/doc/CustomizingEigen.dox b/doc/CustomizingEigen.dox
index 0850863aa..cb25f4ec9 100644
--- a/doc/CustomizingEigen.dox
+++ b/doc/CustomizingEigen.dox
@@ -157,6 +157,67 @@ inline adouble abs2(const adouble& x)  { return x*x; }
 #endif // ADOLCSUPPORT_H
 \endcode
 
+This other example adds support for the \c mpq_class type from <a href="https://gmplib.org/">GMP</a>. It shows in particular how to change the way Eigen picks the best pivot during LU factorization. It selects the coefficient with the highest score, where the score is by default the absolute value of a number, but we can define a different score, for instance to prefer pivots with a more compact representation (this is an example, not a recommendation). Note that the scores should always be non-negative and only zero is allowed to have a score of zero. Also, this can interact badly with thresholds for inexact scalar types.
+
+\code
+#include <gmpxx.h>
+#include <Eigen/Core>
+#include <boost/operators.hpp>
+
+namespace Eigen {
+  template<class> struct NumTraits;
+  template<> struct NumTraits<mpq_class>
+  {
+    typedef mpq_class Real;
+    typedef mpq_class NonInteger;
+    typedef mpq_class Nested;
+
+    static inline Real epsilon() { return 0; }
+    static inline Real dummy_precision() { return 0; }
+
+    enum {
+      IsInteger = 0,
+      IsSigned = 1,
+      IsComplex = 0,
+      RequireInitialization = 1,
+      ReadCost = 6,
+      AddCost = 150,
+      MulCost = 100
+    };
+  };
+
+  namespace internal {
+    template<>
+      struct significant_decimals_impl<mpq_class>
+      {
+	// Infinite precision when printing
+	static inline int run() { return 0; }
+      };
+
+    template<> struct scalar_score_coeff_op<mpq_class> {
+      struct result_type : boost::totally_ordered1<result_type> {
+	std::size_t len;
+	result_type(int i = 0) : len(i) {} // Eigen uses Score(0) and Score()
+	result_type(mpq_class const& q) :
+	  len(mpz_size(q.get_num_mpz_t())+
+	      mpz_size(q.get_den_mpz_t())-1) {}
+	friend bool operator<(result_type x, result_type y) {
+	  // 0 is the worst possible pivot
+	  if (x.len == 0) return y.len > 0;
+	  if (y.len == 0) return false;
+	  // Prefer a pivot with a small representation
+	  return x.len > y.len;
+	}
+	friend bool operator==(result_type x, result_type y) {
+	  // Only used to test if the score is 0
+	  return x.len == y.len;
+	}
+      };
+      result_type operator()(mpq_class const& x) const { return x; }
+    };
+  }
+}
+\endcode
 
 \sa \ref TopicPreprocessorDirectives