## Seeking advice on choice of sparse matrix format

Open discussion for MAGMA library (Matrix Algebra on GPU and Multicore Architectures)
fletchjp
Posts: 203
Joined: Mon Dec 27, 2010 7:29 pm

### Seeking advice on choice of sparse matrix format

I have two complex matrices A and B which are generated in memory in setting up a problem.

For a case about 4000 square, A averages 850 elements per row and B has 16 elements per row.

At present I store all the elements of both A and B and for B this is clearly wasteful.

I need to do two operations involving B

1. Add a multiple of it to A
2. Do repeated multiplications Bx with different vectors x.

I want to move to using the sparsity of B, for the moment tackling operation 2 by storing B in a sparse form and using a sparse multiplication. I need to chose a sparse format for B and also find a sparse BLAS which will do the multiplication. I want to start by doing this on the CPU, but have the possibility to move to using Magma sparse routines. I am working in Fortran 90.

I am currently digging into how Magma sparse matrices and routines are defined. I would welcome some advice on the choice of sparse matrix format and sparse BLAS. It would also be useful if this included a routine for adding a sparse matrix to a dense matrix.

Thank you.

John

hartwig anzt
Posts: 90
Joined: Tue Sep 02, 2014 5:44 pm

### Re: Seeking advice on choice of sparse matrix format

Dear John,

a 4000x4000 matrix is still relatively small. Sparse iterative systems typically target problems with more than 100,000 unknowns. So it might be an option to handle this system as dense. However, if you want, you can handle it as sparse. The easiest way is to convert it to csr, and then use the MAGMA-sparse functions for addition and multiplication.

Assuming the elements of your dense matrix is stored in M, you should be able to use the following code:

// pass it to MAGMA-sparse:
// create a MAGMA_Z_MATRIX

magma_z_matrix A, A_CSR, A_dev;
A.num_rows = 4000;
A.num_cols = 4000;
A.storage_type = Magma_DENSE;
A.val = M;

// convert to CSR
magma_zmconvert( A, &A_CSR, Magma_DENSE, Magma_CSR, queue );

// copy to device
magma_zmtransfer( A_CSR, A_dev, Magma_CPU, Magma_DEV, queue );

at this point you are able to throw all linear algebra tools at it, e.g.
magma_z_matrix A,
magma_z_matrix B,
magma_z_matrix * AB,
magma_queue_t queue
)

* sparse matrix sparsematrix multiplication:
magma_z_matrix x,
magma_z_matrix y,
magma_queue_t queue
)

* sparse matrix vector product
magma_z_spmv (magmaDoubleComplex alpha, magma_z_matrix A, magma_z_matrix x, magmaDoubleComplex beta, magma_z_matrix y, magma_queue_t queue)

Thanks, Hartwig

fletchjp
Posts: 203
Joined: Mon Dec 27, 2010 7:29 pm

### Re: Seeking advice on choice of sparse matrix format

Hartwig

Thank you very much for the quick reply. Indeed, my 4000 by 4000 is just a small case, they get much bigger. I use this case for development. I don't see much advantage using Magma at that size, but it can then be used on larger cases.

I think you have dealt with most of the things. There is an argument called queue on some of the calls. Do I have to set that to something or does it have a convenient default?

My main problem with Magma sparse at the moment is that I cannot run with CUDA until I get CUDA 6.5 running and for that I need to update my NVIDIA drivers. I am still figuring out how to do that without crashing my computer, which I have done once already.

I am hoping I can at least use magma_zmconvert to move the matrix around. Is there a FORTRAN wrapper for it?

I then need to find a CPU sparse BLAS to work on it.

Thanks again

John

hartwig anzt
Posts: 90
Joined: Tue Sep 02, 2014 5:44 pm

### Re: Seeking advice on choice of sparse matrix format

John,

are you working on your personal Desktop computer or do you have access to some cluster? In the second case, I could probably help getting the code running - if granted access. For the CPU side, Intel's MKL provides some sparse support. Unfortunately, I have to admit that the FORTRAN support is very limited... Is it a complex code or is switching to C an option?

Thanks, Hartwig

fletchjp
Posts: 203
Joined: Mon Dec 27, 2010 7:29 pm

Hartwig