Inconsistent Results with cusp::multiply using coo_matrix

[GoogleCodeExporter]

What steps will reproduce the problem?
Using this code:

#include <cusp/multiply.h>
#include <cusp/print.h>

#define MAX_SIZE 8
#define ELEMENTS 2 * MAX_SIZE

int main(void)
{
  // initialize matrix
  cusp::coo_matrix<int, double, cusp::host_memory> A_host(MAX_SIZE,MAX_SIZE,ELEMENTS);
  for(uint32_t j = 0; j < ELEMENTS; j++) {
    A_host.row_indices[j] = j % MAX_SIZE;
    A_host.column_indices[j] = j / MAX_SIZE;
    A_host.values[j] = 1;
  }
  cusp::coo_matrix<int, double, cusp::device_memory> A = A_host;
  cusp::print(A_host);
  cusp::print(A);

  // allocate output vector
  cusp::coo_matrix<int, double, cusp::device_memory>  y      (MAX_SIZE, MAX_SIZE, MAX_SIZE*MAX_SIZE);
  cusp::coo_matrix<int, double, cusp::host_memory>    y_host (MAX_SIZE, MAX_SIZE, MAX_SIZE*MAX_SIZE);

  // compute y = A * A
  cusp::multiply(A, A, y);
  cusp::multiply(A_host, A_host, y_host);

  // print y
  cusp::print(y);
  cusp::print(y_host);

  return 0;
}

What is the expected output? What do you see instead?
Should create an 8x8 matrix and fill the first two columns with 1's.  Result of 
multiply should be 8x8 matrix with first two columns filled with 2's.  Results 
do not seem correct and are not consistent across the CPU / GPU.

Result from Code:
sparse matrix <8, 8> with 16 entries
              0              0              1
              1              0              1
              2              0              1
              3              0              1
              4              0              1
              5              0              1
              6              0              1
              7              0              1
              0              1              1
              1              1              1
              2              1              1
              3              1              1
              4              1              1
              5              1              1
              6              1              1
              7              1              1
sparse matrix <8, 8> with 16 entries
              0              0              1
              1              0              1
              2              0              1
              3              0              1
              4              0              1
              5              0              1
              6              0              1
              7              0              1
              0              1              1
              1              1              1
              2              1              1
              3              1              1
              4              1              1
              5              1              1
              6              1              1
              7              1              1
sparse matrix <8, 8> with 16 entries
              0              0              8
              0              1              2
              1              0              8
              1              1              2
              2              0              8
              2              1              2
              3              0              8
              3              1              2
              4              0              8
              4              1              2
              5              0              8
              5              1              2
              6              0              8
              6              1              2
              7              0              8
              7              1              2
sparse matrix <8, 8> with 9 entries
              0              1              9
              0              0             64
              1              1              1
              2              1              1
              3              1              1
              4              1              1
              5              1              1
              6              1              1
              7              1              1

What version of the product are you using? On what operating system?
Using CUSP v0.3.0 and CUDA v4.1 on Ubuntu 11.10 with a GTX 560.

NOTE: It is quite possible that I am using the API incorrectly.  If so, sorry 
for wasting your time - feel free to mock me :-)

Original issue reported on code.google.com by robert.c...@gmail.com on 2 Apr 2012 at 10:31

[GoogleCodeExporter]

Hi,

The problem with the example is that the row indices have to be ordered, for 
example:

 for(uint32_t j = 0; j < ELEMENTS; j++) {
   A_host.row_indices[j] = j / 2;
   A_host.column_indices[j] = j % 2;
   A_host.values[j] = 1;
 }

This will give the correct result.
Maybe it would be good to have a way to automatically check for the correctness 
of COO matrices as this problem has happened several times.

Cheers,
Filipe

Original comment by filipe.c...@gmail.com on 2 Apr 2012 at 11:26

[GoogleCodeExporter]

Checking format correctness has non-trivial cost, so we couldn't have it 
enabled all the time.  However, one thing we could do is apply 
assert_is_valid_matrix [1] to the appropriate arguments when some CUSP_DEBUG 
compilation mode was enabled.  The format verification in cusp/verify.h is 
fairly rigorous [2], so it could lessen the number of such reports.

[1] http://code.google.com/p/cusp-library/source/browse/cusp/verify.h
[2] 
http://code.google.com/p/cusp-library/source/browse/cusp/detail/verify.inl#96

Original comment by wnbell on 3 Apr 2012 at 1:44

[GoogleCodeExporter]

Aaaah... sorry about that! Right there in the doxygen too: "The matrix entries 
must be sorted by row index."  Thanks for the help.

Original comment by robert.c...@gmail.com on 3 Apr 2012 at 1:56