From: Michael Orlitzky <michael@orlitzky.com>
Date: Sat, 21 Nov 2020 02:32:57 +0000 (-0500)
Subject: eja: echelonize a matrix manually to ensure a decent algorithm is used.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=e78216245caff4f11de66433f28fd995f4670b78;p=sage.d.git

eja: echelonize a matrix manually to ensure a decent algorithm is used.
---

diff --git a/mjo/eja/eja_element_subalgebra.py b/mjo/eja/eja_element_subalgebra.py
index a26381b..d2b54fe 100644
--- a/mjo/eja/eja_element_subalgebra.py
+++ b/mjo/eja/eja_element_subalgebra.py
@@ -1,4 +1,5 @@
 from sage.matrix.constructor import matrix
+from sage.rings.all import QQ
 
 from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
 
@@ -18,6 +19,19 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         P = matrix(field, power_vectors)
 
         if orthonormalize_basis == False:
+            # Echelonize the matrix ourselves, because otherwise the
+            # call to P.pivot_rows() below can choose a non-optimal
+            # row-reduction algorithm. In particular, scaling can
+            # help over AA because it avoids the RecursionError that
+            # gets thrown when we have to look too hard for a root.
+            #
+            # Beware: QQ supports an entirely different set of "algorithm"
+            # keywords than do AA and RR.
+            algo = None
+            if field is not QQ:
+                algo = "scaled_partial_pivoting"
+            P.echelonize(algorithm=algo)
+
             # In this case, we just need to figure out which elements
             # of the "powers" list are redundant... First compute the
             # vector subspace spanned by the powers of the given