From 45ca93fc48c14f733de73ffbbfef0834be731b08 Mon Sep 17 00:00:00 2001
From: Vitaly Magerya <vmagerya@gmail.com>
Date: Wed, 20 Jun 2018 01:13:07 +0200
Subject: [PATCH] Consider solve_algo::markowitz in automatic elimination
 algorithm selection.

Cf. <https://www.ginac.de/pipermail/ginac-list/2018-June/002211.html>.
---
 ginac/flags.h    |  8 +++++---
 ginac/matrix.cpp | 37 +++++++++++++++++++++++++++----------
 2 files changed, 32 insertions(+), 13 deletions(-)

diff --git a/ginac/flags.h b/ginac/flags.h
index cebc3fd4..f759dcd4 100644
--- 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 6464ed8d..f5c99a37 100644
--- 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);
-- 
2.44.0