[GiNaC-list] [patch] Improve automatic algorithm selection in matrix::solve()

[Vitaly Magerya]

On 06/17/2018 11:57 PM, Richard B. Kreckel wrote:
> Two comments on your patch.
> 
> 1) This line
> +		for (auto && r : m) {
> looks unconventional due to it's using rvalue refs. Is this intentional?
> It doesn't make a difference in GCC (same asm code).

What I actually meant to write is "const auto &"; "auto &&" should
be the same, but shorter. I think the "&&" form was introduced
so that templates could seamlessly support both const and non-
const usage, but I found it useful in normal code too. In any
case, I'll make the "const" explicit (not that it matters in
this particular example).

> 2) Be careful when computing the matrix' density by counting non-zero
> elements: That loop breaks out on the first non-numeric element. The
> later decision on what yo use with density might be wrong. Maybe you
> want to compute numeric_flag and density in two separate loops?

Dammit. My mistake, sorry.

Here's the updated patch. I've additionally included a minor
update to the comment near solve_algo::markowitz definition.
-------------- next part --------------
From c4b8475e7b35891223009aa36610413164e7aac6 Mon Sep 17 00:00:00 2001
From: Vitaly Magerya <magv at tx97.net>
Date: Thu, 14 Jun 2018 14:31:28 +0200
Subject: [PATCH] Improve automatic algorithm selection for matrix::solve()

---
 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 3a01da2d..fea52c62 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.13.6