Re: [eigen] Specializing max_coeff_visitor for some number types

[ Thread Index | Date Index | More Archives ]

On Tue, 8 Apr 2014, Christoph Hertzberg wrote:

On 07.04.2014 20:33, Marc Glisse wrote:
The attached patch doesn't change the result of make check (adolc is
still broken, there is a missing -lGL somewhere, and parmetis doesn't
chose any MPI so finds none).

The adolc test is passing at my machine (that is with current devel, I have not tested your patch yet):

Building CXX object unsupported/test/CMakeFiles/forward_adolc.dir/forward_adolc.cpp.o
In file included from /data/repos/eigen-pristine/unsupported/test/forward_adolc.cpp:14:0:
/data/repos/eigen-pristine/unsupported/Eigen/AdolcForward:81:29: error: ʽadoubleʼ is not a member of ʽadtlʼ
 template<> struct NumTraits<adtl::adouble>
/data/repos/eigen-pristine/unsupported/Eigen/AdolcForward:81:29: note: suggested alternative:
In file included from /data/repos/eigen-pristine/unsupported/Eigen/AdolcForward:28:0,
                 from /data/repos/eigen-pristine/unsupported/test/forward_adolc.cpp:14:
/usr/include/adolc/adouble.h:350:24: note:   ʽadoubleʼ
 class ADOLC_DLL_EXPORT adouble:public badouble {

I have libadolc-dev 2.5.0-3 provided by debian. I ran cmake /path/to/eigen then make -j8 -k check.

I don't have METIS installed and strangely, OpenGL is not detected here.

PASTIX is the only thing that is missing here, but since I am only working with small dense matrices...

As you can see in LU, there is a subtlety I hadn't thought of: max(abs)
was returning a RealScalar, not a Scalar. That's not a big problem when
we only want to compare it with 0.

I think that m_maxpivot would technically not be needed for rational types and for interval types, because in neither case we need any kind of thresholding to determine if a pivot is exactly zero (rational), or is possibly zero (interval). But that's actually another problem to specialize rank(), and this likely again requires different specializations for all decompositions.

When I ask for the rank with intervals, I only mean it as a quick check whether it is obviously full rank, but I fully expect to need rationals for any other answer.

I think I didn't break maxPivot() in FullPivLU.

maxPivot wouldn't really make sense for rationals, because you are taking a 'best' pivot using some other criterion -- of course you can still calculate your maxPivot, but it has no significant meaning.

I just want to avoid breaking documented functions, even if they don't make much sense in this particular scenario.

householder, SVD, etc would likely need a similar treatment.

Of course, but let us do that after we agreed on how to finish this.


That said, maybe there is a need to let the user provide different pivoting strategies, depending on his needs? Unfortunately, that risks bloating the interface.

For now providing a different pivot comparison function seems good enough to me.

Is the patch roughly the right approach? Where do we go from there?

I would consider this the right approach especially as long as we don't change existing behavior. But I'd like to hear the opinion of other developers before going further.

Well, it is changing slightly the situation with abs: twice as many calls, no packet support, etc.

Of course, to make your patch worth the effort, we'd also need an example of a specialized scalar_better_coeff_op, demonstrating the benefits.

Oh, did I forget to say? I tested it on a branch of CGAL I am working on. All this test is doing with eigen is compute exactly the sign of the determinant of 9x9 matrices of double. First, it converts to intervals and computes that determinant. If it contains 0, it tries again with a rational type. I was surprised that the determinant computation was almost always failing (this number type throws an exception when we do (bool)(x<y) and it doesn't know the answer) and it ended up using rationals all the time. With the patch and the following, the interval version now succeeds all the time on random data.

--- a/Number_types/include/CGAL/Interval_nt.h
+++ b/Number_types/include/CGAL/Interval_nt.h
@@ -1446,6 +1446,18 @@ namespace Eigen {
       MulCost = 10
+  namespace internal {
+    template<class> struct scalar_better_coeff_op;
+    template<bool b>
+      struct scalar_better_coeff_op <CGAL::Interval_nt<b> > {
+	typedef CGAL::Interval_nt<b> Scalar;
+	typedef bool result_type;
+	bool operator()(Scalar const&x, Scalar const&y)const{
+	  return abs(x).inf()>abs(y).inf();
+	}
+      };
+  }


Marc Glisse

Mail converted by MHonArc 2.6.19+