[GiNaC-list] matrix::solve()-related problems

[Vitaly Magerya]

On 01/18/2018 03:10 PM, I wrote:
> [...]
> Now, there are multiple ways of dealing with these problems. The
> practical workaround for the moment is to:
> a) use solve_algo::gauss for all sparse-ish matrices
> b) normal()ize the whole matrix before passing it into solve()
> 
> There are places where this workaround is not available. One of
> such places is matrix::inverse()

The second place where a better algorithm is needed is matrix::rank().
Currently it blindly fraction_free_elimination(), which should never be
used with sparse matrices.

BTW, I noticed that matrix::determinant() has some logic to determine
matrix sparseness, and it would actually switch away from
fraction_free_elimination() for sparse matrices (it will use minor
expansion method instead), where sparseness is determined by comparing
the number of zero cells to the number of all cells divided by 5.

I think it would be beneficial if the same sort of logic would be used
in matrix::solve() too. (See the attached patch for my proposal). This
would additionally alleviate the need for an 'algo' argument in
matrix::inverse(), since solve_algo::automatic would now perform the
sensible thing for sparse matrices (which is Gauss elimination).

This of course doesn't solve the problem with matrix::rank().

I think there are two good ways to proceed:
1) Apply the attached patch, and make matrix::rank() use similar logic
to determine which *_elimination() to call.
2) Apply the attached patch, and add 'algo' arguments to both
matrix::inverse() and matrix::rank() (probably using separate
enumerations for each).
-------------- next part --------------
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;