From 53ed991b5deef876e9a5757a3595b36e0052daac Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Thu, 27 Aug 2020 15:39:08 +0200
Subject: [PATCH] Fall back to LU instead of LDLT with old Eigen

---
 .../include/CGAL/NewKernel_d/Types/Weighted_point.h | 13 ++++++++++++-
 .../CGAL/NewKernel_d/function_objects_cartesian.h   | 13 ++++++++++++-
 2 files changed, 24 insertions(+), 2 deletions(-)

diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h
index e83af0db167..6e970a5f13a 100644
--- a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h
@@ -236,13 +236,24 @@ template <class R_> struct Power_center : Store_kernel<R_> {
       }
       // Only need to fill the lower half
       for(int i = 0; i < k-1; ++i){
-        for(int j = i; j < k-1; ++j)
+        for(int j = i; j < k-1; ++j){
           m(j, i) = sp(vecs[i], vecs[j]);
+#if ! EIGEN_VERSION_AT_LEAST(3, 3, 5)
+          m(i, j) = m(j, i);
+#endif
+        }
         b[i] += m(i, i);
         b[i] /= 2;
       }
       // Assumes Eigen...
+#if EIGEN_VERSION_AT_LEAST(3, 3, 5)
       Vec res = m.ldlt().solve(b);
+#else
+      // Older versions of Eigen use 1/highest as tolerance,
+      // which we have no way to set to 0 for exact types.
+      // Use something slow but that should work.
+      Vec res = m.fullPivLu().solve(b);
+#endif
       Vector to_center = sv(vecs[0], res[0]);
       for(int i=1;i<k-1;++i)
         to_center = pv(to_center, sv(vecs[i],res[i]));
diff --git a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h
index 3ca5fb3638e..3ab117131dc 100644
--- a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h
@@ -693,12 +693,23 @@ template <class R_> struct Construct_circumcenter : Store_kernel<R_> {
         vecs.emplace_back(dp(*f,p0));
       // Only need to fill the lower half
       for(int i=0;i<k-1;++i){
-        for(int j=i;j<k-1;++j)
+        for(int j=i;j<k-1;++j) {
           m(j,i)=sp(vecs[i],vecs[j]);
+#if ! EIGEN_VERSION_AT_LEAST(3, 3, 5)
+          m(i,j)=m(j,i);
+#endif
+        }
         b[i]=m(i,i)/2;
       }
       // Assumes Eigen...
+#if EIGEN_VERSION_AT_LEAST(3, 3, 5)
       Vec res=m.ldlt().solve(b);
+#else
+      // Older versions of Eigen use 1/highest as tolerance,
+      // which we have no way to set to 0 for exact types.
+      // Use something slow but that should work.
+      Vec res=m.fullPivLu().solve(b);
+#endif
       Point center=p0;
       // Wasteful if we only want the radius
       for(int i=0;i<k-1;++i)