Talk:Reduced row echelon form
"Break" vs. "return" bug
The original author of the Python example mistakenly translated the keyword stop that appears in the Wikipedia pseudocode as break rather than the correct return. This created a control-flow bug that didn't manifest itself when the program was run on the example matrix given in the task description, but did cause an exception if the program was run on, e.g.,
1 2 3 4 3 1 2 4 6 2 6 2 3 6 18 9 9 -6 4 8 12 10 12 4 5 10 24 11 15 -4
I noticed and fixed the bug a couple of days ago, but it seems that several of the other examples written before then (being, by and large, translations from the Python) copied the bug. Hence, I've marked all the examples that looked like they might have this bug with the needs-review template. Note that I erred on the side of false positives. That is, I'm pretty sure the examples I didn't mark are bug-free, but some of the ones I did mark may be fine. —Underscore 00:25, 5 May 2009 (UTC)
I ran the ALGOL 68 version and got:
(( 1.0000, 2.0000, 0.0000, 0.0000, 3.0000, 4.0000), ( 0.0000, 0.0000, 1.0000, 0.0000, 0.0000, -1.0000), ( 0.0000, 0.0000, 0.0000, 1.0000, 0.0000, 0.0000), ( 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000), ( 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000))
Vs the current python version which got:
1, 2, 0, 0, 3, 4 0, 0, 1, 0, 0, -1 0, 0, 0, 1, 0, 0 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0
It looks like - currently- they are both getting the same answer... next step is it to get out a pencil.
NevilleDNZ 01:46, 5 May 2009 (UTC)
Bug in Common Lisp code
There is a bug in find-pivot that causes an infinite loop on some input. Here is a corrected version along with the matrix that results in an infinite loop with the current code. I tested the code with SBCL 1.0.40.0.debian on Ubuntu 10.10 64bit.
<lang lisp>(defun convert-to-row-echelon-form (matrix) (defun convert-to-row-echelon-form (matrix)
(let* ((dimensions (array-dimensions matrix))
(row-count (first dimensions)) (column-count (second dimensions)) (lead 0))
(labels ((find-pivot (start lead)
(let ((i start)) (loop :while (zerop (aref matrix i lead)) :do (progn (incf i) (when (= i row-count) (setf i start) (incf lead) (when (= lead column-count) (return-from convert-to-row-echelon-form matrix)))) :finally (return (values i lead))))) (swap-rows (r1 r2) (loop :for c :upfrom 0 :below column-count :do (rotatef (aref matrix r1 c) (aref matrix r2 c)))) (divide-row (r value) (loop :for c :upfrom 0 :below column-count :do (setf (aref matrix r c) (/ (aref matrix r c) value)))))
(loop
:for r :upfrom 0 :below row-count :when (<= column-count lead) :do (return matrix) :do (multiple-value-bind (i nlead) (find-pivot r lead) (setf lead nlead) (swap-rows i r) (divide-row r (aref matrix r lead)) (loop :for i :upfrom 0 :below row-count :when (/= i r) :do (let ((scale (aref matrix i lead))) (loop :for c :upfrom 0 :below column-count :do (decf (aref matrix i c) (* scale (aref matrix r c)))))) (incf lead)) :finally (return matrix)))))
(defvar *M* '(( 1 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 0)
( 1 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0) ( 1 0 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0) ( 0 1 0 0 0 0 1 0 0 0 0 0 0 0 -1 0 0 0) ( 0 1 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0) ( 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 -1 0) ( 0 0 1 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0) ( 0 0 1 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0) ( 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 -1 0 0) ( 0 0 0 1 0 0 0 0 0 1 0 0 -1 0 0 0 0 0) ( 0 0 0 0 1 0 0 1 0 0 0 0 0 -1 0 0 0 0) ( 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 -1 0) ( 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 -1 0 0) ( 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0) ( 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0) ( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1) ( 0 0 0 0 0 1 0 0 0 0 1 0 0 0 -1 0 0 0)))
(defvar *M-array* (make-array (list (length *M*) (length (first *M*))) :initial-contents *M*))
(print (convert-to-row-echelon-form *M-array*)) </lang>
Carlo Hamalainen 3:30PM, 6 Dec 2010 (GMT+10)
Fixed on main page.
Carlo Hamalainen 9:39AM, 16 Mar 2011 (GMT+10)
Code doesn't work for singular matrices
I tried putting in a singular matrix into the C# code and it is giving divide by zero errors. When I told it to skip when it is dividing by zero, the results keep changing. Also I corrected an error in the C# code that caused the lead value to be out of bounds.
- I won't be able to test it myself, but can you provide a sample matrix that reproduces the problem? --Michael Mol 03:08, 10 October 2010 (UTC)
- this gives a divide by zero:
- 1,0,1,0,1,0
- 1,0,1,0,0,1
- 1,0,0,1,1,0
- 1,0,0,1,0,1
- 0,1,0,1,1,0
- 0,1,0,1,0,1
- 0,1,1,0,1,0
- 0,1,1,0,0,1
- I haven't looked into it in detail, but adding an "if (div != 0)" line before the offending division seems to make it behave reasonably.
- edit: looking into it, the lead-- line seems to deviate from the original algorithm, replacing that with a "return matrix;" seems to fix it.
New test cases
Recently, I found that the earlier Java version did not work as expected due to floating point errors.
I recommend including these in your respecting languages: double matrix_2 [][] = { {2, 0, -1, 0, 0}, {1, 0, 0, -1, 0}, {3, 0, 0, -2, -1}, {0, 1, 0, 0, -2}, {0, 1, -1, 0, 0} };
solution: [[1, 0, 0, 0, -1] [0, 1, 0, 0, -2] [0, 0, 1, 0, -2] [0, 0, 0, 1, -1] [0, 0, 0, 0, 0]]
double matrix_3 [][] = { {1, 2, 3, 4, 3, 1}, {2, 4, 6, 2, 6, 2}, {3, 6, 18, 9, 9, -6}, {4, 8, 12, 10, 12, 4}, {5, 10, 24, 11, 15, -4} };
solution: [[1, 2, 0, 0, 3, 4] [0, 0, 1, 0, 0, -1] [0, 0, 0, 1, 0, 0] [0, 0, 0, 0, 0, 0] [0, 0, 0, 0, 0, 0]]
double matrix_4 [][] = { {0, 1}, {1, 2}, {0,5} };
solution: [[1, 0] [0, 1] [0, 0]]