Consider solve_algo::markowitz in automatic elimination algorithm selection.

Cf. <https://www.ginac.de/pipermail/ginac-list/2018-June/002211.html>.

diff --git a/ginac/flags.h b/ginac/flags.h
index cebc3fd449966d647c6e745081f71bd286ca956b..f759dcd411bce1edc740a94cfbf1a2e10bc63ad5 100644 (file)
--- a/ginac/flags.h
+++ b/ginac/flags.h
@@ -185,9 +185,11 @@ public:
                bareiss,
                /** Markowitz-ordered Gaussian elimination. Same as the usual
                 *  Gaussian elimination, but with additional effort spent on
-                *  selection pivots that minimize fill-in. Much faster than the
-                *  methods above for large sparse matrices, marginally slower
-                *  than Gaussian elimination otherwise. */
+                *  selecting pivots that minimize fill-in. Faster than the
+                *  methods above for large sparse matrices (particularly with
+                *  symbolic coefficients), otherwise slightly slower than
+                *  Gaussian elimination.
+                */
                markowitz
        };
 };

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 6464ed8d52cd90731511ddf04e976d915fb734b9..f5c99a37b2d98fab176b97562e83c587ac3895b5 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -1211,21 +1211,38 @@ matrix::echelon_form(unsigned algo, int n)
        if (algo == solve_algo::automatic) {
                // Gather some statistical information about the augmented matrix:
                bool numeric_flag = true;
-               for (auto & r : m) {
+               for (const auto & r : m) {
                        if (!r.info(info_flags::numeric)) {
                                numeric_flag = false;
                                break;
                        }
                }
-               // Bareiss (fraction-free) elimination is generally a good guess:
-               algo = solve_algo::bareiss;
-               // For row<3, Bareiss elimination is equivalent to division free
-               // elimination but has more logistic overhead
-               if (row<3)
-                       algo = solve_algo::divfree;
-               // This overrides any prior decisions.
-               if (numeric_flag)
-                       algo = solve_algo::gauss;
+               unsigned density = 0;
+               for (const auto & r : m) {
+                       density += !r.is_zero();
+               }
+               unsigned ncells = col*row;
+               if (numeric_flag) {
+                       // For numerical matrices Gauss is good, but Markowitz becomes
+                       // better for large sparse matrices.
+                       if ((ncells > 200) && (density < ncells/2)) {
+                               algo = solve_algo::markowitz;
+                       } else {
+                               algo = solve_algo::gauss;
+                       }
+               } else {
+                       // For symbolic matrices Markowitz is good, but Bareiss/Divfree
+                       // is better for small and dense matrices.
+                       if ((ncells < 120) && (density*5 > ncells*3)) {
+                               if (ncells <= 12) {
+                                       algo = solve_algo::divfree;
+                               } else {
+                                       algo = solve_algo::bareiss;
+                               }
+                       } else {
+                               algo = solve_algo::markowitz;
+                       }
+               }
        }
        // Eliminate the augmented matrix:
        std::vector<unsigned> colid(col);