rref returns wrong results for sparse matrices

[John W. Eaton]

On  2-Mar-2008, Nils Bluthgen wrote:

| Bug report for Octave 3.0.0 configured for i386-apple-darwin8.9.1
| 
| Description:
| -----------
| 
| Calculation of the reduced row echelon form with rref.m gives  
| different results for sparse and full matrices.
| 
| Repeat-By:
| ---------
|   A=[1 2 3 4; 2 3 4 5; 2 3 4 2]
| rref(A)
| full(rref(sparse(A)))

In watching what rref is doing with sparse matrices, I noticed this
strange behaivor:

  octave:11> A = sparse ([1, 1.5, 2, 2.5; 1, 2, 3, 4; 2, 3, 4, 2])
  A =

  Compressed Column Sparse (rows = 3, cols = 4, nnz = 12)

    (1, 1) ->  1
    (2, 1) ->  1
    (3, 1) ->  2
    (1, 2) ->  1.5000
    (2, 2) ->  2
    (3, 2) ->  3
    (1, 3) ->  2
    (2, 3) ->  3
    (3, 3) ->  4
    (1, 4) ->  2.5000
    (2, 4) ->  4
    (3, 4) ->  2

  octave:12> A([2,3],1:4) = A([2,3],1:4) - A([2,3],1) * A(1,1:4)
  A =

  Compressed Column Sparse (rows = 3, cols = 4, nnz = 8)

    (1, 1) ->  1
    (1, 2) ->  1.5000
    (2, 2) ->  0.50000
    (1, 3) ->  2
    (2, 3) ->  1
    (1, 4) ->  2.5000
    (2, 4) ->  1.5000
    (3, 4) -> -3


Is that normal?  I would expect the output to always be column-major
ordering.

I also see nnz being greater than the number of elements in the
matrix, so that doesn't seem right.

I'm not all that familiar with the sparse matrix code in Octave, so I
think someone else will have to fix this.

jwe