Sparse Matrix (gpuArray) LU Decomposition

How can I do LU decomposition on GPU device to a sparse-matrix (H)? Suppose H is a sparse matrix, I did

>> H=gpuArray(H);
>> [L,U,P]=lu(H);

Error using gpuArray/lu
Sparse gpuArrays are not supported for this function.

Q1: Could you please let me know which function allows me to apply LU decomposition to a sparse matrix on GPU?

Later, I will use the L, U and P to process a series of data (d1,d2, ... dn) individually as below. Each data is a matrix where each column is d.

>>for id=1:nd

>>result(:,id)=L\U\P*data(:,id);

>>end

In this case, I can recycle the L, U and P in each inversion for a given data(:,id).

Thank you.

[Matt J]

What do you plan to do with the decomposition once you have it?

[Ted]

Later, I will use the L, U and P to process a series of data (d1,d2, ... dn) individually as below. Each data is a matrix where each column is d.

>>for id=1:nd

>>result(:,id)=L\U\P*data(:,id);

>>end

In this case, I can recycle the L, U and P in each inversion for a given data(:,id).

Thank you.

[Bruno Luong]

On sparse matrix LU use some graph technique to limits the band of the matrix so that avoid to fill L/U with too many non-zero elements. I suspect the graph method used by LU on sparse is not parallelizable and not suitable for GPU calculation.

Just stick with CPU solution for now.

Answer 1

[Joss Knight]

It looks like you're just after the backslash operator, so why not use that, and forget about looping over the RHS.

result = gpuArray(H) \ data;