diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 915d4543..0c3a305a 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -1014,21 +1014,23 @@ matrix matrix::solve(const matrix & vars,
 	
 	// Gather some statistical information about the augmented matrix:
 	bool numeric_flag = true;
+	unsigned sparse_count = 0;  // counts non-zero elements
 	for (auto & r : aug.m) {
-		if (!r.info(info_flags::numeric)) {
+		if (!r.info(info_flags::numeric))
 			numeric_flag = false;
-			break;
-		}
+		if (r.is_zero())
+			++sparse_count;
 	}
 	
 	// Here is the heuristics in case this routine has to decide:
 	if (algo == solve_algo::automatic) {
-		// Bareiss (fraction-free) elimination is generally a good guess:
-		algo = solve_algo::bareiss;
-		// For m<3, Bareiss elimination is equivalent to division free
-		// elimination but has more logistic overhead
-		if (m<3)
-			algo = solve_algo::divfree;
+		// Gauss elimination is generally a good guess:
+		algo = solve_algo::gauss;
+		// For dense matrices Bareiss elimination is faster than Gauss,
+		// and for m<3 division-free elimination is equivalent to Bareiss
+		// but has more logistic overhead.
+		if (5*sparse_count<=row*col)
+			algo = (m<3) ? solve_algo::divfree : solve_algo::bareiss;
 		// This overrides any prior decisions.
 		if (numeric_flag)
 			algo = solve_algo::gauss;