How to interface Cusp's solvers with ell_matrix_view?

[Tao Han]

Hi,

I'm a new user of Cusp. I plan to solve a high dimensional sparse linear system (230400 * 230400 matrix A with 11 nonzeros per row) based on ell_matrix_view format. But I am troubled by how to interface Cusp's solvers with ell_matrix_view.

To make the problem simple, I wrote a simple example for testing. I construct a 4X4 using the data stored in device memory:

typedef cusp::ell_matrix_view<DeviceArray2dViewInt, DeviceArray2dViewFloat> EllMatView;

EllMatView A(4, 4, 10, col_indices_array2d, mat_values_array2d);

But when I use the following solver, the residual norm always diverge:

cusp::krylov::cg(A, x, b, monitor);

......

Iteration Number  | Residual Norm

                0       3.162278e+00

                1              inf

                2              nan

                3              nan

                4              nan

                5              nan

                6              nan

However, if I save the ell_matrix_view on disk and reload it as an ell_matrix, then everything works well:

cusp::io::write_matrix_market_file(A, "A.mtx");

cusp::ell_matrix<int, float, cusp::device_memory> A_ell(4, 4, 10, 3);

cusp::io::read_matrix_market_file(A_ell, "A.mtx");

.......

cusp::krylov::cg(A_ell, x, b, monitor);

......

  Iteration Number  | Residual Norm

                0       3.162278e+00

                1       1.581139e+00

                2       2.356080e-08

And I also print the ell_matrix_view and the ell_matrix. The results are identical. So I am a little confused here.

The source code of the example can be found in the attachment.

Thanks!

Tao

[Steven Dalton]

Hello,

  I think there are a few things going on here. All functions assume ELL indices and values are stored in column major format but you are creating a few indicating row_major which should either be a compile error or it should adapt the internal access patterns [1] . Another issue is that the alignment of entries for ELL should be propagated throughout but those settings seem to be ignored at the moment when a ELL matrix or view is created. I created an issue to track this bug since it could cause a lot of issues for people creating their own views.

You can see the assumptions of the internal methods if you try the following:

#include <cusp/coo_matrix.h>

#include <cusp/ell_matrix.h>

#include <cusp/print.h>

int main (int argc, char* argv[])

{

  cusp::array2d<float,cusp::host_memory> A(4,4);

  A(0,0) =  2; A(0,1) = -1; A(0,2) = 0; A(0,3) = 0;

  A(1,0) = -1; A(1,1) =  2; A(1,2) =-1; A(1,3) = 0;

  A(2,0) =  0; A(2,1) = -1; A(2,2) = 2; A(2,3) =-1;

  A(3,0) =  0; A(3,1) =  0; A(3,2) =-1; A(3,3) = 2;

  cusp::ell_matrix<int,float,cusp::host_memory> B(4,4,10,4);

  B = A;

  cusp::print(B.column_indices.values);

        return 0;

}

[1] https://github.com/cusplibrary/cusplibrary/blob/develop/cusp/ell_matrix.h#L132

[2] https://github.com/cusplibrary/cusplibrary/issues/74

--
You received this message because you are subscribed to the Google Groups "cusp-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cusp-users+unsubscribe@googlegroups.com.
To post to this group, send email to cusp-...@googlegroups.com.
Visit this group at https://groups.google.com/group/cusp-users.
For more options, visit https://groups.google.com/d/optout.

[Tao Han]

Hello,

Thank you for your suggestions. I have changed my code to make the ELL indices and values stored in column major. My simulation system works well now. :)

To make my system more effective, I have an additional problem:

In my code, I define a diagonal preconditioner and call the CUSP PCG solver like follow:

cusp::precond::diagonal<float, cusp::device_memory> M(A);

cusp::krylov::cg(A, x, b, monitor, M);

The result is good according to my experimental results. However, I notice that the matrix A of the linear system Ax = b is not only a sparse matrix but also a symmetric matrix in my simulation problem. So I wonder it is possible to use this feature to make the PCG solver converge more effectively? 

Thanks!

Tao

To unsubscribe from this group and stop receiving emails from it, send an email to cusp-users+...@googlegroups.com.