[GiNaC-list] matrix::solve()-related problems
Vitaly Magerya
vmagerya at gmail.com
Mon Jan 22 17:09:37 CET 2018
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;
More information about the GiNaC-list
mailing list