From 7c3c5284314b4c0edc4d23ca3ee4db11aa0839c9 Mon Sep 17 00:00:00 2001
From: Vitaly Magerya <magv@tx97.net>
Date: Thu, 14 Jun 2018 14:31:28 +0200
Subject: [PATCH] Improve automatic algorithm selection in matrix::solve()

---
 ginac/matrix.cpp | 35 +++++++++++++++++++++++++----------
 1 file changed, 25 insertions(+), 10 deletions(-)

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 3a01da2d..1bb33e0c 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -1211,21 +1211,36 @@ 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) {
+		unsigned ncells = col*row;
+		unsigned density = 0;
+		for (auto && r : m) {
 			if (!r.info(info_flags::numeric)) {
 				numeric_flag = false;
 				break;
 			}
+			density += !r.is_zero();
+		}
+		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;
+			}
 		}
-		// 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;
 	}
 	// Eliminate the augmented matrix:
 	std::vector<unsigned> colid(col);
-- 
2.13.6