diff --git a/src/LinearAlgebra/Reduction/RowEchelonForm.php b/src/LinearAlgebra/Reduction/RowEchelonForm.php index 510c6d7d7..f2be68c8d 100644 --- a/src/LinearAlgebra/Reduction/RowEchelonForm.php +++ b/src/LinearAlgebra/Reduction/RowEchelonForm.php @@ -127,6 +127,19 @@ public static function gaussianElimination(NumericMatrix $A): array $swaps = 0; $ε = $A->getError(); + // Scale the matrix by its smallest non-zero element + $min = PHP_INT_MAX; + for ($i = 0; $i < $m; $i++) { + for ($j = 0; $j < $n; $j++) { + $elem = \abs($A[$i][$j]); + if (Support::isNotZero($elem, $ε) and $elem < $min) { + $min = $elem; + } + } + } + + $R = $A->scalarMultiply(1/$min)->getMatrix(); + for ($k = 0; $k < $size; $k++) { // Find column max $i_max = $k; diff --git a/tests/LinearAlgebra/Reduction/ReducedRowEchelonFormTest.php b/tests/LinearAlgebra/Reduction/ReducedRowEchelonFormTest.php index bdb426758..240dba643 100644 --- a/tests/LinearAlgebra/Reduction/ReducedRowEchelonFormTest.php +++ b/tests/LinearAlgebra/Reduction/ReducedRowEchelonFormTest.php @@ -1045,6 +1045,44 @@ public function dataProviderForRref(): array [0, 0, 0, 0, 0], ], ], + // Large magnitude case + [ + [ + [-2.82878905765*10**6, 0, 0, 0, 0, 0], + [0, -2.82878905765*10**6, 0, 0, 0, 1729.7], + [0, 0, -2.82878905765*10**6, 0, -1729.7, 0], + [0, 0, 0, -2.82879105765*10**6, 0, 0], + [0, 0, -1729.7, 0, -1.05765, 0], + [0, 1729.7, 0, 0, 0, -1.05765] + ], + [ + [1, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 1] + ] + ], + // Scaled version of large magnitude case + [ + [ + [-0.99999966687087, 0, 0, 0, 0, 0], + [0, -0.99999966687087, 0, 0, 0, 0.00061146285160793], + [0, 0, -0.99999966687087, 0, -0.00061146285160793, 0], + [0, 0, 0, -1.0000003738869, 0, 0], + [0, 0, -0.00061146285160793, 0, -3.7388694340557*10**-7, 0], + [0, 0.00061146285160793, 0, 0, 0, -3.7388694340557*10**-7] + ], + [ + [1, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 1] + ] + ] ]; }