rref returns wrong results for sparse matrices

[David Bateman]

John W. Eaton wrote:
> 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.
>   

The above is in column major order. Doing

 A = sparse ([1, 1.5, 2, 2.5; 1, 2, 3, 4; 2, 3, 4, 2]);
 A([2,3],1:4) = A([2,3],1:4) - A([2,3],1) * A(1,1:4);
 BS = full(A);
 A = [1, 1.5, 2, 2.5; 1, 2, 3, 4; 2, 3, 4, 2];
 A([2,3],1:4) = A([2,3],1:4) - A([2,3],1) * A(1,1:4);
 BF = A;
 BS == BF

I see,

ans =

   1   1   1   1
   1   1   1   1
   1   1   1   1

so I don't see the issue with full/sparse matrices here.
> I also see nnz being greater than the number of elements in the
> matrix, so that doesn't seem right.
>   
Given the assumptions I made in the code nnz should not be bigger than
the number of elements, though isn't necessarily an issue as such. I'd
consider it a bug at this point unless we implement a real spalloc function.

> I'm not all that familiar with the sparse matrix code in Octave, so I
> think someone else will have to fix this.
>   
I'll take a look. Just for reference the problem is in the sparse assign
functions (1 and 2 index value versions) when there are the same values
in the index. IE A([1,1,1]) = ..., A([1,1],;0 = ..., etc


D.

> jwe
> _______________________________________________
> Bug-octave mailing list
> Bug-octave at octave.org
> https://www.cae.wisc.edu/mailman/listinfo/bug-octave
>
>   


-- 
David Bateman                                David.Bateman at motorola.com
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary