[Octave-bug-tracker] [bug #32364] problem multiplying permutation and sparse matrices

[John W. Eaton]

Follow-up Comment #1, bug #32364 (project octave):

[Reposting because I screwed up a verbatim tag and the report on the web was
incomplete; I wish this thing had a preview button...]

While investigating the following bug report

https://mailman.cae.wisc.edu/pipermail/bug-octave/2011-February/017300.html

I uncovered a problem with multiplying permutation and sparse matrices.

Given


A = sparse ([1   1   1   0   0   0   0   0   0
             0   0   0   1   1   1   0   0   0
             0   0   0   0   0   0   1   1   1
             1   0   0   1   0   0   1   0   0
             0   1   0   0   1   0   0   1   0
             0   0   1   0   0   1   0   0   1])
[L, U, P, Q, R] = lu (A);


I was trying to verify that


P*(RA)*Q - L*U


was numerically zero, but found


octave:8> P*(RA)*Q - L*U
ans =

Compressed Column Sparse (rows = 6, cols = 9, nnz = 4 [7.4%])

  (4, 5) -> -0.33333
  (4, 5) ->  0.33333
  (3, 7) -> -0.33333
  (3, 7) ->  0.33333


Note the strange result that some elements apparently have two values!

However,


full (P)*(RA)*Q - L*U


returns a matrix of zeros.  Playing around a bit more, I noticed 


octave:10> P*A
ans =

Compressed Column Sparse (rows = 6, cols = 9, nnz = 18 [33%])

  (1, 1) ->  1
  (3, 1) ->  1
  (1, 2) ->  1
  (4, 2) ->  1
  (1, 3) ->  1
  (6, 3) ->  1
  (2, 4) ->  1
  (3, 4) ->  1
  (2, 5) ->  1
  (4, 5) ->  1
  (2, 6) ->  1
  (6, 6) ->  1
  (5, 7) ->  1
  (3, 7) ->  1
  (5, 8) ->  1
  (4, 8) ->  1
  (5, 9) ->  1
  (6, 9) ->  1


Note the elements in columns 7 and 8 are out of order.  I assume that the
sparse matrix multiplication algorithm assumes that the elements in each
column are ordered by ascending row number.

I've looked at the code in Sparse-perm-op-defs.h, but I think it will take me
a while to figure this out, so I would really appreciate it if someone who
understands this code better could take a look at this problem.

I'm adding the following people to the CC list for this report: Jason Reidy
since he is listed as the author of Sparse-perm-op-defs.h, Jaroslav Hajek
because he is the author of the permutation matrix code, and David Bateman
because he has good knowledge of the sparse matrix code.


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?32364>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/